diff --git a/Project.toml b/Project.toml
index 44e6b30..a649907 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,6 +11,7 @@ ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StringDistances = "88034a9c-02f8-509d-84a9-84ec65e18404"
 
 [compat]
diff --git a/src/TestImages.jl b/src/TestImages.jl
index e7b5770..d6882b4 100644
--- a/src/TestImages.jl
+++ b/src/TestImages.jl
@@ -5,6 +5,7 @@ using StringDistances
 using ColorTypes
 using ColorTypes.FixedPointNumbers
 using NDTools: reorient
+using StaticArrays
 const artifacts_toml = abspath(joinpath(@__DIR__, "..", "Artifacts.toml"))
 
 export testimage, testimage_dip3e
@@ -243,40 +244,42 @@ function shepp_logan(N::Int, M::Int, O::Int; high_contrast::Bool=true, highContr
         # a faster case when ϕ == 0.0
         abs2(dx * b) + abs2(dy * a) < (a * b)^2
     end
-    
-    function _ellipse3d(dx, dy, dz, a, b, c, R)
 
-        # Apply rotation to the point 
-        tx, ty, tz = R * [dx, dy, dz]
+    # multiplying the transformation matrix
+    @inline function mat_mul(t::SVector{N, T}, matrix_c::SMatrix{N,N,T})::SVector{N,T} where {N,T}
+        return matrix_c * t
+    end
 
-        # 3D ellipse equation 
-        (abs2(tx / a) + abs2(ty / b) + abs2(tz / c)) <= 1.0
+    @inline function mat_div(t::SVector{N, T}, v::SVector{N, T})::SVector{N,T} where {N,T}
+        return t ./ v
+    end
+    @inline function sum_abs2(t::SVector{N, T})::T where {N,T}
+        return sum(abs2.(t))
     end
 
     P = zeros(Float64, N, M, O)
     for l = 1:length(ϕ)
         if ϕ[l] == 0.0 && θ[l] == 0.0 && ψ[l] == 0.0
             # Rotation matrix using Z-Y-Z Euler angles
-            R = [
+            R = SMatrix{3, 3, Float64}([
                 1.0  0.0  0.0;
                 0.0  1.0  0.0;
                 0.0  0.0  1.0 
-            ]
-            P .+= A[l] .* _ellipse3d.(x .- x₀[l], y .- y₀[l], z .- z₀[l], a[l], b[l], c[l], Ref(R))
+            ])
         else
             sinϕ, cosϕ = sincosd(ϕ[l])
             sinθ, cosθ = sincosd(θ[l])
             sinψ, cosψ = sincosd(ψ[l])
             # Rotation matrix using Z-Y-Z Euler angles
-            R = [
+            R = SMatrix{3, 3, Float64}([
                 (cosψ*cosϕ-cosθ*sinψ*sinϕ)  (cosψ*sinϕ+cosθ*sinψ*cosϕ)  (sinψ*sinθ);
                 (-sinψ*cosϕ-cosθ*cosψ*sinϕ) (-sinψ*sinϕ+cosθ*cosψ*cosϕ) (cosψ*sinθ);
                 (sinθ*sinϕ)                     (-sinθ*cosϕ)                   (cosθ) 
-            ]
-            P .+= A[l] .* _ellipse3d.(x .- x₀[l], y .- y₀[l], z .- z₀[l], a[l], b[l], c[l], Ref(R))
+            ])
         end
-    end
+        P .+= A[l] .* (sum_abs2.(mat_div.(mat_mul.(SVector{3, Float64}.(x .- x₀[l], y .- y₀[l], z .- z₀[l]), Ref(R)) , Ref(SVector{3, Float64}(a[l], b[l], c[l])))) .<= 1.0)
 
+    end
     return P
 end
 shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Int(N), Int(M), 1; kwargs...)