From ebce189ece508321e02a4db3511adf0d33e5efea Mon Sep 17 00:00:00 2001 From: hzarei4 Date: Thu, 30 May 2024 19:05:11 +0200 Subject: [PATCH 1/4] added 3D Shepp-Logan support --- Project.toml | 3 ++- src/TestImages.jl | 66 ++++++++++++++++++++++++++++++++++++----------- 2 files changed, 53 insertions(+), 16 deletions(-) diff --git a/Project.toml b/Project.toml index 6d75b0c..44e6b30 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" +NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" StringDistances = "88034a9c-02f8-509d-84a9-84ec65e18404" @@ -15,9 +16,9 @@ StringDistances = "88034a9c-02f8-509d-84a9-84ec65e18404" [compat] AxisArrays = "0.3, 0.4" ColorTypes = "0.8, 0.9, 0.10, 0.11" +FileIO = "1" ImageIO = "0.3, 0.4, 0.5, 0.6" ImageMagick = "1" -FileIO = "1" OffsetArrays = "1" StringDistances = "0.10, 0.11" julia = "1.3" diff --git a/src/TestImages.jl b/src/TestImages.jl index 00e5075..e7b5770 100644 --- a/src/TestImages.jl +++ b/src/TestImages.jl @@ -4,6 +4,7 @@ using Pkg.Artifacts using StringDistances using ColorTypes using ColorTypes.FixedPointNumbers +using NDTools: reorient const artifacts_toml = abspath(joinpath(@__DIR__, "..", "Artifacts.toml")) export testimage, testimage_dip3e @@ -197,14 +198,22 @@ MRI version [2] is calculated. [3] Jain, Anil K. Fundamentals of digital image processing. _Prentice-Hall, Inc._, (1989): 439. """ -function shepp_logan(N::Int, M::Int; high_contrast::Bool=true, highContrast=nothing) +function shepp_logan(N::Int, M::Int, O::Int; high_contrast::Bool=true, highContrast=nothing) + #println("shepp_logan function called") if !isnothing(highContrast) # compatbitity to Images.shepp_logan # remove this when we remove Images.shepp_logan Base.depwarn("keyword `highContrast` is deprecated, use `high_contrast` instead.", :shepp_logan) end - x = Array(range(-1, stop=1, length=M)') - y = Array(range(1, stop=-1, length=N)) + x = reorient(Array(range(-1, stop=1, length=M)), Val(1)) + y = reorient(Array(range(-1, stop=1, length=N)), Val(2)) + if O == 1 + z = [0.0] + else + z = reorient(Array(range(-1, stop=1, length=O)), Val(3)) + #println("z: ", z) + end + #z = Array(range(-1, stop=1, length=O)) # follow the notation in [2] A = high_contrast ? @@ -215,11 +224,15 @@ function shepp_logan(N::Int, M::Int; high_contrast::Bool=true, highContrast=noth # [3] p.439 uses the following setting for the CT version # and is used by MATLAB's built-in `phantom` with method `Shepp-Logan` # (1.0 , -0.98 , -0.02 , -0.02 , 0.01, 0.01 , 0.01 , 0.01 , 0.01 , 0.01 ) - x₀ = (0.0 , 0.0 , 0.22 , -0.22 , 0.0 , 0.0 , 0.0 , -0.08 , 0.0 , 0.06 ) - y₀ = (0.0 , -0.0184, 0.0 , 0.0 , 0.35, 0.1 , -0.1 , -0.605, -0.605, -0.605) - a = (0.69, 0.6624, 0.11 , 0.16 , 0.21, 0.046, 0.046, 0.046, 0.023, 0.023) - b = (0.92, 0.874 , 0.31 , 0.41 , 0.25, 0.046, 0.046, 0.023, 0.023, 0.046) - ϕ = (0.0 , 0.0 , -18.0 , 18.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ) + x₀ = (0.0 , 0.0 , 0.22 , -0.22 , 0.0 , 0.0 , 0.0 , -0.08 , 0.0 , 0.06 ) + y₀ = (0.0 , -0.0184, 0.0 , 0.0 , 0.35 , 0.1 , -0.1 , -0.605, -0.605, -0.605) + z₀ = (0.0 , 0.0 , 0.0 , 0.0 , -0.15 , 0.25 , 0.25 , 0.0 , 0.0 , 0.0) + a = (0.69, 0.6624, 0.11 , 0.16 , 0.21 , 0.046, 0.046, 0.046, 0.023, 0.023) + b = (0.92, 0.874 , 0.31 , 0.41 , 0.25 , 0.046, 0.046, 0.023, 0.023, 0.046) + c = (0.81, 0.780 , 0.22 , 0.28 , 0.41 , 0.050, 0.050, 0.050, 0.020, 0.020) + ϕ = (0.0 , 0.0 , -18.0 , 18.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ) + θ = (0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ) + ψ = (0.0 , 0.0 , 10.0 , 10.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ) function _ellipse(dx, dy, a, b, sin_ϕ, cos_ϕ) tx = cos_ϕ * dx + sin_ϕ * dy @@ -230,21 +243,44 @@ function shepp_logan(N::Int, M::Int; high_contrast::Bool=true, highContrast=noth # 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] + + # 3D ellipse equation + (abs2(tx / a) + abs2(ty / b) + abs2(tz / c)) <= 1.0 + end - P = zeros(Gray{Float64}, N, M) + P = zeros(Float64, N, M, O) for l = 1:length(ϕ) - if ϕ[l] == 0.0 - @. P = gray(P) + A[l] * _ellipse(x - x₀[l], y - y₀[l], a[l], b[l]) + if ϕ[l] == 0.0 && θ[l] == 0.0 && ψ[l] == 0.0 + # Rotation matrix using Z-Y-Z Euler angles + R = [ + 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]) - @. P = gray(P) + A[l] * _ellipse(x - x₀[l], y - y₀[l], a[l], b[l], sin_ϕ, cos_ϕ) + sinϕ, cosϕ = sincosd(ϕ[l]) + sinθ, cosθ = sincosd(θ[l]) + sinψ, cosψ = sincosd(ψ[l]) + # Rotation matrix using Z-Y-Z Euler angles + R = [ + (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 return P end -shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Int(N), Int(M); kwargs...) -shepp_logan(N::Integer; kwargs...) = shepp_logan(Int(N), Int(N); kwargs...) +shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Int(N), Int(M), 1; kwargs...) +shepp_logan(N::Integer; kwargs...) = shepp_logan(Int(N), Int(N), 1; kwargs...) function _precompile_() ccall(:jl_generating_output, Cint, ()) == 1 || return nothing From ae70d2693e0247afeb23133e321bc4a85a618012 Mon Sep 17 00:00:00 2001 From: hzarei4 Date: Mon, 10 Jun 2024 16:35:17 +0200 Subject: [PATCH 2/4] optimized 3D shepp-logan --- Project.toml | 1 + src/TestImages.jl | 29 ++++++++++++++++------------- 2 files changed, 17 insertions(+), 13 deletions(-) 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...) From c965784fbe962720cca44da0ee6c5b26905de333 Mon Sep 17 00:00:00 2001 From: hzarei4 Date: Tue, 11 Jun 2024 14:18:30 +0200 Subject: [PATCH 3/4] added type modification --- src/TestImages.jl | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/TestImages.jl b/src/TestImages.jl index d6882b4..0dfd374 100644 --- a/src/TestImages.jl +++ b/src/TestImages.jl @@ -199,7 +199,7 @@ MRI version [2] is calculated. [3] Jain, Anil K. Fundamentals of digital image processing. _Prentice-Hall, Inc._, (1989): 439. """ -function shepp_logan(N::Int, M::Int, O::Int; high_contrast::Bool=true, highContrast=nothing) +function shepp_logan(::Type{T}, N::Int, M::Int, O::Int; high_contrast::Bool=true, highContrast=nothing) where {T} #println("shepp_logan function called") if !isnothing(highContrast) # compatbitity to Images.shepp_logan @@ -235,17 +235,6 @@ function shepp_logan(N::Int, M::Int, O::Int; high_contrast::Bool=true, highContr θ = (0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ) ψ = (0.0 , 0.0 , 10.0 , 10.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ) - function _ellipse(dx, dy, a, b, sin_ϕ, cos_ϕ) - tx = cos_ϕ * dx + sin_ϕ * dy - ty = sin_ϕ * dx - cos_ϕ * dy - abs2(tx * b) + abs2(ty * a) < (a * b)^2 - end - function _ellipse(dx, dy, a, b) - # a faster case when ϕ == 0.0 - abs2(dx * b) + abs2(dy * a) < (a * b)^2 - end - - # 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 @@ -257,20 +246,20 @@ function shepp_logan(N::Int, M::Int, O::Int; high_contrast::Bool=true, highContr return sum(abs2.(t)) end - P = zeros(Float64, N, M, O) + P = zeros(T, 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 = SMatrix{3, 3, Float64}([ - 1.0 0.0 0.0; - 0.0 1.0 0.0; - 0.0 0.0 1.0 + 1.0 0.0 0.0; + 0.0 1.0 0.0; + 0.0 0.0 1.0 ]) else sinϕ, cosϕ = sincosd(ϕ[l]) sinθ, cosθ = sincosd(θ[l]) sinψ, cosψ = sincosd(ψ[l]) - # Rotation matrix using Z-Y-Z Euler angles + 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θ); @@ -282,8 +271,9 @@ function shepp_logan(N::Int, M::Int, O::Int; high_contrast::Bool=true, highContr end return P end -shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Int(N), Int(M), 1; kwargs...) -shepp_logan(N::Integer; kwargs...) = shepp_logan(Int(N), Int(N), 1; kwargs...) +shepp_logan(N::Integer, M::Integer, O::Integer; kwargs...) = shepp_logan(Float64, Int(N), Int(M), Int(O); kwargs...) +shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Float64, Int(N), Int(M), 1; kwargs...)[:, :, 1] +shepp_logan(N::Integer; kwargs...) = shepp_logan(Float64, Int(N), Int(N), 1; kwargs...)[:, :, 1] function _precompile_() ccall(:jl_generating_output, Cint, ()) == 1 || return nothing From 2fac62180fbeb9da47d02fe76a3192783d5fa72c Mon Sep 17 00:00:00 2001 From: hzarei4 Date: Tue, 11 Jun 2024 15:21:31 +0200 Subject: [PATCH 4/4] N, M size correction --- src/TestImages.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/TestImages.jl b/src/TestImages.jl index 0dfd374..5afedab 100644 --- a/src/TestImages.jl +++ b/src/TestImages.jl @@ -206,8 +206,8 @@ function shepp_logan(::Type{T}, N::Int, M::Int, O::Int; high_contrast::Bool=true # remove this when we remove Images.shepp_logan Base.depwarn("keyword `highContrast` is deprecated, use `high_contrast` instead.", :shepp_logan) end - x = reorient(Array(range(-1, stop=1, length=M)), Val(1)) - y = reorient(Array(range(-1, stop=1, length=N)), Val(2)) + x = reorient(Array(range(-1, stop=1, length=N)), Val(1)) + y = reorient(Array(range(-1, stop=1, length=M)), Val(2)) if O == 1 z = [0.0] else @@ -265,15 +265,16 @@ function shepp_logan(::Type{T}, N::Int, M::Int, O::Int; high_contrast::Bool=true (-sinψ*cosϕ-cosθ*cosψ*sinϕ) (-sinψ*sinϕ+cosθ*cosψ*cosϕ) (cosψ*sinθ); (sinθ*sinϕ) (-sinθ*cosϕ) (cosθ) ]) + 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) + P .+= T(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, O::Integer; kwargs...) = shepp_logan(Float64, Int(N), Int(M), Int(O); kwargs...) -shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Float64, Int(N), Int(M), 1; kwargs...)[:, :, 1] -shepp_logan(N::Integer; kwargs...) = shepp_logan(Float64, Int(N), Int(N), 1; kwargs...)[:, :, 1] +shepp_logan(N::Integer, M::Integer, O::Integer; kwargs...) = shepp_logan(Gray{Float64}, Int(N), Int(M), Int(O); kwargs...) +shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Gray{Float64}, Int(N), Int(M), 1; kwargs...)[:, :, 1] +shepp_logan(N::Integer; kwargs...) = shepp_logan(Gray{Float64}, Int(N), Int(N), 1; kwargs...)[:, :, 1] function _precompile_() ccall(:jl_generating_output, Cint, ()) == 1 || return nothing