diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl index 192b8ddd1..3fd03d493 100644 --- a/KomaMRIBase/src/KomaMRIBase.jl +++ b/KomaMRIBase/src/KomaMRIBase.jl @@ -10,10 +10,11 @@ using Parameters using Interpolations #Reconstruction using MRIBase -@reexport using MRIBase: Profile, RawAcquisitionData, AcquisitionData, AcquisitionHeader, EncodingCounters, Limit +@reexport using MRIBase: + Profile, RawAcquisitionData, AcquisitionData, AcquisitionHeader, EncodingCounters, Limit using MAT # For loading example phantoms -global γ = 42.5774688e6; # Hz/T gyromagnetic constant for H1, JEMRIS uses 42.5756 MHz/T +global γ = 42.5774688e6 # Hz/T gyromagnetic constant for H1, JEMRIS uses 42.5756 MHz/T # Hardware include("datatypes/Scanner.jl") @@ -43,7 +44,15 @@ export times, ampls, freqs # This are also used for simulation export kfoldperm, trapz, cumtrapz # Phantom -export brain_phantom2D, brain_phantom3D, pelvis_phantom2D +export brain_phantom2D, brain_phantom3D, pelvis_phantom2D, heart_phantom +# Motion +export MotionModel +export NoMotion, SimpleMotion, ArbitraryMotion +export SimpleMotionType +export Translation, Rotation, HeartBeat +export PeriodicTranslation, PeriodicRotation, PeriodicHeartBeat +export get_spin_coords, sort_motions! +export LinearInterpolator # Secondary export get_kspace, rotx, roty, rotz # Additionals @@ -57,6 +66,8 @@ export PulseDesigner #Package version, KomaMRIBase.__VERSION__ using Pkg -__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]) +__VERSION__ = VersionNumber( + Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"] +) end # module KomaMRIBase diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl index a125cb438..62c1ac81f 100644 --- a/KomaMRIBase/src/datatypes/Phantom.jl +++ b/KomaMRIBase/src/datatypes/Phantom.jl @@ -1,5 +1,12 @@ +abstract type MotionModel{T<:Real} end + +#Motion models: +include("phantom/motion/SimpleMotion.jl") +include("phantom/motion/ArbitraryMotion.jl") +include("phantom/motion/NoMotion.jl") + """ - obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, ux, uy, uz) + obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, motion) The Phantom struct. Most of its field names are vectors, with each element associated with a property value representing a spin. This struct serves as an input for the simulation. @@ -17,9 +24,7 @@ a property value representing a spin. This struct serves as an input for the sim - `Dλ1`: (`::AbstractVector{T<:Real}`) spin Dλ1 (diffusion) parameter vector - `Dλ2`: (`::AbstractVector{T<:Real}`) spin Dλ2 (diffusion) parameter vector - `Dθ`: (`::AbstractVector{T<:Real}`) spin Dθ (diffusion) parameter vector -- `ux`: (`::Function`) displacement field in the x-axis -- `uy`: (`::Function`) displacement field in the y-axis -- `uz`: (`::Function`) displacement field in the z-axis +- `motion`: (`::MotionModel{T<:Real}`) motion model # Returns - `obj`: (`::Phantom`) Phantom struct @@ -31,27 +36,25 @@ julia> obj = Phantom(x=[0.0]) julia> obj.ρ ``` """ - @with_kw mutable struct Phantom{T<:Real} - name::String = "spins" - x::AbstractVector{T} - y::AbstractVector{T} = zeros(size(x)) - z::AbstractVector{T} = zeros(size(x)) - ρ::AbstractVector{T} = ones(size(x)) - T1::AbstractVector{T} = ones(size(x)) * 1_000_000 - T2::AbstractVector{T} = ones(size(x)) * 1_000_000 - T2s::AbstractVector{T} = ones(size(x)) * 1_000_000 +@with_kw mutable struct Phantom{T<:Real} + name::String = "spins" + x :: AbstractVector{T} + y::AbstractVector{T} = zeros(eltype(x), size(x)) + z::AbstractVector{T} = zeros(eltype(x), size(x)) + ρ::AbstractVector{T} = ones(eltype(x), size(x)) + T1::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000 + T2::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000 + T2s::AbstractVector{T} = ones(eltype(x), size(x)) * 1_000_000 #Off-resonance related - Δw::AbstractVector{T} = zeros(size(x)) + Δw::AbstractVector{T} = zeros(eltype(x), size(x)) #χ::Vector{SusceptibilityModel} #Diffusion - Dλ1::AbstractVector{T} = zeros(size(x)) - Dλ2::AbstractVector{T} = zeros(size(x)) - Dθ::AbstractVector{T} = zeros(size(x)) + Dλ1::AbstractVector{T} = zeros(eltype(x), size(x)) + Dλ2::AbstractVector{T} = zeros(eltype(x), size(x)) + Dθ::AbstractVector{T} = zeros(eltype(x), size(x)) #Diff::Vector{DiffusionModel} #Diffusion map #Motion - ux::Function = (x,y,z,t)->0 - uy::Function = (x,y,z,t)->0 - uz::Function = (x,y,z,t)->0 + motion::MotionModel{T} = NoMotion{eltype(x)}() end """Size and length of a phantom""" @@ -59,117 +62,144 @@ size(x::Phantom) = size(x.ρ) Base.length(x::Phantom) = length(x.ρ) # To enable to iterate and broadcast over the Phantom Base.iterate(x::Phantom) = (x[1], 2) -Base.iterate(x::Phantom, i::Integer) = (i <= length(x)) ? (x[i], i+1) : nothing +Base.iterate(x::Phantom, i::Integer) = (i <= length(x)) ? (x[i], i + 1) : nothing Base.lastindex(x::Phantom) = length(x) Base.getindex(x::Phantom, i::Integer) = x[i:i] """Compare two phantoms""" -Base.isapprox(obj1::Phantom, obj2::Phantom) = begin - obj1.x ≈ obj2.x && - obj1.y ≈ obj2.y && - obj1.z ≈ obj2.z && - obj1.ρ ≈ obj2.ρ && - obj1.T1 ≈ obj2.T1 && - obj1.T2 ≈ obj2.T2 && - obj1.T2s ≈ obj2.T2s && - obj1.Δw ≈ obj2.Δw && - obj1.Dλ1 ≈ obj2.Dλ1 && - obj1.Dλ2 ≈ obj2.Dλ2 && - obj1.Dθ ≈ obj2.Dθ -end +Base.:(==)(obj1::Phantom, obj2::Phantom) = reduce( + &, + [ + getfield(obj1, field) == getfield(obj2, field) for + field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom)) + ], +) +Base.:(≈)(obj1::Phantom, obj2::Phantom) = reduce(&, [getfield(obj1, field) ≈ getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))]) +Base.:(==)(m1::MotionModel, m2::MotionModel) = false +Base.:(≈)(m1::MotionModel, m2::MotionModel) = false """ Separate object spins in a sub-group """ -Base.getindex(obj::Phantom, p::AbstractRange) = begin - Phantom(name=obj.name, - x=obj.x[p], - y=obj.y[p], - z=obj.z[p], - ρ=obj.ρ[p], - T1=obj.T1[p], - T2=obj.T2[p], - T2s=obj.T2s[p], - Δw=obj.Δw[p], - #Diff=obj.Diff[p], #TODO! - Dλ1=obj.Dλ1[p], - Dλ2=obj.Dλ2[p], - Dθ=obj.Dθ[p], - #Χ=obj.Χ[p], #TODO! - ux=obj.ux, - uy=obj.uy, - uz=obj.uz - ) +Base.getindex(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begin + fields = [] + for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom)) + push!(fields, (field, getfield(obj, field)[p])) + end + return Phantom(; name=obj.name, fields...) end """Separate object spins in a sub-group (lightweigth).""" -Base.view(obj::Phantom, p::AbstractRange) = begin - @views Phantom(name=obj.name, - x=obj.x[p], - y=obj.y[p], - z=obj.z[p], - ρ=obj.ρ[p], - T1=obj.T1[p], - T2=obj.T2[p], - T2s=obj.T2s[p], - Δw=obj.Δw[p], - #Diff=obj.Diff[p], #TODO! - Dλ1=obj.Dλ1[p], - Dλ2=obj.Dλ2[p], - Dθ=obj.Dθ[p], - #Χ=obj.Χ[p], #TODO! - ux=obj.ux, - uy=obj.uy, - uz=obj.uz - ) -end +Base.view(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = @views obj[p] """Addition of phantoms""" -+(s1::Phantom,s2::Phantom) = begin - Phantom(name=s1.name*"+"*s2.name, - x=[s1.x;s2.x], - y=[s1.y;s2.y], - z=[s1.z;s2.z], - ρ=[s1.ρ;s2.ρ], - T1=[s1.T1;s2.T1], - T2=[s1.T2;s2.T2], - T2s=[s1.T2s;s2.T2s], - Δw=[s1.Δw;s2.Δw], - #Diff=obj.Diff[p], #TODO! - Dλ1=[s1.Dλ1;s2.Dλ1], - Dλ2=[s1.Dλ2;s2.Dλ2], - Dθ=[s1.Dθ;s2.Dθ], - #Χ=obj.Χ[p], #TODO! - ux=s1.ux, - uy=s1.uy, - uz=s1.uz - ) ++(obj1::Phantom, obj2::Phantom) = begin + fields = [] + for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom)) + push!(fields, (field, [getfield(obj1, field); getfield(obj2, field)])) + end + Nmaxchars = 50 + name = first(obj1.name * "+" * obj2.name, Nmaxchars) + return Phantom(; name=name, fields...) end """Scalar multiplication of a phantom""" -*(α::Real,obj::Phantom) = begin - Phantom(name=obj.name, - x=obj.x, - y=obj.y, - z=obj.z, - ρ=α*obj.ρ, #Only affects the proton density - T1=obj.T1, - T2=obj.T2, - T2s=obj.T2s, - Δw=obj.Δw, - #Diff=obj.Diff[p], #TODO! - Dλ1=obj.Dλ1, - Dλ2=obj.Dλ2, - Dθ=obj.Dθ, - #Χ=obj.Χ[p], #TODO! - ux=obj.ux, - uy=obj.uy, - uz=obj.uz +*(α::Real, obj::Phantom) = begin + obj1 = copy(obj) + obj1.ρ .*= α + return obj1 +end + +""" +dims = get_dims(obj) +""" +function get_dims(obj::Phantom) + dims = Bool[] + push!(dims, any(x -> x != 0, obj.x)) + push!(dims, any(x -> x != 0, obj.y)) + push!(dims, any(x -> x != 0, obj.z)) + return dims +end + +function sort_motions!(motion::MotionModel) + return nothing +end + +""" + obj = heart_phantom(...) + +Heart-like LV phantom. The variable `circumferential_strain` and `radial_strain` are for streching (if positive) +or contraction (if negative). `rotation_angle` is for rotation. + +# Arguments +- `circumferential_strain`: (`::Real`, `=-0.3`) contraction parameter +- `radial_strain`: (`::Real`, `=-0.3`) contraction parameter +- `rotation_angle`: (`::Real`, `=1`) rotation parameter + +# Returns +- `phantom`: (`::Phantom`) Heart-like LV phantom struct +""" +function heart_phantom( + circumferential_strain=-0.3, + radial_strain=-0.3, + rotation_angle=15.0; + heart_rate=60, + asymmetry=0.2, +) + #PARAMETERS + FOV = 10e-2 # [m] Diameter ventricule + N = 21 + Δxr = FOV / (N - 1) #Aprox rec resolution, use Δx_pix and Δy_pix + Ns = 50 #number of spins per voxel + Δx = Δxr / sqrt(Ns) #spin separation + #POSITIONS + x = y = (-FOV / 2):Δx:(FOV / 2 - Δx) #spin coordinates + x, y = x .+ y' * 0, x * 0 .+ y' #grid points + #PHANTOM + ⚪(R) = (x .^ 2 .+ y .^ 2 .<= R^2) * 1.0 #Circle of radius R + period = 60 / heart_rate # [s] Period + # Water spins + R = 9 / 10 * FOV / 2 + r = 6 / 11 * FOV / 2 + ring = ⚪(R) .- ⚪(r) + ρ = ⚪(r) .+ 0.9 * ring #proton density + # Diffusion tensor model + D = 2e-9 #Diffusion of free water m2/s + D1, D2 = D, D / 20 + Dλ1 = D1 * ⚪(R) #Diffusion map + Dλ2 = D1 * ⚪(r) .+ D2 * ring #Diffusion map + Dθ = atan.(x, -y) .* ring #Diffusion map + T1 = (1400 * ⚪(r) .+ 1026 * ring) * 1e-3 #Myocardial T1 + T2 = (308 * ⚪(r) .+ 42 * ring) * 1e-3 #T2 map [s] + # Generating Phantoms + phantom = Phantom(; + name="LeftVentricle", + x=x[ρ .!= 0], + y=y[ρ .!= 0], + ρ=ρ[ρ .!= 0], + T1=T1[ρ .!= 0], + T2=T2[ρ .!= 0], + Dλ1=Dλ1[ρ .!= 0], + Dλ2=Dλ2[ρ .!= 0], + Dθ=Dθ[ρ .!= 0], + motion=SimpleMotion([ + PeriodicHeartBeat(; + period=period, + asymmetry=asymmetry, + circumferential_strain=circumferential_strain, + radial_strain=radial_strain, + longitudinal_strain=0.0, + ), + PeriodicRotation(; + period=period, asymmetry=asymmetry, yaw=rotation_angle, pitch=0.0, roll=0.0 + ), + ]), ) + return phantom end """ - obj = brain_phantom2D(; axis="axial", ss=4, us=1) + phantom = brain_phantom2D(;axis="axial", ss=4) Creates a two-dimensional brain Phantom struct. Default ss=4 sample spacing is 2 mm. Original file (ss=1) sample spacing is .5 mm. @@ -200,91 +230,95 @@ julia> plot_phantom_map(obj, :ρ) ``` """ function brain_phantom2D(; axis="axial", ss=4, us=1) - # check and filter input ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(2, ss, us) # Get data from .mat file path = @__DIR__ - data = MAT.matread(path*"/phantom/brain2D.mat") + data = MAT.matread(path * "/phantom/brain2D.mat") # subsample or upsample the phantom data - class = repeat(data[axis][1:ssx:end,1:ssy:end], inner=[usx, usy]) - + class = repeat(data[axis][1:ssx:end, 1:ssy:end]; inner=[usx, usy]) + # Define spin position vectors - Δx = .5e-3*ssx/usx - Δy = .5e-3*ssy/usy + Δx = .5e-3 * ssx / usx + Δy = .5e-3 * ssy / usy M, N = size(class) - FOVx = (M-1)*Δx #[m] - FOVy = (N-1)*Δy #[m] - x = -FOVx/2:Δx:FOVx/2 #spin coordinates - y = -FOVy/2:Δy:FOVy/2 #spin coordinates - x, y = x .+ y'*0, x*0 .+ y' #grid points + FOVx = (M - 1) * Δx #[m] + FOVy = (N - 1) * Δy #[m] + x = (-FOVx / 2):Δx:(FOVx / 2) #spin coordinates + y = (-FOVy / 2):Δy:(FOVy / 2) #spin coordinates + x, y = x .+ y' * 0, x * 0 .+ y' #grid points # Define spin property vectors - T2 = (class.==23)*329 .+ #CSF - (class.==46)*83 .+ #GM - (class.==70)*70 .+ #WM - (class.==93)*70 .+ #FAT1 - (class.==116)*47 .+ #MUSCLE - (class.==139)*329 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*70 .+ #FAT2 - (class.==232)*329 .+ #DURA - (class.==255)*70 #MARROW - T2s = (class.==23)*58 .+ #CSF - (class.==46)*69 .+ #GM - (class.==70)*61 .+ #WM - (class.==93)*58 .+ #FAT1 - (class.==116)*30 .+ #MUSCLE - (class.==139)*58 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*61 .+ #FAT2 - (class.==232)*58 .+ #DURA - (class.==255)*61 .+#MARROW - (class.==255)*70 #MARROW - T1 = (class.==23)*2569 .+ #CSF - (class.==46)*833 .+ #GM - (class.==70)*500 .+ #WM - (class.==93)*350 .+ #FAT1 - (class.==116)*900 .+ #MUSCLE - (class.==139)*569 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*500 .+ #FAT2 - (class.==232)*2569 .+ #DURA - (class.==255)*500 #MARROW - ρ = (class.==23)*1 .+ #CSF - (class.==46)*.86 .+ #GM - (class.==70)*.77 .+ #WM - (class.==93)*1 .+ #FAT1 - (class.==116)*1 .+ #MUSCLE - (class.==139)*.7 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*.77 .+ #FAT2 - (class.==232)*1 .+ #DURA - (class.==255)*.77 #MARROW - Δw_fat = -220*2π - Δw = (class.==93)*Δw_fat .+ #FAT1 - (class.==209)*Δw_fat #FAT2 - T1 = T1*1e-3 - T2 = T2*1e-3 - T2s = T2s*1e-3 + T2 = + (class .== 23) * 329 .+ #CSF + (class .== 46) * 83 .+ #GM + (class .== 70) * 70 .+ #WM + (class .== 93) * 70 .+ #FAT1 + (class .== 116) * 47 .+ #MUSCLE + (class .== 139) * 329 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 70 .+ #FAT2 + (class .== 232) * 329 .+ #DURA + (class .== 255) * 70 #MARROW + T2s = + (class .== 23) * 58 .+ #CSF + (class .== 46) * 69 .+ #GM + (class .== 70) * 61 .+ #WM + (class .== 93) * 58 .+ #FAT1 + (class .== 116) * 30 .+ #MUSCLE + (class .== 139) * 58 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 61 .+ #FAT2 + (class .== 232) * 58 .+ #DURA + (class .== 255) * 61 .+#MARROW + (class .== 255) * 70 #MARROW + T1 = + (class .== 23) * 2569 .+ #CSF + (class .== 46) * 833 .+ #GM + (class .== 70) * 500 .+ #WM + (class .== 93) * 350 .+ #FAT1 + (class .== 116) * 900 .+ #MUSCLE + (class .== 139) * 569 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 500 .+ #FAT2 + (class .== 232) * 2569 .+ #DURA + (class .== 255) * 500 #MARROW + ρ = + (class .== 23) * 1 .+ #CSF + (class .== 46) * 0.86 .+ #GM + (class .== 70) * 0.77 .+ #WM + (class .== 93) * 1 .+ #FAT1 + (class .== 116) * 1 .+ #MUSCLE + (class .== 139) * 0.7 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 0.77 .+ #FAT2 + (class .== 232) * 1 .+ #DURA + (class .== 255) * 0.77 #MARROW + Δw_fat = -220 * 2π + Δw = + (class .== 93) * Δw_fat .+ #FAT1 + (class .== 209) * Δw_fat #FAT2 + T1 = T1 * 1e-3 + T2 = T2 * 1e-3 + T2s = T2s * 1e-3 # Define and return the Phantom struct - obj = Phantom{Float64}( - name = "brain2D_"*axis, - x = y[ρ.!=0], - y = x[ρ.!=0], - z = 0*x[ρ.!=0], - ρ = ρ[ρ.!=0], - T1 = T1[ρ.!=0], - T2 = T2[ρ.!=0], - T2s = T2s[ρ.!=0], - Δw = Δw[ρ.!=0], + obj = Phantom{Float64}(; + name="brain2D_" * axis, + x=y[ρ .!= 0], + y=x[ρ .!= 0], + z=0 * x[ρ .!= 0], + ρ=ρ[ρ .!= 0], + T1=T1[ρ .!= 0], + T2=T2[ρ .!= 0], + T2s=T2s[ρ .!= 0], + Δw=Δw[ρ .!= 0], ) return obj end @@ -319,97 +353,104 @@ julia> obj = brain_phantom3D(; us=[2, 2, 1]) julia> plot_phantom_map(obj, :ρ) ``` """ -function brain_phantom3D(;ss=4, us=1, start_end=[160, 200]) - +function brain_phantom3D(; ss=4, us=1, start_end=[160, 200]) # check and filter input ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(3, ss, us) # Get data from .mat file path = @__DIR__ - data = MAT.matread(path*"/phantom/brain3D.mat") + data = MAT.matread(path * "/phantom/brain3D.mat") # subsample or upsample the phantom data - class = repeat(data["data"][1:ssx:end,1:ssy:end, start_end[1]:ssz:start_end[2]], inner=[usx, usy, usz]) + class = repeat( + data["data"][1:ssx:end, 1:ssy:end, start_end[1]:ssz:start_end[2]]; + inner=[usx, usy, usz], + ) # Define spin position vectors - Δx = .5e-3*ssx/usx - Δy = .5e-3*ssy/usy - Δz = .5e-3*ssz/usz + Δx = .5e-3 * ssx / usx + Δy = .5e-3 * ssy / usy + Δz = .5e-3 * ssz / usz M, N, Z = size(class) - FOVx = (M-1)*Δx #[m] - FOVy = (N-1)*Δy #[m] - FOVz = (Z-1)*Δz #[m] - xx = reshape(-FOVx/2:Δx:FOVx/2,M,1,1) #spin coordinates - yy = reshape(-FOVy/2:Δy:FOVy/2,1,N,1) #spin coordinates - zz = reshape(-FOVz/2:Δz:FOVz/2,1,1,Z) #spin coordinates - x = 1*xx .+ 0*yy .+ 0*zz - y = 0*xx .+ 1*yy .+ 0*zz - z = 0*xx .+ 0*yy .+ 1*zz + FOVx = (M - 1) * Δx #[m] + FOVy = (N - 1) * Δy #[m] + FOVz = (Z - 1) * Δz #[m] + xx = reshape((-FOVx / 2):Δx:(FOVx / 2), M, 1, 1) #spin coordinates + yy = reshape((-FOVy / 2):Δy:(FOVy / 2), 1, N, 1) #spin coordinates + zz = reshape((-FOVz / 2):Δz:(FOVz / 2), 1, 1, Z) #spin coordinates + x = 1 * xx .+ 0 * yy .+ 0 * zz + y = 0 * xx .+ 1 * yy .+ 0 * zz + z = 0 * xx .+ 0 * yy .+ 1 * zz # Define spin property vectors - T2 = (class.==23)*329 .+ #CSF - (class.==46)*83 .+ #GM - (class.==70)*70 .+ #WM - (class.==93)*70 .+ #FAT1 - (class.==116)*47 .+ #MUSCLE - (class.==139)*329 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*70 .+ #FAT2 - (class.==232)*329 .+ #DURA - (class.==255)*70 #MARROW - T2s = (class.==23)*58 .+ #CSF - (class.==46)*69 .+ #GM - (class.==70)*61 .+ #WM - (class.==93)*58 .+ #FAT1 - (class.==116)*30 .+ #MUSCLE - (class.==139)*58 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*61 .+ #FAT2 - (class.==232)*58 .+ #DURA - (class.==255)*61 .+#MARROW - (class.==255)*70 #MARROW - T1 = (class.==23)*2569 .+ #CSF - (class.==46)*833 .+ #GM - (class.==70)*500 .+ #WM - (class.==93)*350 .+ #FAT1 - (class.==116)*900 .+ #MUSCLE - (class.==139)*569 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*500 .+ #FAT2 - (class.==232)*2569 .+ #DURA - (class.==255)*500 #MARROW - ρ = (class.==23)*1 .+ #CSF - (class.==46)*.86 .+ #GM - (class.==70)*.77 .+ #WM - (class.==93)*1 .+ #FAT1 - (class.==116)*1 .+ #MUSCLE - (class.==139)*.7 .+ #SKIN/MUSCLE - (class.==162)*0 .+ #SKULL - (class.==185)*0 .+ #VESSELS - (class.==209)*.77 .+ #FAT2 - (class.==232)*1 .+ #DURA - (class.==255)*.77 #MARROW - Δw_fat = -220*2π - Δw = (class.==93)*Δw_fat .+ #FAT1 - (class.==209)*Δw_fat #FAT2 - T1 = T1*1e-3 - T2 = T2*1e-3 - T2s = T2s*1e-3 + T2 = + (class .== 23) * 329 .+ #CSF + (class .== 46) * 83 .+ #GM + (class .== 70) * 70 .+ #WM + (class .== 93) * 70 .+ #FAT1 + (class .== 116) * 47 .+ #MUSCLE + (class .== 139) * 329 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 70 .+ #FAT2 + (class .== 232) * 329 .+ #DURA + (class .== 255) * 70 #MARROW + T2s = + (class .== 23) * 58 .+ #CSF + (class .== 46) * 69 .+ #GM + (class .== 70) * 61 .+ #WM + (class .== 93) * 58 .+ #FAT1 + (class .== 116) * 30 .+ #MUSCLE + (class .== 139) * 58 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 61 .+ #FAT2 + (class .== 232) * 58 .+ #DURA + (class .== 255) * 61 .+#MARROW + (class .== 255) * 70 #MARROW + T1 = + (class .== 23) * 2569 .+ #CSF + (class .== 46) * 833 .+ #GM + (class .== 70) * 500 .+ #WM + (class .== 93) * 350 .+ #FAT1 + (class .== 116) * 900 .+ #MUSCLE + (class .== 139) * 569 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 500 .+ #FAT2 + (class .== 232) * 2569 .+ #DURA + (class .== 255) * 500 #MARROW + ρ = + (class .== 23) * 1 .+ #CSF + (class .== 46) * 0.86 .+ #GM + (class .== 70) * 0.77 .+ #WM + (class .== 93) * 1 .+ #FAT1 + (class .== 116) * 1 .+ #MUSCLE + (class .== 139) * 0.7 .+ #SKIN/MUSCLE + (class .== 162) * 0 .+ #SKULL + (class .== 185) * 0 .+ #VESSELS + (class .== 209) * 0.77 .+ #FAT2 + (class .== 232) * 1 .+ #DURA + (class .== 255) * 0.77 #MARROW + Δw_fat = -220 * 2π + Δw = + (class .== 93) * Δw_fat .+ #FAT1 + (class .== 209) * Δw_fat #FAT2 + T1 = T1 * 1e-3 + T2 = T2 * 1e-3 + T2s = T2s * 1e-3 # Define and return the Phantom struct - obj = Phantom{Float64}( - name = "brain3D", - x = y[ρ.!=0], - y = x[ρ.!=0], - z = z[ρ.!=0], - ρ = ρ[ρ.!=0], - T1 = T1[ρ.!=0], - T2 = T2[ρ.!=0], - T2s = T2s[ρ.!=0], - Δw = Δw[ρ.!=0], + obj = Phantom{Float64}(; + name="brain3D", + x=y[ρ .!= 0], + y=x[ρ .!= 0], + z=z[ρ .!= 0], + ρ=ρ[ρ .!= 0], + T1=T1[ρ .!= 0], + T2=T2[ρ .!= 0], + T2s=T2s[ρ .!= 0], + Δw=Δw[ρ .!= 0], ) return obj end @@ -437,61 +478,64 @@ julia> pelvis_phantom2D(obj, :ρ) ``` """ function pelvis_phantom2D(; ss=4, us=1) - # check and filter input ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(2, ss, us) # Get data from .mat file path = @__DIR__ - data = MAT.matread(path*"/phantom/pelvis2D.mat") - + data = MAT.matread(path * "/phantom/pelvis2D.mat") + # subsample or upsample the phantom data - class = repeat(data["pelvis3D_slice"][1:ssx:end,1:ssy:end], inner=[usx, usy]) + class = repeat(data["pelvis3D_slice"][1:ssx:end, 1:ssy:end]; inner=[usx, usy]) # Define spin position vectors - Δx = .5e-3*ssx/usx - Δy = .5e-3*ssy/usy + Δx = .5e-3 * ssx / usx + Δy = .5e-3 * ssy / usy M, N = size(class) - FOVx = (M-1)*Δx # [m] - FOVy = (N-1)*Δy # [m] - x = -FOVx/2:Δx:FOVx/2 # spin coordinates - y = -FOVy/2:Δy:FOVy/2 # spin coordinates - x, y = x .+ y'*0, x*0 .+ y' # grid points + FOVx = (M - 1) * Δx # [m] + FOVy = (N - 1) * Δy # [m] + x = (-FOVx / 2):Δx:(FOVx / 2) # spin coordinates + y = (-FOVy / 2):Δy:(FOVy / 2) # spin coordinates + x, y = x .+ y' * 0, x * 0 .+ y' # grid points # Define spin property vectors - ρ = (class.==102)*.86 .+ # Fat - (class.==153)*.9 .+ # SoftTissue - (class.==204)*.4 .+ # SpongyBone - (class.==255)*.2 # CorticalBone - T1 = (class.==102)*366 .+ # Fat - (class.==153)*1200 .+ # SoftTissue - (class.==204)*381 .+ # SpongyBone - (class.==255)*100 # CorticalBone - T2 = (class.==102)*70 .+ # Fat - (class.==153)*80 .+ # SoftTissue - (class.==204)*52 .+ # SpongyBone - (class.==255)*.3 # CorticalBone - T2s = (class.==102)*70 .+ # Fat - (class.==153)*80 .+ # SoftTissue - (class.==204)*52 .+ # SpongyBone - (class.==255)*.3 # CorticalBone + ρ = + (class .== 102) * 0.86 .+ # Fat + (class .== 153) * 0.9 .+ # SoftTissue + (class .== 204) * 0.4 .+ # SpongyBone + (class .== 255) * 0.2 # CorticalBone + T1 = + (class .== 102) * 366 .+ # Fat + (class .== 153) * 1200 .+ # SoftTissue + (class .== 204) * 381 .+ # SpongyBone + (class .== 255) * 100 # CorticalBone + T2 = + (class .== 102) * 70 .+ # Fat + (class .== 153) * 80 .+ # SoftTissue + (class .== 204) * 52 .+ # SpongyBone + (class .== 255) * 0.3 # CorticalBone + T2s = + (class .== 102) * 70 .+ # Fat + (class .== 153) * 80 .+ # SoftTissue + (class .== 204) * 52 .+ # SpongyBone + (class .== 255) * 0.3 # CorticalBone Δw_fat = -220 * 2π - Δw = (class.==102) * Δw_fat # FAT1 - T1 = T1*1e-3 - T2 = T2*1e-3 - T2s = T2s*1e-3 + Δw = (class .== 102) * Δw_fat # FAT1 + T1 = T1 * 1e-3 + T2 = T2 * 1e-3 + T2s = T2s * 1e-3 # Define and return the Phantom struct - obj = Phantom{Float64}( - name = "pelvis2D", - x = y[ρ.!=0], - y = x[ρ.!=0], - z = 0*x[ρ.!=0], - ρ = ρ[ρ.!=0], - T1 = T1[ρ.!=0], - T2 = T2[ρ.!=0], - T2s = T2s[ρ.!=0], - Δw = Δw[ρ.!=0], + obj = Phantom{Float64}(; + name="pelvis2D", + x=y[ρ .!= 0], + y=x[ρ .!= 0], + z=0 * x[ρ .!= 0], + ρ=ρ[ρ .!= 0], + T1=T1[ρ .!= 0], + T2=T2[ρ .!= 0], + T2s=T2s[ρ .!= 0], + Δw=Δw[ρ .!= 0], ) return obj end @@ -516,8 +560,7 @@ julia> ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(2, 1, 1) julia> ssx, ssy, ssz, usx, usy, usz = check_phantom_arguments(3, 4, [2, 2, 2]) ``` """ -function check_phantom_arguments( nd, ss, us) - +function check_phantom_arguments(nd, ss, us) # check for valid input ssz = -9999 usz = -9999 @@ -529,33 +572,49 @@ function check_phantom_arguments( nd, ss, us) @assert length(ss) <= 3 "ss=$(ss) invalid, ss can have up to three components [ssx, ssy, ssz] for a 3D phantom" @assert length(us) <= 3 "us=$(us) invalid, us can have up to three components [usx, usy, usz] for a 3D phantom" if length(us) == 1 - usx = us[1]; usy = us[1]; usz = us[1] - elseif length( us) == 2 - usx = us[1]; usy = us[2]; usz = us[2] + usx = us[1] + usy = us[1] + usz = us[1] + elseif length(us) == 2 + usx = us[1] + usy = us[2] + usz = us[2] @warn "Using us=$([usx, usy, usz]) in place of us=$([usx, usy])." else - usx = us[1]; usy = us[2]; usz = us[3] + usx = us[1] + usy = us[2] + usz = us[3] end if length(ss) == 1 - ssx = ss[1]; ssy = ss[1]; ssz = ss[1] - elseif length( ss) == 2 - ssx = ss[1]; ssy = ss[2]; ssz = ss[2] + ssx = ss[1] + ssy = ss[1] + ssz = ss[1] + elseif length(ss) == 2 + ssx = ss[1] + ssy = ss[2] + ssz = ss[2] @warn "Using ss=$([ssx, ssy, ssz]) in place of ss=$([ssx, ssy])." else - ssx = ss[1]; ssy = ss[2]; ssz = ss[3] - end - elseif nd == 2 + ssx = ss[1] + ssy = ss[2] + ssz = ss[3] + end + elseif nd == 2 @assert length(ss) <= 2 "ss=$(ss) invalid, ss can have up to two components [ssx, ssy] for a 2D phantom" @assert length(us) <= 2 "us=$(us) invalid, us can have up to two components [usx, usy] for a 2D phantom" - if length(us) == 1 - usx = us[1]; usy = us[1] + if length(us) == 1 + usx = us[1] + usy = us[1] else - usx = us[1]; usy = us[2] + usx = us[1] + usy = us[2] end if length(ss) == 1 - ssx = ss[1]; ssy = ss[1] + ssx = ss[1] + ssy = ss[1] else - ssx = ss[1]; ssy = ss[2] + ssx = ss[1] + ssy = ss[2] end end return ssx, ssy, ssz, usx, usy, usz diff --git a/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl new file mode 100644 index 000000000..a3c3bc644 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl @@ -0,0 +1,120 @@ +# TODO: Consider different Extrapolations apart from periodic LinerInterpolator{T,ETPType} +# Interpolator{T,Degree,ETPType}, +# Degree = Linear,Cubic.... +# ETPType = Periodic, Flat... +const LinearInterpolator = Interpolations.Extrapolation{ + T, + 1, + Interpolations.GriddedInterpolation{T,1,V,Gridded{Linear{Throw{OnGrid}}},Tuple{V}}, + Gridded{Linear{Throw{OnGrid}}}, + Periodic{Nothing}, +} where {T<:Real,V<:AbstractVector{T}} + +""" +Arbitrary Motion + +x = x + ux +""" +struct ArbitraryMotion{T<:Real,V<:AbstractVector{T}} <: MotionModel{T} + period_durations::Vector{T} + dx::Array{T,2} + dy::Array{T,2} + dz::Array{T,2} + ux::Vector{LinearInterpolator{T,V}} + uy::Vector{LinearInterpolator{T,V}} + uz::Vector{LinearInterpolator{T,V}} +end + +function ArbitraryMotion( + period_durations::AbstractVector{T}, + Δx::AbstractArray{T,2}, + Δy::AbstractArray{T,2}, + Δz::AbstractArray{T,2}, +) where {T<:Real} + @warn "Note that ArbitraryMotion is under development so it is not optimized so far" maxlog = 1 + Ns = size(Δx)[1] + num_pieces = size(Δx)[2] + 1 + limits = times(period_durations, num_pieces) + + #! format: off + Δ = zeros(Ns,length(limits),4) + Δ[:,:,1] = hcat(repeat(hcat(zeros(Ns,1),Δx),1,length(period_durations)),zeros(Ns,1)) + Δ[:,:,2] = hcat(repeat(hcat(zeros(Ns,1),Δy),1,length(period_durations)),zeros(Ns,1)) + Δ[:,:,3] = hcat(repeat(hcat(zeros(Ns,1),Δz),1,length(period_durations)),zeros(Ns,1)) + + etpx = [extrapolate(interpolate((limits,), Δ[i,:,1], Gridded(Linear())), Periodic()) for i in 1:Ns] + etpy = [extrapolate(interpolate((limits,), Δ[i,:,2], Gridded(Linear())), Periodic()) for i in 1:Ns] + etpz = [extrapolate(interpolate((limits,), Δ[i,:,3], Gridded(Linear())), Periodic()) for i in 1:Ns] + #! format: on + + return ArbitraryMotion(period_durations, Δx, Δy, Δz, etpx, etpy, etpz) +end + +function Base.getindex( + motion::ArbitraryMotion, p::Union{AbstractRange,AbstractVector,Colon} +) + fields = [] + for field in fieldnames(ArbitraryMotion) + if field in (:dx, :dy, :dz) + push!(fields, getfield(motion, field)[p, :]) + elseif field in (:ux, :uy, :uz) + push!(fields, getfield(motion, field)[p]) + else + push!(fields, getfield(motion, field)) + end + end + return ArbitraryMotion(fields...) +end + +Base.:(==)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) == getfield(m2, field) for field in fieldnames(ArbitraryMotion)]) +Base.:(≈)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) ≈ getfield(m2, field) for field in fieldnames(ArbitraryMotion)]) + +function Base.vcat(m1::ArbitraryMotion, m2::ArbitraryMotion) + fields = [] + @assert m1.period_durations == m2.period_durations "period_durations of both ArbitraryMotions must be the same" + for field in + Iterators.filter(x -> !(x == :period_durations), fieldnames(ArbitraryMotion)) + push!(fields, [getfield(m1, field); getfield(m2, field)]) + end + return ArbitraryMotion(m1.period_durations, fields...) +end + +""" + limits = times(obj.motion) +""" +function times(motion::ArbitraryMotion) + period_durations = motion.period_durations + num_pieces = size(motion.dx)[2] + 1 + return times(period_durations, num_pieces) +end + +function times(period_durations::AbstractVector, num_pieces::Int) + # Pre-allocating memory + limits = zeros(eltype(period_durations), num_pieces * length(period_durations) + 1) + + idx = 1 + for i in 1:length(period_durations) + segment_increment = period_durations[i] / num_pieces + cumulative_sum = limits[idx] # Start from the last computed value in limits + for j in 1:num_pieces + cumulative_sum += segment_increment + limits[idx + 1] = cumulative_sum + idx += 1 + end + end + return limits +end + +# TODO: Calculate interpolation functions "on the fly" +function get_spin_coords( + motion::ArbitraryMotion{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + xt = x .+ reduce(vcat, [etp.(t) for etp in motion.ux]) + yt = y .+ reduce(vcat, [etp.(t) for etp in motion.uy]) + zt = z .+ reduce(vcat, [etp.(t) for etp in motion.uz]) + return xt, yt, zt +end diff --git a/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl new file mode 100644 index 000000000..3a453d98c --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl @@ -0,0 +1,28 @@ +""" +No Motion + +x = x +""" + +struct NoMotion{T<:Real} <: MotionModel{T} end + +Base.getindex(motion::NoMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion + +Base.:(==)(m1::NoMotion, m2::NoMotion) = true +Base.:(≈)(m1::NoMotion, m2::NoMotion) = true + +Base.vcat(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = NoMotion{T}() + +function get_spin_coords( + motion::NoMotion, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + return x, y, z +end + +function times(motion::NoMotion) + return [0.0] +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl new file mode 100644 index 000000000..c9387f337 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl @@ -0,0 +1,120 @@ +# ------ SimpleMotionType +abstract type SimpleMotionType{T<:Real} end + +is_composable(motion_type::SimpleMotionType{T}) where {T<:Real} = false + +""" + motion = SimpleMotion(types) + +SimpleMotion model + +# Arguments +- `types`: (`::Vector{<:SimpleMotionType{T}}`) vector of simple motion types + +# Returns +- `motion`: (`::SimpleMotion`) SimpleMotion struct +""" +struct SimpleMotion{T<:Real} <: MotionModel{T} + types::Vector{<:SimpleMotionType{T}} +end + +Base.getindex(motion::SimpleMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion +function Base.getindex( + motion::SimpleMotion, + p::Union{AbstractRange,AbstractVector,Colon}, + q::Union{AbstractRange,AbstractVector,Colon}, +) + return motion +end + +Base.vcat(m1::SimpleMotion, m2::SimpleMotion) = SimpleMotion([m1.types; m2.types]) + +Base.:(==)(m1::SimpleMotion, m2::SimpleMotion) = reduce(&, [m1.types[i] == m2.types[i] for i in 1:length(m1.types)]) +Base.:(≈)(m1::SimpleMotion, m2::SimpleMotion) = reduce(&, [m1.types[i] ≈ m2.types[i] for i in 1:length(m1.types)]) +# When they are the same type +Base.:(==)(t1::T, t2::T) where {T<:SimpleMotionType} = reduce(&, [getfield(t1, field) == getfield(t2, field) for field in fieldnames(T)]) +Base.:(≈)(t1::T, t2::T) where {T<:SimpleMotionType} = reduce(&, [getfield(t1, field) ≈ getfield(t2, field) for field in fieldnames(T)]) +# When they are not (default) +Base.:(==)(t1::SimpleMotionType, t2::SimpleMotionType) = false +Base.:(≈)(t1::SimpleMotionType, t2::SimpleMotionType) = false + +""" + x, y, z = het_spin_coords(motion, x, y, z, t') + +Calculates the position of each spin at a set of arbitrary time instants, i.e. the time steps of the simulation. + +# Arguments +- `motion`: (`::MotionModel`) phantom motion +- `x`: (`::AbstractVector{T<:Real}`, `[m]`) spin x-position vector +- `y`: (`::AbstractVector{T<:Real}`, `[m]`) spin y-position vector +- `z`: (`::AbstractVector{T<:Real}`, `[m]`) spin z-position vector +- `t`: (`::AbstractArray{T<:Real}`) horizontal array of time instants + +# Returns +- `z, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time +""" +function get_spin_coords( + motion::SimpleMotion{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + xt, yt, zt = x .+ 0*t, y .+ 0*t, z .+ 0*t + # Composable motions: they need to be run sequentially + for motion in Iterators.filter(is_composable, motion.types) + aux = xt .+ 0, yt .+ 0, zt .+ 0 + xt .+= displacement_x(motion, aux..., t) + yt .+= displacement_y(motion, aux..., t) + zt .+= displacement_z(motion, aux..., t) + end + # Additive motions: these motions can be run in parallel + for motion in Iterators.filter(!is_composable, motion.types) + xt .+= displacement_x(motion, x, y, z, t) + yt .+= displacement_y(motion, x, y, z, t) + zt .+= displacement_z(motion, x, y, z, t) + end + return xt, yt, zt +end + +function times(motion::SimpleMotion) + nodes = reduce(vcat, [times(type) for type in motion.types]) + nodes = unique(sort(nodes)) + return nodes +end + +function sort_motions!(motion::SimpleMotion) + return sort!(motion.types; by=m -> times(m)[1]) +end + +# --------- Simple Motion Types: ------------- +# Non-periodic types: defined by an initial time (t_start), an end time (t_end) and a displacement +include("simplemotion/Translation.jl") +include("simplemotion/Rotation.jl") +include("simplemotion/HeartBeat.jl") + +function unit_time(t::AbstractArray{T}, t_start::T, t_end::T) where {T<:Real} + if t_start == t_end + return (t .>= t_start) .* one(T) # The problem with this is that it returns a BitVector (type stability issues) + else + return min.(max.((t .- t_start) ./ (t_end - t_start), zero(T)), one(T)) + end +end + +# Periodic types: defined by the period, the temporal symmetry and a displacement (amplitude) +include("simplemotion/PeriodicTranslation.jl") +include("simplemotion/PeriodicRotation.jl") +include("simplemotion/PeriodicHeartBeat.jl") + +function unit_time_triangular(t, period, asymmetry) + t_rise = period * asymmetry + t_fall = period * (1 - asymmetry) + t_relative = mod.(t, period) + t_unit = + ifelse.( + t_relative .< t_rise, + t_relative ./ t_rise, + 1 .- (t_relative .- t_rise) ./ t_fall, + ) + return t_unit +end diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl new file mode 100644 index 000000000..b9e6b30b1 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl @@ -0,0 +1,81 @@ +@doc raw""" + heartbeat = HeartBeat(t_start, t_end, dx, dy, dz) + +HeartBeat struct. It produces a heartbeat-like motion, characterised by three types of strain: +Circumferential, Radial and Longitudinal + +# Arguments +- `circumferential_strain`: (`::Real`, `=-0.3`) contraction parameter +- `radial_strain`: (`::Real`, `=-0.3`) contraction parameter +- `longitudinal_strain`: (`::Real`, `=1`) contraction parameter +- `t_start`: (`::Real`, `[s]`) initial time +- `t_end`: (`::Real`, `[s]`) final time + +# Returns +- `heartbeat`: (`::HeartBeat`) HeartBeat struct +""" +@with_kw struct HeartBeat{T<:Real} <: SimpleMotionType{T} + circumferential_strain :: T + radial_strain :: T + longitudinal_strain::T = typeof(circumferential_strain)(0.0) + t_start::T = typeof(circumferential_strain)(0.0) + t_end::T = typeof(circumferential_strain)(0.0) + @assert t_end >= t_start "t_end must be major or equal than t_start" +end + +is_composable(motion_type::HeartBeat) = true + +function displacement_x( + motion_type::HeartBeat{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + r = sqrt.(x .^ 2 + y .^ 2) + θ = atan.(y, x) + Δ_circunferential = motion_type.circumferential_strain * maximum(r) + Δ_radial = -motion_type.radial_strain * (maximum(r) .- r) + Δr = t_unit .* (Δ_circunferential .+ Δ_radial) + # Map negative radius to r=0 + neg = (r .+ Δr) .< 0 + Δr = (.!neg) .* Δr + Δr .-= neg .* r + return Δr .* cos.(θ) +end + +function displacement_y( + motion_type::HeartBeat{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + r = sqrt.(x .^ 2 + y .^ 2) + θ = atan.(y, x) + Δ_circunferential = motion_type.circumferential_strain * maximum(r) + Δ_radial = -motion_type.radial_strain * (maximum(r) .- r) + Δr = t_unit .* (Δ_circunferential .+ Δ_radial) + # Map negative radius to r=0 + neg = (r .+ Δr) .< 0 + Δr = (.!neg) .* Δr + Δr .-= neg .* r + return Δr .* sin.(θ) +end + +function displacement_z( + motion_type::HeartBeat{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + return t_unit .* (z .* motion_type.longitudinal_strain) +end + +times(motion_type::HeartBeat) = begin + return [motion_type.t_start, motion_type.t_end] +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl new file mode 100644 index 000000000..17c201848 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl @@ -0,0 +1,80 @@ +@doc raw""" + periodic_heartbeat = PeriodicHeartBeat(t_start, t_end, dx, dy, dz) + +HeartBeat struct. It produces a heartbeat-like motion, characterised by three types of strain: +Circumferential, Radial and Longitudinal + +# Arguments +- `circumferential_strain`: (`::Real`, `=-0.3`) contraction parameter +- `radial_strain`: (`::Real`, `=-0.3`) contraction parameter +- `longitudinal_strain`: (`::Real`, `=1`) contraction parameter +- `period`: (`::Real`, `[s]`) period +- `asymmetry`: (`::Real`) asymmetry factor, between 0 and 1 + +# Returns +- `periodic_heartbeat`: (`::PeriodicHeartBeat`) PeriodicHeartBeat struct +""" +@with_kw struct PeriodicHeartBeat{T<:Real} <: SimpleMotionType{T} + circumferential_strain :: T + radial_strain :: T + longitudinal_strain::T = typeof(circumferential_strain)(0.0) + period::T = typeof(circumferential_strain)(0.0) + asymmetry::T = typeof(circumferential_strain)(0.5) +end + +is_composable(motion_type::PeriodicHeartBeat) = true + +function displacement_x( + motion_type::PeriodicHeartBeat{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + r = sqrt.(x .^ 2 + y .^ 2) + θ = atan.(y, x) + Δ_circunferential = motion_type.circumferential_strain * maximum(r) + Δ_radial = -motion_type.radial_strain * (maximum(r) .- r) + Δr = t_unit .* (Δ_circunferential .+ Δ_radial) + # Map negative radius to r=0 + neg = (r .+ Δr) .< 0 + Δr = (.!neg) .* Δr + Δr .-= neg .* r + return Δr .* cos.(θ) +end + +function displacement_y( + motion_type::PeriodicHeartBeat{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + r = sqrt.(x .^ 2 + y .^ 2) + θ = atan.(y, x) + Δ_circunferential = motion_type.circumferential_strain * maximum(r) + Δ_radial = -motion_type.radial_strain * (maximum(r) .- r) + Δr = t_unit .* (Δ_circunferential .+ Δ_radial) + # Map negative radius to r=0 + neg = (r .+ Δr) .< 0 + Δr = (.!neg) .* Δr + Δr .-= neg .* r + return Δr .* sin.(θ) +end + +function displacement_z( + motion_type::PeriodicHeartBeat{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + return t_unit .* (z .* motion_type.longitudinal_strain) +end + +function times(motion_type::PeriodicHeartBeat) + return [0, motion_type.period * motion_type.asymmetry, motion_type.period] +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl new file mode 100644 index 000000000..b80af6463 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl @@ -0,0 +1,76 @@ +@doc raw""" + periodic_rotation = PeriodicRotation(period, asymmetry, pitch, roll, yaw) + +PeriodicRotation motion struct. It produces a rotation of the phantom in the three axes: +x (pitch), y (roll), and z (yaw) + +# Arguments +- `pitch`: (`::Real`, `[º]`) rotation in x +- `roll`: (`::Real`, `[º]`) rotation in y +- `yaw`: (`::Real`, `[º]`) rotation in z +- `period`: (`::Real`, `[s]`) period +- `asymmetry`: (`::Real`) asymmetry factor, between 0 and 1 + +# Returns +- `periodic_rotation`: (`::PeriodicRotation`) PeriodicRotation struct + +""" +@with_kw struct PeriodicRotation{T<:Real} <: SimpleMotionType{T} + pitch :: T + roll :: T + yaw :: T + period::T = typeof(pitch)(0.0) + asymmetry::T = typeof(pitch)(0.5) +end + +is_composable(motion_type::PeriodicRotation) = true + +function displacement_x( + motion_type::PeriodicRotation{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + α = t_unit .* motion_type.pitch + β = t_unit .* motion_type.roll + γ = t_unit .* motion_type.yaw + return cosd.(γ) .* cosd.(β) .* x + + (cosd.(γ) .* sind.(β) .* sind.(α) .- sind.(γ) .* cosd.(α)) .* y + + (cosd.(γ) .* sind.(β) .* cosd.(α) .+ sind.(γ) .* sind.(α)) .* z .- x +end + +function displacement_y( + motion_type::PeriodicRotation{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + α = t_unit .* motion_type.pitch + β = t_unit .* motion_type.roll + γ = t_unit .* motion_type.yaw + return sind.(γ) .* cosd.(β) .* x + + (sind.(γ) .* sind.(β) .* sind.(α) .+ cosd.(γ) .* cosd.(α)) .* y + + (sind.(γ) .* sind.(β) .* cosd.(α) .- cosd.(γ) .* sind.(α)) .* z .- y +end + +function displacement_z( + motion_type::PeriodicRotation{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + α = t_unit .* motion_type.pitch + β = t_unit .* motion_type.roll + γ = t_unit .* motion_type.yaw + return -sind.(β) .* x + cosd.(β) .* sind.(α) .* y .- z +end + +function times(motion_type::PeriodicRotation) + return [0, motion_type.period * motion_type.asymmetry, motion_type.period] +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl new file mode 100644 index 000000000..4347ad33d --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl @@ -0,0 +1,61 @@ +@doc raw""" + periodic_translation = PeriodicTranslation(period, asymmetry, dx, dy, dz) + +PeriodicTranslation motion struct. It produces a periodic translation of the phantom in the three directions x, y and z. +The amplitude of the oscillation will be defined by dx, dy and dz + +# Arguments +- `dx`: (`::Real`, `[m]`) translation in x +- `dy`: (`::Real`, `[m]`) translation in y +- `dz`: (`::Real`, `[m]`) translation in z +- `period`: (`::Real`, `[s]`) period +- `asymmetry`: (`::Real`) asymmetry factor, between 0 and 1 + +# Returns +- `periodic_translation`: (`::PeriodicTranslation`) PeriodicTranslation struct + +""" +@with_kw struct PeriodicTranslation{T<:Real} <: SimpleMotionType{T} + dx :: T + dy :: T + dz :: T + period::T = typeof(dx)(0.0) + asymmetry::T = typeof(dx)(0.5) +end + +function displacement_x( + motion_type::PeriodicTranslation{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + return t_unit .* motion_type.dx +end + +function displacement_y( + motion_type::PeriodicTranslation{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + return t_unit .* motion_type.dy +end + +function displacement_z( + motion_type::PeriodicTranslation{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) + return t_unit .* motion_type.dz +end + +function times(motion_type::PeriodicTranslation) + return [0, motion_type.period * motion_type.asymmetry, motion_type.period] +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl new file mode 100644 index 000000000..a81920d6d --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl @@ -0,0 +1,126 @@ +@doc raw""" + rotation = Rotation(t_start, t_end, pitch, roll, yaw) + +Rotation motion struct. It produces a rotation of the phantom in the three axes: +x (pitch), y (roll), and z (yaw) + +```math +\begin{align*} +ux &= cos(\gamma)cos(\beta)x \\ + &+ (cos(\gamma)sin(\beta)sin(\alpha) - sin(\gamma)cos(\alpha))y\\ + &+ (cos(\gamma)sin(\beta)cos(\alpha) + sin(\gamma)sin(\alpha))z\\ + &- x +\end{align*} +``` + +```math +\begin{align*} +uy &= sin(\gamma)cos(\beta)x \\ + &+ (sin(\gamma)sin(\beta)sin(\alpha) + cos(\gamma)cos(\alpha))y\\ + &+ (sin(\gamma)sin(\beta)cos(\alpha) - cos(\gamma)sin(\alpha))z\\ + &- y +\end{align*} +``` + +```math +\begin{align*} +uz &= -sin(\beta)x \\ + &+ cos(\beta)sin(\alpha)y\\ + &+ cos(\beta)cos(\alpha)z\\ + &- z +\end{align*} +``` + +where: + +```math +\alpha = \left\{\begin{matrix} +0, & t <= t_start \\ +\frac{pitch}{t_end-t_start}(t-t_start), & t_start < t < t_end \\ +pitch, & t >= t_end +\end{matrix}\right. +,\qquad +\beta = \left\{\begin{matrix} +0, & t <= t_start \\ +\frac{roll}{t_end-t_start}(t-t_start), & t_start < t < t_end \\ +roll, & t >= t_end +\end{matrix}\right. +,\qquad +\gamma = \left\{\begin{matrix} +0, & t <= t_start \\ +\frac{yaw}{t_end-t_start}(t-t_start), & t_start < t < t_end \\ +yaw, & t >= t_end +\end{matrix}\right. +``` + +# Arguments +- `pitch`: (`::Real`, `[º]`) rotation in x +- `roll`: (`::Real`, `[º]`) rotation in y +- `yaw`: (`::Real`, `[º]`) rotation in z +- `t_start`: (`::Real`, `[s]`) initial time +- `t_end`: (`::Real`, `[s]`) final time + +# Returns +- `rotation`: (`::Rotation`) Rotation struct + +""" +@with_kw struct Rotation{T<:Real} <: SimpleMotionType{T} + pitch :: T + roll :: T + yaw :: T + t_start::T = typeof(pitch)(0.0) + t_end = typeof(pitch)(0.0) + @assert t_end >= t_start "t_end must be major or equal than t_start" +end + +is_composable(motion_type::Rotation) = true + +function displacement_x( + motion_type::Rotation{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + α = t_unit .* (motion_type.pitch) + β = t_unit .* (motion_type.roll) + γ = t_unit .* (motion_type.yaw) + return cosd.(γ) .* cosd.(β) .* x + + (cosd.(γ) .* sind.(β) .* sind.(α) .- sind.(γ) .* cosd.(α)) .* y + + (cosd.(γ) .* sind.(β) .* cosd.(α) .+ sind.(γ) .* sind.(α)) .* z .- x +end + +function displacement_y( + motion_type::Rotation{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + α = t_unit .* motion_type.pitch + β = t_unit .* motion_type.roll + γ = t_unit .* motion_type.yaw + return sind.(γ) .* cosd.(β) .* x + + (sind.(γ) .* sind.(β) .* sind.(α) .+ cosd.(γ) .* cosd.(α)) .* y + + (sind.(γ) .* sind.(β) .* cosd.(α) .- cosd.(γ) .* sind.(α)) .* z .- y +end + +function displacement_z( + motion_type::Rotation{T}, + x::AbstractArray{T}, + y::AbstractArray{T}, + z::AbstractArray{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + α = t_unit .* motion_type.pitch + β = t_unit .* motion_type.roll + γ = t_unit .* motion_type.yaw + return -sind.(β) .* x + cosd.(β) .* sind.(α) .* y .- z +end + +times(motion_type::Rotation) = begin + return [motion_type.t_start, motion_type.t_end] +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl new file mode 100644 index 000000000..c002be489 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl @@ -0,0 +1,85 @@ +@doc raw""" + translation = Translation(t_start, t_end, dx, dy, dz) + +Translation motion struct. It produces a translation of the phantom in the three directions x, y and z. + +```math +ux=\left\{\begin{matrix} +0, & t <= t_start\\ +\frac{dx}{t_end-t_start}(t-t_start), & t_start < t < t_end\\ +dx, & t >= t_end +\end{matrix}\right. +``` + +```math +uy=\left\{\begin{matrix} +0, & t <= t_start\\ +\frac{dy}{t_end-t_start}(t-t_start), & t_start < t < t_end\\ +dy, & t >= t_end +\end{matrix}\right. +``` + +```math +uz=\left\{\begin{matrix} +0, & t <= t_start\\ +\frac{dz}{t_end-t_start}(t-t_start), & t_start < t < t_end\\ +dz, & t >= t_end +\end{matrix}\right. +``` + +# Arguments +- `dx`: (`::Real`, `[m]`) translation in x +- `dy`: (`::Real`, `[m]`) translation in y +- `dz`: (`::Real`, `[m]`) translation in z +- `t_start`: (`::Real`, `[s]`) initial time +- `t_end`: (`::Real`, `[s]`) final time + +# Returns +- `translation`: (`::Translation`) Translation struct + +""" +@with_kw struct Translation{T<:Real} <: SimpleMotionType{T} + dx :: T + dy :: T + dz :: T + t_start::T = typeof(dx)(0.0) + t_end::T = typeof(dx)(0.0) + @assert t_end >= t_start "t_end must be major or equal than t_start" +end + +function displacement_x( + motion_type::Translation{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + return t_unit .* motion_type.dx +end + +function displacement_y( + motion_type::Translation{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + return t_unit .* motion_type.dy +end + +function displacement_z( + motion_type::Translation{T}, + x::AbstractVector{T}, + y::AbstractVector{T}, + z::AbstractVector{T}, + t::AbstractArray{T}, +) where {T<:Real} + t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) + return t_unit .* motion_type.dz +end + +times(motion_type::Translation) = begin + return [motion_type.t_start, motion_type.t_end] +end diff --git a/KomaMRIBase/src/datatypes/sequence/Grad.jl b/KomaMRIBase/src/datatypes/sequence/Grad.jl index 600ef2c83..456d6694f 100644 --- a/KomaMRIBase/src/datatypes/sequence/Grad.jl +++ b/KomaMRIBase/src/datatypes/sequence/Grad.jl @@ -9,9 +9,11 @@ Rotates vector counter-clockwise with respect to the x-axis. # Returns - `Rx`: (`::Matrix{Int64}`) rotation matrix """ -rotx(θ::Real) = [1 0 0 - 0 cos(θ) -sin(θ); - 0 sin(θ) cos(θ)] +rotx(θ::Real) = [ + 1 0 0 + 0 cos(θ) -sin(θ) + 0 sin(θ) cos(θ) +] """ Ry = roty(θ::Real) @@ -24,9 +26,11 @@ Rotates vector counter-clockwise with respect to the y-axis. # Returns - `Ry`: (`::Matrix{Int64}`) rotation matrix """ -roty(θ::Real) = [cos(θ) 0 sin(θ); - 0 1 0; - -sin(θ) 0 cos(θ)] +roty(θ::Real) = [ + cos(θ) 0 sin(θ) + 0 1 0 + -sin(θ) 0 cos(θ) +] """ Rz = rotz(θ::Real) @@ -39,9 +43,11 @@ Rotates vector counter-clockwise with respect to the z-axis. # Returns - `Rz`: (`::Matrix{Int64}`) rotation matrix """ -rotz(θ::Real) = [cos(θ) -sin(θ) 0; - sin(θ) cos(θ) 0; - 0 0 1] +rotz(θ::Real) = [ + cos(θ) -sin(θ) 0 + sin(θ) cos(θ) 0 + 0 0 1 +] """ gr = Grad(A, T) @@ -70,28 +76,43 @@ julia> seq = Sequence([gr]); plot_seq(seq) ``` """ mutable struct Grad - A - T - rise::Real - fall::Real - delay::Real - first - last + A + T + rise::Real + fall::Real + delay::Real + first + last function Grad(A, T, rise, fall, delay) - all(T .< 0) || rise < 0 || fall < 0 || delay < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, fall, delay, 0.0, 0.0) + return if all(T .< 0) || rise < 0 || fall < 0 || delay < 0 + error("Gradient timings must be positive.") + else + new(A, T, rise, fall, delay, 0.0, 0.0) + end end - function Grad(A, T, rise, delay) - all(T .< 0) < 0 || rise < 0 || delay < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, rise, delay, 0.0, 0.0) + function Grad(A, T, rise, delay) + return if all(T .< 0) < 0 || rise < 0 || delay < 0 + error("Gradient timings must be positive.") + else + new(A, T, rise, rise, delay, 0.0, 0.0) + end end - function Grad(A, T, rise) - all(T .< 0) < 0 || rise < 0 ? error("Gradient timings must be positive.") : new(A, T, rise, rise, 0.0, 0.0, 0.0) + function Grad(A, T, rise) + return if all(T .< 0) < 0 || rise < 0 + error("Gradient timings must be positive.") + else + new(A, T, rise, rise, 0.0, 0.0, 0.0) + end end - function Grad(A, T) - all(T .< 0) < 0 ? error("Gradient timings must be positive.") : new(A, T, 0.0, 0.0, 0.0, 0.0, 0.0) + function Grad(A, T) + return if all(T .< 0) < 0 + error("Gradient timings must be positive.") + else + new(A, T, 0.0, 0.0, 0.0, 0.0, 0.0) + end end end - """ gr = Grad(f::Function, T::Real, N::Integer; delay::Real) @@ -117,12 +138,11 @@ julia> seq = Sequence([gx]); plot_seq(seq) ``` """ Grad(f::Function, T::Real, N::Integer=300; delay::Real=0) = begin - t = range(0.0, T; length=N) - G = f.(t) - return Grad(G, T, 0.0, 0.0, delay) + t = range(0.0, T; length=N) + G = f.(t) + return Grad(G, T, 0.0, 0.0, delay) end - """ str = show(io::IO, x::Grad) @@ -135,19 +155,30 @@ Displays information about the Grad struct `x` in the julia REPL. - `str` (`::String`) output string message """ Base.show(io::IO, x::Grad) = begin - r(x) = round.(x,digits=4) - compact = get(io, :compact, false) - if !compact - wave = length(x.A) == 1 ? r(x.A*1e3) : "∿" - if x.rise == x.fall == 0. - print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(sum(x.T)*1e3)) ms)") - else - print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"Grad($(wave) mT, $(r(sum(x.T)*1e3)) ms, ↑$(r(x.rise*1e3)) ms, ↓$(r(x.fall*1e3)) ms)") - end - else - wave = length(x.A) == 1 ? "⊓" : "∿" - print(io, (is_on(x) ? wave : "⇿")*"($(r((x.delay+x.rise+x.fall+sum(x.T))*1e3)) ms)") - end + r(x) = round.(x, digits=4) + compact = get(io, :compact, false) + if !compact + wave = length(x.A) == 1 ? r(x.A * 1e3) : "∿" + if x.rise == x.fall == 0.0 + print( + io, + (x.delay > 0 ? "←$(r(x.delay*1e3)) ms→ " : "") * + "Grad($(wave) mT, $(r(sum(x.T)*1e3)) ms)", + ) + else + print( + io, + (x.delay > 0 ? "←$(r(x.delay*1e3)) ms→ " : "") * + "Grad($(wave) mT, $(r(sum(x.T)*1e3)) ms, ↑$(r(x.rise*1e3)) ms, ↓$(r(x.fall*1e3)) ms)", + ) + end + else + wave = length(x.A) == 1 ? "⊓" : "∿" + print( + io, + (is_on(x) ? wave : "⇿") * "($(r((x.delay+x.rise+x.fall+sum(x.T))*1e3)) ms)", + ) + end end """ @@ -167,45 +198,58 @@ directly without the need to iterate elementwise. - `y`: (`::Vector{Any}` or `::Matrix{Any}`) vector or matrix with the property defined by the symbol `f` for all elements of the Grad vector or matrix `x` """ -getproperty(x::Vector{Grad}, f::Symbol) = getproperty.(x,f) +getproperty(x::Vector{Grad}, f::Symbol) = getproperty.(x, f) getproperty(x::Matrix{Grad}, f::Symbol) = begin - if f == :x - x[1,:] - elseif f == :y && size(x,1) >= 2 - x[2,:] - elseif f == :z && size(x,1) >= 3 - x[3,:] - elseif f == :dur - maximum(dur.(x), dims=1)[:] - else - getproperty.(x,f) - end + if f == :x + x[1, :] + elseif f == :y && size(x, 1) >= 2 + x[2, :] + elseif f == :z && size(x, 1) >= 3 + x[3, :] + elseif f == :dur + maximum(dur.(x); dims=1)[:] + else + getproperty.(x, f) + end end # Gradient comparison -Base.isapprox(gr1::Grad, gr2::Grad) = begin - return all(length(getfield(gr1, k)) ≈ length(getfield(gr2, k)) for k ∈ fieldnames(Grad)) && - all(getfield(gr1, k) ≈ getfield(gr2, k) for k ∈ fieldnames(Grad)) +function Base.isapprox(gr1::Grad, gr2::Grad) + return all( + length(getfield(gr1, k)) ≈ length(getfield(gr2, k)) for k in fieldnames(Grad) + ) && all(getfield(gr1, k) ≈ getfield(gr2, k) for k in fieldnames(Grad)) end # Gradient operations -*(x::Grad,α::Real) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay) -*(α::Real,x::Grad) = Grad(α*x.A,x.T,x.rise,x.fall,x.delay) -*(A::Matrix{Float64},GR::Matrix{Grad}) = begin - N, M = size(GR) - [sum(A[i,1:N] .* GR[:,j]) for i=1:N, j=1:M] +*(x::Grad, α::Real) = Grad(α * x.A, x.T, x.rise, x.fall, x.delay) +*(α::Real, x::Grad) = Grad(α * x.A, x.T, x.rise, x.fall, x.delay) +function *(A::Matrix{Float64}, GR::Matrix{Grad}) + N, M = size(GR) + return [sum(A[i, 1:N] .* GR[:, j]) for i in 1:N, j in 1:M] end -Base.zero(::Grad) = Grad(0.0,0.0) +Base.zero(::Grad) = Grad(0.0, 0.0) size(g::Grad, i::Int64) = 1 #To fix [g;g;;] concatenation of Julia 1.7.3 -/(x::Grad,α::Real) = Grad(x.A/α,x.T,x.rise,x.fall,x.delay) -+(x::Grad,y::Grad) = Grad(x.A.+y.A,x.T,x.rise,x.fall,x.delay) -+(x::Array{Grad,1}, y::Array{Grad,1}) = [x[i]+y[i] for i=1:length(x)] --(x::Grad) = -1*x --(x::Grad,y::Grad) = Grad(x.A.-y.A,x.T,x.rise,x.fall,x.delay) +/(x::Grad, α::Real) = Grad(x.A / α, x.T, x.rise, x.fall, x.delay) ++(x::Grad, y::Grad) = Grad(x.A .+ y.A, x.T, x.rise, x.fall, x.delay) ++(x::Array{Grad,1}, y::Array{Grad,1}) = [x[i] + y[i] for i in 1:length(x)] +-(x::Grad) = -1 * x +-(x::Grad, y::Grad) = Grad(x.A .- y.A, x.T, x.rise, x.fall, x.delay) # Gradient functions -vcat(x::Array{Grad,1},y::Array{Grad,1}) = [i==1 ? x[j] : y[j] for i=1:2,j=1:length(x)] -vcat(x::Array{Grad,1},y::Array{Grad,1},z::Array{Grad,1}) = [i==1 ? x[j] : i==2 ? y[j] : z[j] for i=1:3,j=1:length(x)] +function vcat(x::Array{Grad,1}, y::Array{Grad,1}) + return [i == 1 ? x[j] : y[j] for i in 1:2, j in 1:length(x)] +end +function vcat(x::Array{Grad,1}, y::Array{Grad,1}, z::Array{Grad,1}) + return [ + if i == 1 + x[j] + elseif i == 2 + y[j] + else + z[j] + end for i in 1:3, j in 1:length(x) + ] +end """ y = dur(x::Grad) @@ -221,4 +265,4 @@ the duration is the maximum duration of all the elements of the gradient vector. - `y`: (`::Float64`, `[s]`) duration of the RF struct or RF array """ dur(x::Grad) = x.delay + x.rise + sum(x.T) + x.fall -dur(x::Vector{Grad}) = maximum(dur.(x), dims=1)[:] +dur(x::Vector{Grad}) = maximum(dur.(x); dims=1)[:] \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/sequence/RF.jl b/KomaMRIBase/src/datatypes/sequence/RF.jl index 7d2268505..9d303e5bb 100644 --- a/KomaMRIBase/src/datatypes/sequence/RF.jl +++ b/KomaMRIBase/src/datatypes/sequence/RF.jl @@ -24,18 +24,22 @@ julia> seq = Sequence(); seq += rf; plot_seq(seq) ``` """ mutable struct RF - A - T - Δf - delay::Real - function RF(A, T, Δf, delay) - any(T .< 0) || delay < 0 ? error("RF timings must be non-negative.") : new(A, T, Δf, delay) + A + T + Δf + delay::Real + function RF(A, T, Δf, delay) + return if any(T .< 0) || delay < 0 + error("RF timings must be non-negative.") + else + new(A, T, Δf, delay) + end end - function RF(A, T, Δf) - any(T .< 0) ? error("RF timings must be non-negative.") : new(A, T, Δf, 0.) + function RF(A, T, Δf) + return any(T .< 0) ? error("RF timings must be non-negative.") : new(A, T, Δf, 0.0) end - function RF(A, T) - any(T .< 0) ? error("RF timings must be non-negative.") : new(A, T, 0., 0.) + function RF(A, T) + return any(T .< 0) ? error("RF timings must be non-negative.") : new(A, T, 0.0, 0.0) end end @@ -51,15 +55,19 @@ Displays information about the RF struct `x` in the julia REPL. - `str`: (`::String`) output string message """ Base.show(io::IO, x::RF) = begin - r(x) = round.(x,digits=4) - compact = get(io, :compact, false) - if !compact - wave = length(x.A) == 1 ? r(x.A*1e6) : "∿" - print(io, (x.delay>0 ? "←$(r(x.delay*1e3)) ms→ " : "")*"RF($(wave) uT, $(r(sum(x.T)*1e3)) ms, $(r(x.Δf)) Hz)") - else - wave = length(x.A) == 1 ? "⊓" : "∿" - print(io, (sum(abs.(x.A)) > 0 ? wave : "⇿")*"($(r((x.delay+sum(x.T))*1e3)) ms)") - end + r(x) = round.(x, digits=4) + compact = get(io, :compact, false) + if !compact + wave = length(x.A) == 1 ? r(x.A * 1e6) : "∿" + print( + io, + (x.delay > 0 ? "←$(r(x.delay*1e3)) ms→ " : "") * + "RF($(wave) uT, $(r(sum(x.T)*1e3)) ms, $(r(x.Δf)) Hz)", + ) + else + wave = length(x.A) == 1 ? "⊓" : "∿" + print(io, (sum(abs.(x.A)) > 0 ? wave : "⇿") * "($(r((x.delay+sum(x.T))*1e3)) ms)") + end end """ @@ -78,34 +86,35 @@ directly without the need to iterate elementwise. - `y`: (`::Vector{Any}` or `::Matrix{Any}`) vector with the property defined by the symbol `f` for all elements of the RF vector or matrix `x` """ -getproperty(x::Vector{RF}, f::Symbol) = getproperty.(x,f) +getproperty(x::Vector{RF}, f::Symbol) = getproperty.(x, f) getproperty(x::Matrix{RF}, f::Symbol) = begin - if f == :Bx - real.(getproperty.(x,:A)) - elseif f == :By - imag.(getproperty.(x,:A)) - elseif f == :Δf - getproperty.(x,:Δf) - elseif f == :T || f == :delay - getproperty.(x[1,:],f) - elseif f == :dur - T, delay = x.T, x.delay - ΔT = [sum(t) for t=T] .+ delay - ΔT - else - getproperty.(x,f) - end + if f == :Bx + real.(getproperty.(x, :A)) + elseif f == :By + imag.(getproperty.(x, :A)) + elseif f == :Δf + getproperty.(x, :Δf) + elseif f == :T || f == :delay + getproperty.(x[1, :], f) + elseif f == :dur + T, delay = x.T, x.delay + ΔT = [sum(t) for t in T] .+ delay + ΔT + else + getproperty.(x, f) + end end # RF comparison -Base.isapprox(rf1::RF, rf2::RF) = begin - return all(length(getfield(rf1, k)) == length(getfield(rf2, k)) for k ∈ fieldnames(RF)) - all(≈(getfield(rf1, k), getfield(rf2, k), atol=1e-9) for k ∈ fieldnames(RF)) +function Base.isapprox(rf1::RF, rf2::RF) + return all(length(getfield(rf1, k)) == length(getfield(rf2, k)) for k in fieldnames(RF)) + return all(≈(getfield(rf1, k), getfield(rf2, k); atol=1e-9) for k in fieldnames(RF)) end # Properties size(r::RF, i::Int64) = 1 #To fix [r;r;;] concatenation of Julia 1.7.3 -*(α::Complex{T}, x::RF) where {T<:Real} = RF(α*x.A,x.T,x.Δf,x.delay) +*(α::Complex{T}, x::RF) where {T<:Real} = RF(α * x.A, x.T, x.Δf, x.delay) +*(α::Real, x::RF) = RF(α * x.A, x.T, x.Δf, x.delay) """ y = dur(x::RF) @@ -121,8 +130,10 @@ Duration time in [s] of RF struct or RF array. - `y`: (`::Float64`, [`s`]) duration of the RF struct or RF array """ dur(x::RF) = sum(x.T) -dur(x::Array{RF,1}) = sum(sum(x[i].T) for i=1:size(x,1)) -dur(x::Array{RF,2}) = maximum(sum([sum(x[i,j].T) for i=1:size(x,1),j=1:size(x,2)],dims=2)) +dur(x::Array{RF,1}) = sum(sum(x[i].T) for i in 1:size(x, 1)) +function dur(x::Array{RF,2}) + return maximum(sum([sum(x[i, j].T) for i in 1:size(x, 1), j in 1:size(x, 2)]; dims=2)) +end """ rf = RF_fun(f::Function, T::Real, N::Int64) @@ -141,9 +152,9 @@ Generate an RF sequence with amplitudes sampled from a function waveform. - `rf`:(`::RF`) RF struct with amplitud defined by the function `f` """ RF(f::Function, T::Real, N::Int64=301; delay::Real=0, Δf=0) = begin - t = range(0,T;length=N) - A = f.(t) - RF(A,T,Δf,delay) + t = range(0, T; length=N) + A = f.(t) + RF(A, T, Δf, delay) end """ diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl index 38368f46a..f15bfaca2 100644 --- a/KomaMRIBase/test/runtests.jl +++ b/KomaMRIBase/test/runtests.jl @@ -1,6 +1,6 @@ using TestItems, TestItemRunner -@run_package_tests filter=ti->!(:skipci in ti.tags)&&(:base in ti.tags) #verbose=true +@run_package_tests filter=t_start->!(:skipci in t_start.tags)&&(:base in t_start.tags) #verbose=true @testitem "Sequence" tags=[:base] begin @testset "Init" begin @@ -359,66 +359,188 @@ end end end -@testitem "Phantom" tags=[:base] begin - +@testitem "Phantom" tags = [:base] begin + using Suppressor # Test phantom struct creation name = "Bulks" - x = [-2e-3; -1e-3; 0.; 1e-3; 2e-3] - y = [-4e-3; -2e-3; 0.; 2e-3; 4e-3] - z = [-6e-3; -3e-3; 0.; 3e-3; 6e-3] - ρ = [.2; .4; .6; .8; 1.] - T1 = [.9; .9; .5; .25; .4] - T2 = [.09; .05; .04; .07; .005] - T2s = [.1; .06; .05; .08; .015] - Δw = [-2e-6; -1e-6; 0.; 1e-6; 2e-6] - Dλ1 = [-4e-6; -2e-6; 0.; 2e-6; 4e-6] - Dλ2 = [-6e-6; -3e-6; 0.; 3e-6; 6e-6] - Dθ = [-8e-6; -4e-6; 0.; 4e-6; 8e-6] - u = (x,y,z,t)->0 - obj = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ, ux=u, uy=u, uz=u) - obj2 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ, ux=u, uy=u, uz=u) - @test obj ≈ obj2 + x = [-2e-3; -1e-3; 0.0; 1e-3; 2e-3] + y = [-4e-3; -2e-3; 0.0; 2e-3; 4e-3] + z = [-6e-3; -3e-3; 0.0; 3e-3; 6e-3] + ρ = [0.2; 0.4; 0.6; 0.8; 1.0] + T1 = [0.9; 0.9; 0.5; 0.25; 0.4] + T2 = [0.09; 0.05; 0.04; 0.07; 0.005] + T2s = [0.1; 0.06; 0.05; 0.08; 0.015] + Δw = [-2e-6; -1e-6; 0.0; 1e-6; 2e-6] + Dλ1 = [-4e-6; -2e-6; 0.0; 2e-6; 4e-6] + Dλ2 = [-6e-6; -3e-6; 0.0; 3e-6; 6e-6] + Dθ = [-8e-6; -4e-6; 0.0; 4e-6; 8e-6] + obj1 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ) + obj2 = Phantom(name=name, x=x, y=y, z=z, ρ=ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ) + @test obj1 == obj2 # Test size and length definitions of a phantom - @test size(obj) == size(ρ) - #@test length(obj) == length(ρ) + @test size(obj1) == size(ρ) + @test length(obj1) == length(ρ) + + # Test obtaining spin psositions + @testset "SimpleMotion" begin + ph = Phantom(x=[1.0], y=[1.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) + period = 2.0 + asymmetry = 0.5 + # Translation + dx, dy, dz = [1.0, 0.0, 0.0] + vx, vy, vz = [dx, dy, dz] ./ (t_end - t_start) + translation = SimpleMotion([Translation(dx, dy, dz, t_start, t_end)]) + xt, yt, zt = get_spin_coords(translation, ph.x, ph.y, ph.z, t') + @test xt == ph.x .+ vx.*t' + @test yt == ph.y .+ vy.*t' + @test zt == ph.z .+ vz.*t' + # PeriodicTranslation + periodictranslation = SimpleMotion([PeriodicTranslation(dx, dy, dz, period, asymmetry)]) + xt, yt, zt = get_spin_coords(periodictranslation, ph.x, ph.y, ph.z, t') + @test xt == ph.x .+ vx.*t' + @test yt == ph.y .+ vy.*t' + @test zt == ph.z .+ vz.*t' + # Rotation (2D) + pitch = 0.0 + roll = 0.0 + yaw = 45.0 + rotation = SimpleMotion([Rotation(pitch, roll, yaw, t_start, t_end)]) + xt, yt, zt = get_spin_coords(rotation, ph.x, ph.y, ph.z, t') + @test xt[:,end] == ph.x .* cosd(yaw) - ph.y .* sind(yaw) + @test yt[:,end] == ph.x .* sind(yaw) + ph.y .* cosd(yaw) + @test zt[:,end] == ph.z + # PeriodicRotation (2D) + periodicrotation = SimpleMotion([PeriodicRotation(pitch, roll, yaw, period, asymmetry)]) + xt, yt, zt = get_spin_coords(periodicrotation, ph.x, ph.y, ph.z, t') + @test xt[:,end] == ph.x .* cosd(yaw) - ph.y .* sind(yaw) + @test yt[:,end] == ph.x .* sind(yaw) + ph.y .* cosd(yaw) + @test zt[:,end] == ph.z + # HeartBeat + circumferential_strain = -0.1 + radial_strain = 0.0 + longitudinal_strain = -0.1 + heartbeat = SimpleMotion([HeartBeat(circumferential_strain, radial_strain, longitudinal_strain, t_start, t_end)]) + xt, yt, zt = get_spin_coords(heartbeat, ph.x, ph.y, ph.z, t') + r = sqrt.(ph.x .^ 2 + ph.y .^ 2) + θ = atan.(ph.y, ph.x) + @test xt[:,end] == ph.x .* (1 .+ circumferential_strain * maximum(r) .* cos.(θ)) + @test yt[:,end] == ph.y .* (1 .+ circumferential_strain * maximum(r) .* sin.(θ)) + @test zt[:,end] == ph.z .* (1 .+ longitudinal_strain) + # PeriodicHeartBeat + periodicheartbeat = SimpleMotion([PeriodicHeartBeat(circumferential_strain, radial_strain, longitudinal_strain, period, asymmetry)]) + xt, yt, zt = get_spin_coords(heartbeat, ph.x, ph.y, ph.z, t') + @test xt[:,end] == ph.x .* (1 .+ circumferential_strain * maximum(r) .* cos.(θ)) + @test yt[:,end] == ph.y .* (1 .+ circumferential_strain * maximum(r) .* sin.(θ)) + @test zt[:,end] == ph.z .* (1 .+ longitudinal_strain) + end + @testset "ArbitraryMotion" begin + ph = Phantom(x=[1.0], y=[1.0]) + Ns = length(ph) + period_durations = [1.0] + num_pieces = 10 + dx = dy = dz = rand(Ns, num_pieces - 1) + arbitrarymotion = @suppress ArbitraryMotion(period_durations, dx, dy, dz) + t = times(arbitrarymotion) + xt, yt, zt = get_spin_coords(arbitrarymotion, ph.x, ph.y, ph.z, t') + @test xt[:,2:end-1] == ph.x .+ dx + @test yt[:,2:end-1] == ph.y .+ dy + @test zt[:,2:end-1] == ph.z .+ dz + end + + simplemotion = SimpleMotion([ + PeriodicTranslation(dx=0.05, dy=0.05, dz=0.0, period=0.5, asymmetry=0.5), + Rotation(pitch=0.0, roll=0.0, yaw=π / 2, t_start=0.05, t_end=0.5), + ]) - # Test phantom comparison - ue = (x,y,z,t)->1 - obe = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, ue, ue, ue) - @test obj ≈ obe + Ns = length(obj1) + K = 10 + arbitrarymotion = @suppress ArbitraryMotion([1.0], 0.01 .* rand(Ns, K - 1), 0.01 .* rand(Ns, K - 1), 0.01 .* rand(Ns, K - 1)) # Test phantom subset + obs1 = Phantom( + name, + x, + y, + z, + ρ, + T1, + T2, + T2s, + Δw, + Dλ1, + Dλ2, + Dθ, + simplemotion + ) rng = 1:2:5 - ug = (x,y,z,t)->-1 - obg = Phantom(name, x[rng], y[rng], z[rng], ρ[rng], T1[rng], T2[rng], T2s[rng], Δw[rng], - Dλ1[rng], Dλ2[rng], Dθ[rng], ug, ug, ug) - @test obj[rng] ≈ obg - @test @view(obj[rng]) ≈ obg + obs2 = Phantom( + name, + x[rng], + y[rng], + z[rng], + ρ[rng], + T1[rng], + T2[rng], + T2s[rng], + Δw[rng], + Dλ1[rng], + Dλ2[rng], + Dθ[rng], + simplemotion[rng], + ) + @test obs1[rng] == obs2 + @test @view(obs1[rng]) == obs2 + + obs1.motion = arbitrarymotion + obs2.motion = arbitrarymotion[rng] + @test obs1[rng] == obs2 + # @test @view(obs1[rng]) == obs2 # Test addition of phantoms - ua = (x,y,z,t)->2 - oba = Phantom(name, [x; x[rng]], [y; y[rng]], [z; z[rng]], [ρ; ρ[rng]], - [T1; T1[rng]], [T2; T2[rng]], [T2s; T2s[rng]], [Δw; Δw[rng]], - [Dλ1; Dλ1[rng]], [Dλ2; Dλ2[rng]], [Dθ; Dθ[rng]], ua, ua, ua) - @test obj + obg ≈ oba + oba = Phantom( + name, + [x; x[rng]], + [y; y[rng]], + [z; z[rng]], + [ρ; ρ[rng]], + [T1; T1[rng]], + [T2; T2[rng]], + [T2s; T2s[rng]], + [Δw; Δw[rng]], + [Dλ1; Dλ1[rng]], + [Dλ2; Dλ2[rng]], + [Dθ; Dθ[rng]], + [obs1.motion; obs2.motion] + ) + @test obs1 + obs2 == oba # Test scalar multiplication of a phantom - c = 7. - obc = Phantom(name, x, y, z, c*ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, u, u, u) - @test c*obj ≈ obc + c = 7 + obc = Phantom(name=name, x=x, y=y, z=z, ρ=c*ρ, T1=T1, T2=T2, T2s=T2s, Δw=Δw, Dλ1=Dλ1, Dλ2=Dλ2, Dθ=Dθ) + @test c * obj1 == obc #Test brain phantom 2D ph = brain_phantom2D() - @test ph.name=="brain2D_axial" + @test ph.name == "brain2D_axial" + @test KomaMRIBase.get_dims(ph) == Bool[1, 1, 0] #Test brain phantom 3D ph = brain_phantom3D() - @test ph.name=="brain3D" + @test ph.name == "brain3D" + @test KomaMRIBase.get_dims(ph) == Bool[1, 1, 1] #Test pelvis phantom 2D ph = pelvis_phantom2D() - @test ph.name=="pelvis2D" + @test ph.name == "pelvis2D" + @test KomaMRIBase.get_dims(ph) == Bool[1, 1, 0] + + #Test heart phantom + ph = heart_phantom() + @test ph.name == "LeftVentricle" + @test KomaMRIBase.get_dims(ph) == Bool[1, 1, 0] end @testitem "Scanner" tags=[:base] begin @@ -432,8 +554,6 @@ end @testitem "TrapezoidalIntegration" tags=[:base] begin dt = Float64[1 1 1 1] x = Float64[0 1 2 1 0] - @test KomaMRIBase.trapz(dt, x)[1] ≈ 4 #Triangle area = bh/2, with b = 4 and h = 2 - @test KomaMRIBase.cumtrapz(dt, x) ≈ [0.5 2 3.5 4] end diff --git a/KomaMRICore/src/other/DiffusionModel.jl b/KomaMRICore/src/other/DiffusionModel.jl index e4fbb2b74..d993840c1 100644 --- a/KomaMRICore/src/other/DiffusionModel.jl +++ b/KomaMRICore/src/other/DiffusionModel.jl @@ -1,10 +1,12 @@ ## UNDER CONSTRUCTION! """Matrix A such as A ⋅ b = n × b""" cross(n) = begin - nx,ny,nz = n - [0 -nz ny; - nz 0 -nx; - -ny nx 0] + nx, ny, nz = n + [ + 0 -nz ny + nz 0 -nx + -ny nx 0 + ] end """Rodrigues' formula: Rotation matrix that when applied rotates with respect to "n" in an angle θ anti clock-wise""" -Un(θ,n) = [1 0 0; 0 1 0; 0 0 1] * cos(θ) + sin(θ) * cross(n) + (1-cos(θ)) * (n * n') +Un(θ, n) = [1 0 0; 0 1 0; 0 0 1] * cos(θ) + sin(θ) * cross(n) + (1 - cos(θ)) * (n * n') diff --git a/KomaMRICore/src/simulation/Bloch/BlochDictSimulationMethod.jl b/KomaMRICore/src/simulation/Bloch/BlochDictSimulationMethod.jl index 27ac99ec4..11d0f57c6 100644 --- a/KomaMRICore/src/simulation/Bloch/BlochDictSimulationMethod.jl +++ b/KomaMRICore/src/simulation/Bloch/BlochDictSimulationMethod.jl @@ -1,13 +1,15 @@ Base.@kwdef struct BlochDict <: SimulationMethod - save_Mz::Bool=false + save_Mz::Bool = false end export BlochDict Base.show(io::IO, s::BlochDict) = begin - print(io, "BlochDict(save_Mz=$(s.save_Mz))") + print(io, "BlochDict(save_Mz=$(s.save_Mz))") end -function sim_output_dim(obj::Phantom{T}, seq::Sequence, sys::Scanner, sim_method::BlochDict) where {T<:Real} +function sim_output_dim( + obj::Phantom{T}, seq::Sequence, sys::Scanner, sim_method::BlochDict +) where {T<:Real} out_state_dim = sim_method.save_Mz ? 2 : 1 return (sum(seq.ADC.N), length(obj), out_state_dim) end @@ -27,15 +29,18 @@ precession. - `S`: (`Vector{ComplexF64}`) raw signal over time - `M0`: (`::Vector{Mag}`) final state of the Mag vector """ -function run_spin_precession!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}}, - M::Mag{T}, sim_method::BlochDict) where {T<:Real} +function run_spin_precession!( + p::Phantom{T}, + seq::DiscreteSequence{T}, + sig::AbstractArray{Complex{T}}, + M::Mag{T}, + sim_method::BlochDict, +) where {T<:Real} #Simulation #Motion - xt = p.x .+ p.ux(p.x, p.y, p.z, seq.t') - yt = p.y .+ p.uy(p.x, p.y, p.z, seq.t') - zt = p.z .+ p.uz(p.x, p.y, p.z, seq.t') + x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') #Effective field - Bz = xt .* seq.Gx' .+ yt .* seq.Gy' .+ zt .* seq.Gz' .+ p.Δw / T(2π * γ) + Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw / T(2π * γ) #Rotation if is_ADC_on(seq) ϕ = T(-2π * γ) .* cumtrapz(seq.Δt', Bz) @@ -47,13 +52,12 @@ function run_spin_precession!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::Abst dur = sum(seq.Δt) # Total length, used for signal relaxation Mxy = [M.xy M.xy .* exp.(1im .* ϕ .- tp' ./ p.T2)] #This assumes Δw and T2 are constant in time M.xy .= Mxy[:, end] - #Acquired signal - sig[:,:,1] .= transpose(Mxy[:, findall(seq.ADC)]) + sig[:, :, 1] .= transpose(Mxy[:, findall(seq.ADC)]) if sim_method.save_Mz Mz = [M.z M.z .* exp.(-tp' ./ p.T1) .+ p.ρ .* (1 .- exp.(-tp' ./ p.T1))] #Calculate intermediate points - sig[:,:,2] .= transpose(Mz[:, findall(seq.ADC)]) #Save state to signal + sig[:, :, 2] .= transpose(Mz[:, findall(seq.ADC)]) #Save state to signal M.z .= Mz[:, end] else M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1)) #Jump to the last point diff --git a/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl b/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl index 0b088d613..65a5fed0d 100644 --- a/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl +++ b/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl @@ -5,16 +5,21 @@ export Bloch include("Magnetization.jl") #Defines Mag <: SpinStateRepresentation @functor Mag #Gives gpu acceleration capabilities, see GPUFunctions.jl -function sim_output_dim(obj::Phantom{T}, seq::Sequence, sys::Scanner, sim_method::SimulationMethod) where {T<:Real} +function sim_output_dim( + obj::Phantom{T}, seq::Sequence, sys::Scanner, sim_method::SimulationMethod +) where {T<:Real} return (sum(seq.ADC.N), 1) #Nt x Ncoils, This should consider the coil info from sys end """Magnetization initialization for Bloch simulation method.""" -function initialize_spins_state(obj::Phantom{T}, sim_method::SimulationMethod) where {T<:Real} +function initialize_spins_state( + obj::Phantom{T}, sim_method::SimulationMethod +) where {T<:Real} Nspins = length(obj) Mxy = zeros(T, Nspins) Mz = obj.ρ Xt = Mag{T}(Mxy, Mz) + sort_motions!(obj.motion) return Xt, obj end @@ -33,15 +38,18 @@ precession. - `S`: (`Vector{ComplexF64}`) raw signal over time - `M0`: (`::Vector{Mag}`) final state of the Mag vector """ -function run_spin_precession!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}}, - M::Mag{T}, sim_method::SimulationMethod) where {T<:Real} +function run_spin_precession!( + p::Phantom{T}, + seq::DiscreteSequence{T}, + sig::AbstractArray{Complex{T}}, + M::Mag{T}, + sim_method::SimulationMethod, +) where {T<:Real} #Simulation #Motion - xt = p.x .+ p.ux(p.x, p.y, p.z, seq.t') - yt = p.y .+ p.uy(p.x, p.y, p.z, seq.t') - zt = p.z .+ p.uz(p.x, p.y, p.z, seq.t') + x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') #Effective field - Bz = xt .* seq.Gx' .+ yt .* seq.Gy' .+ zt .* seq.Gz' .+ p.Δw / T(2π * γ) + Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw / T(2π * γ) #Rotation if is_ADC_on(seq) ϕ = T(-2π * γ) .* cumtrapz(seq.Δt', Bz) @@ -49,13 +57,14 @@ function run_spin_precession!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::Abst ϕ = T(-2π * γ) .* trapz(seq.Δt', Bz) end #Mxy precession and relaxation, and Mz relaxation - tp = cumsum(seq.Δt) # t' = t - t0 - dur = sum(seq.Δt) # Total length, used for signal relaxation - Mxy = [M.xy M.xy .* exp.(1im .* ϕ .- tp' ./ p.T2)] #This assumes Δw and T2 are constant in time + tp = cumsum(seq.Δt) # t' = t - t0 + dur = sum(seq.Δt) # Total length, used for signal relaxation + Mxy = [M.xy M.xy .* exp.(1im .* ϕ .- tp' ./ p.T2)] #This assumes Δw and T2 are constant in time M.xy .= Mxy[:, end] M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1)) #Acquired signal sig .= transpose(sum(Mxy[:, findall(seq.ADC)]; dims=1)) #<--- TODO: add coil sensitivities + return nothing end @@ -74,25 +83,28 @@ It gives rise to a rotation of `M0` with an angle given by the efective magnetic a part of the complete Mag vector and it's a part of the initial state for the next precession simulation step) """ -function run_spin_excitation!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}}, - M::Mag{T}, sim_method::SimulationMethod) where {T<:Real} +function run_spin_excitation!( + p::Phantom{T}, + seq::DiscreteSequence{T}, + sig::AbstractArray{Complex{T}}, + M::Mag{T}, + sim_method::SimulationMethod, +) where {T<:Real} #Simulation - for s ∈ seq #This iterates over seq, "s = seq[i,:]" + for s in seq #This iterates over seq, "s = seq[i,:]" #Motion - xt = p.x .+ p.ux(p.x, p.y, p.z, s.t) - yt = p.y .+ p.uy(p.x, p.y, p.z, s.t) - zt = p.z .+ p.uz(p.x, p.y, p.z, s.t) + x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, s.t) #Effective field - ΔBz = p.Δw ./ T(2π * γ) .- s.Δf ./ T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(xt,yt,zt) - Bz = (s.Gx .* xt .+ s.Gy .* yt .+ s.Gz .* zt) .+ ΔBz + ΔBz = p.Δw ./ T(2π * γ) .- s.Δf ./ T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z) + Bz = (s.Gx .* x .+ s.Gy .* y .+ s.Gz .* z) .+ ΔBz B = sqrt.(abs.(s.B1) .^ 2 .+ abs.(Bz) .^ 2) B[B .== 0] .= eps(T) #Spinor Rotation φ = T(-2π * γ) * (B .* s.Δt) # TODO: Use trapezoidal integration here (?), this is just Forward Euler - mul!( Q(φ, s.B1 ./ B, Bz ./ B), M ) + mul!(Q(φ, s.B1 ./ B, Bz ./ B), M) #Relaxation M.xy .= M.xy .* exp.(-s.Δt ./ p.T2) - M.z .= M.z .* exp.(-s.Δt ./ p.T1) .+ p.ρ .* (1 .- exp.(-s.Δt ./ p.T1)) + M.z .= M.z .* exp.(-s.Δt ./ p.T1) .+ p.ρ .* (1 .- exp.(-s.Δt ./ p.T1)) end #Acquired signal #sig .= -1.4im #<-- This was to test if an ADC point was inside an RF block diff --git a/KomaMRICore/src/simulation/GPUFunctions.jl b/KomaMRICore/src/simulation/GPUFunctions.jl index 79bb978bb..dff9725b4 100644 --- a/KomaMRICore/src/simulation/GPUFunctions.jl +++ b/KomaMRICore/src/simulation/GPUFunctions.jl @@ -11,36 +11,58 @@ Simple function to print the CUDA devices available in the host. function print_gpus() check_use_cuda() if use_cuda[] - cuda_devices = [Symbol("($(i-1)$(i == 1 ? "*" : " "))") => name(d) for (i,d) = enumerate(devices())] - @info "$(length(devices())) CUDA capable device(s)." cuda_devices... + cuda_devices = [ + Symbol("($(i-1)$(i == 1 ? "*" : " "))") => name(d) for + (i, d) in enumerate(devices()) + ] + @info "$(length(devices())) CUDA capable device(s)." cuda_devices... else @info "0 CUDA capable devices(s)." end - return + return nothing end """ Checks if the PC has a functional CUDA installation. Inspired by Flux's `check_use_cuda` funciton. """ function check_use_cuda() - if use_cuda[] === nothing - use_cuda[] = CUDA.functional() - if !(use_cuda[]) - @info """The GPU function is being called but the GPU is not accessible. - Defaulting back to the CPU. (No action is required if you want to run on the CPU).""" maxlog=1 - end - end + if use_cuda[] === nothing + use_cuda[] = CUDA.functional() + if !(use_cuda[]) + @info """The GPU function is being called but the GPU is not accessible. + Defaulting back to the CPU. (No action is required if you want to run on the CPU).""" maxlog = + 1 + end + end end #Aux. funcitons to check if the variable we want to convert to CuArray is numeric _isbitsarray(::AbstractArray{<:Real}) = true -_isbitsarray(::AbstractArray{T}) where T = isbitstype(T) +_isbitsarray(::AbstractArray{T}) where {T} = isbitstype(T) _isbitsarray(x) = false _isleaf(x) = _isbitsarray(x) || isleaf(x) # GPU adaptor struct KomaCUDAAdaptor end adapt_storage(to::KomaCUDAAdaptor, x) = CUDA.cu(x) +adapt_storage(to::KomaCUDAAdaptor, x::NoMotion) = NoMotion{Float32}() +adapt_storage(to::KomaCUDAAdaptor, x::SimpleMotion) = f32(x) +function adapt_storage(to::KomaCUDAAdaptor, x::ArbitraryMotion) + fields = [] + for field in fieldnames(ArbitraryMotion) + if field in (:ux, :uy, :uz) + push!(fields, adapt(KomaCUDAAdaptor(), getfield(x, field))) + else + push!(fields, f32(getfield(x, field))) + end + end + return ArbitraryMotion(fields...) +end +function adapt_storage( + to::KomaCUDAAdaptor, x::Vector{LinearInterpolator{T,V}} +) where {T<:Real,V<:AbstractVector{T}} + return CUDA.cu.(x) +end """ gpu(x) @@ -59,8 +81,8 @@ x = x |> gpu ``` """ function gpu(x) - check_use_cuda() - use_cuda[] ? fmap(x -> adapt(KomaCUDAAdaptor(), x), x; exclude = _isleaf) : x + check_use_cuda() + return use_cuda[] ? fmap(x -> adapt(KomaCUDAAdaptor(), x), x; exclude=_isleaf) : x end #CPU adaptor @@ -87,10 +109,19 @@ cpu(x) = fmap(x -> adapt(KomaCPUAdaptor(), x), x) #Precision paramtype(T::Type{<:Real}, m) = fmap(x -> adapt(T, x), m) - -adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Real}) = convert.(T, xs) #Type piracy -adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Complex}) = convert.(Complex{T}, xs) #Type piracy -adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Bool}) = xs #Type piracy +adapt_storage(T::Type{<:Real}, xs::Real) = convert(T, xs) +adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Real}) = convert.(T, xs) +adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Complex}) = convert.(Complex{T}, xs) +adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Bool}) = xs +adapt_storage(T::Type{<:Real}, xs::SimpleMotion) = SimpleMotion(paramtype(T, xs.types)) +adapt_storage(T::Type{<:Real}, xs::NoMotion) = NoMotion{T}() +function adapt_storage(T::Type{<:Real}, xs::ArbitraryMotion) + fields = [] + for field in fieldnames(ArbitraryMotion) + push!(fields, paramtype(T, getfield(xs, field))) + end + return ArbitraryMotion(fields...) +end """ f32(m) @@ -114,6 +145,14 @@ f64(m) = paramtype(Float64, m) #The functor macro makes it easier to call a function in all the parameters @functor Phantom + +@functor Translation +@functor Rotation +@functor HeartBeat +@functor PeriodicTranslation +@functor PeriodicRotation +@functor PeriodicHeartBeat + @functor Spinor @functor DiscreteSequence diff --git a/KomaMRICore/src/simulation/SimulatorCore.jl b/KomaMRICore/src/simulation/SimulatorCore.jl index e5b906121..31c6b0029 100644 --- a/KomaMRICore/src/simulation/SimulatorCore.jl +++ b/KomaMRICore/src/simulation/SimulatorCore.jl @@ -2,8 +2,8 @@ abstract type SimulationMethod end #get all available types by using subtypes(Ko abstract type SpinStateRepresentation{T<:Real} end #get all available types by using subtypes(KomaMRI.SpinStateRepresentation) #Defined methods: -include("Bloch/BlochSimulationMethod.jl") #Defines Bloch simulation method -include("Bloch/BlochDictSimulationMethod.jl") #Defines BlochDict simulation method +include("Bloch/BlochSimulationMethod.jl") #Defines Bloch simulation method +include("Bloch/BlochDictSimulationMethod.jl") #Defines BlochDict simulation method """ sim_params = default_sim_params(sim_params=Dict{String,Any}()) @@ -43,7 +43,11 @@ allowing the user to define some of them. """ function default_sim_params(sim_params=Dict{String,Any}()) sampling_params = KomaMRIBase.default_sampling_params() - get!(sim_params, "gpu", true); if sim_params["gpu"] check_use_cuda(); sim_params["gpu"] &= use_cuda[] end + get!(sim_params, "gpu", true) + if sim_params["gpu"] + check_use_cuda() + sim_params["gpu"] &= use_cuda[] + end get!(sim_params, "gpu_device", 0) get!(sim_params, "Nthreads", sim_params["gpu"] ? 1 : Threads.nthreads()) get!(sim_params, "Nblocks", 20) @@ -75,15 +79,21 @@ separating the spins of the phantom `obj` in `Nthreads`. next simulation step (the next step can be another precession step or an excitation step)) """ -function run_spin_precession_parallel!(obj::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}}, - Xt::SpinStateRepresentation{T}, sim_method::SimulationMethod; - Nthreads=Threads.nthreads()) where {T<:Real} - +function run_spin_precession_parallel!( + obj::Phantom{T}, + seq::DiscreteSequence{T}, + sig::AbstractArray{Complex{T}}, + Xt::SpinStateRepresentation{T}, + sim_method::SimulationMethod; + Nthreads=Threads.nthreads(), +) where {T<:Real} parts = kfoldperm(length(obj), Nthreads) - dims = [Colon() for i=1:ndims(sig)-1] # :,:,:,... Ndim times + dims = [Colon() for i in 1:(ndims(sig) - 1)] # :,:,:,... Ndim times ThreadsX.foreach(enumerate(parts)) do (i, p) - run_spin_precession!(@view(obj[p]), seq, @view(sig[dims...,i]), @view(Xt[p]), sim_method) + run_spin_precession!( + @view(obj[p]), seq, @view(sig[dims..., i]), @view(Xt[p]), sim_method + ) end return nothing @@ -108,15 +118,21 @@ different number threads to excecute the process. - `M0`: (`::Vector{Mag}`) final state of the Mag vector after a rotation (or the initial state for the next precession simulation step) """ -function run_spin_excitation_parallel!(obj::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}}, - Xt::SpinStateRepresentation{T}, sim_method::SimulationMethod; - Nthreads=Threads.nthreads()) where {T<:Real} - +function run_spin_excitation_parallel!( + obj::Phantom{T}, + seq::DiscreteSequence{T}, + sig::AbstractArray{Complex{T}}, + Xt::SpinStateRepresentation{T}, + sim_method::SimulationMethod; + Nthreads=Threads.nthreads(), +) where {T<:Real} parts = kfoldperm(length(obj), Nthreads) - dims = [Colon() for i=1:ndims(sig)-1] # :,:,:,... Ndim times + dims = [Colon() for i in 1:(ndims(sig) - 1)] # :,:,:,... Ndim times ThreadsX.foreach(enumerate(parts)) do (i, p) - run_spin_excitation!(@view(obj[p]), seq, @view(sig[dims...,i]), @view(Xt[p]), sim_method) + run_spin_excitation!( + @view(obj[p]), seq, @view(sig[dims..., i]), @view(Xt[p]), sim_method + ) end return nothing @@ -147,30 +163,49 @@ take advantage of CPU parallel processing. - `S_interp`: (`::Vector{ComplexF64}`) interpolated raw signal - `M0`: (`::Vector{Mag}`) final state of the Mag vector """ -function run_sim_time_iter!(obj::Phantom, seq::DiscreteSequence, sig::AbstractArray{Complex{T}}, - Xt::SpinStateRepresentation{T}, sim_method::SimulationMethod; - Nblocks=1, Nthreads=Threads.nthreads(), parts=[1:length(seq)], excitation_bool=ones(Bool, size(parts)), w=nothing) where {T<:Real} +function run_sim_time_iter!( + obj::Phantom, + seq::DiscreteSequence, + sig::AbstractArray{Complex{T}}, + Xt::SpinStateRepresentation{T}, + sim_method::SimulationMethod; + Nblocks=1, + Nthreads=Threads.nthreads(), + parts=[1:length(seq)], + excitation_bool=ones(Bool, size(parts)), + w=nothing, +) where {T<:Real} # Simulation rfs = 0 samples = 1 - progress_bar = Progress(Nblocks) - for (block, p) = enumerate(parts) + progress_bar = Progress(Nblocks; desc="Running simulation...") + + for (block, p) in enumerate(parts) seq_block = @view seq[p] # Params # excitation_bool = is_RF_on(seq_block) #&& is_ADC_off(seq_block) #PATCH: the ADC part should not be necessary, but sometimes 1 sample is identified as RF in an ADC block Nadc = sum(seq_block.ADC) - acq_samples = samples:samples+Nadc-1 - dims = [Colon() for i=1:ndims(sig)-1] # :,:,:,... Ndim times + acq_samples = samples:(samples + Nadc - 1) + dims = [Colon() for i in 1:(ndims(sig) - 1)] # :,:,:,... Ndim times # Simulation wrappers if excitation_bool[block] - run_spin_excitation_parallel!(obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method; Nthreads) + run_spin_excitation_parallel!( + obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method; Nthreads + ) rfs += 1 else - run_spin_precession_parallel!(obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method; Nthreads) + run_spin_precession_parallel!( + obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method; Nthreads + ) end samples += Nadc #Update progress - next!(progress_bar, showvalues=[(:simulated_blocks, block), (:rf_blocks, rfs), (:acq_samples, samples-1)]) + next!( + progress_bar; + showvalues=[ + (:simulated_blocks, block), (:rf_blocks, rfs), (:acq_samples, samples - 1) + ], + ) update_blink_window_progress!(w, block, Nblocks) end return nothing @@ -187,72 +222,71 @@ RF-on or RF-off. The function returns the ranges of the discrete sequence blocks with a boolean vector indicating whether each block has RF. """ function get_sim_ranges(seqd::DiscreteSequence; Nblocks) - ranges = UnitRange{Int}[] - ranges_bool = Bool[] - start_idx_rf_block = 0 - start_idx_gr_block = 0 - #Split 1:N into Nblocks like kfoldperm - N = length(seqd.Δt) - k = min(N, Nblocks) - n, r = divrem(N, k) #N >= k, N < k - breaks = collect(1:n:N+1) - for i in eachindex(breaks) - breaks[i] += i > r ? r : i-1 - end - breaks = breaks[2:end-1] #Remove borders, - #Iterate over B1 values to decide the simulation UnitRanges - for i in eachindex(seqd.Δt) - if abs(seqd.B1[i]) > 1e-9 #TODO: This is needed as the function ⏢ in get_rfs is not very accurate - if start_idx_rf_block == 0 #End RF block - start_idx_rf_block = i - end - if start_idx_gr_block > 0 #End of GR block - push!(ranges, start_idx_gr_block:i-1) - push!(ranges_bool, false) - start_idx_gr_block = 0 - end - else - if start_idx_gr_block == 0 #Start GR block - start_idx_gr_block = i - end - if start_idx_rf_block > 0 #End of RF block - push!(ranges, start_idx_rf_block:i-1) - push!(ranges_bool, true) - start_idx_rf_block = 0 - end - end - #More subdivisions - if i in breaks - if start_idx_rf_block > 0 #End of RF block - if length(start_idx_rf_block:i-1) > 1 - push!(ranges, start_idx_rf_block:i-1) - push!(ranges_bool, true) - start_idx_rf_block = i - end - end - if start_idx_gr_block > 0 #End of RF block - if length(start_idx_gr_block:i-1) > 1 - push!(ranges, start_idx_gr_block:i-1) - push!(ranges_bool, false) - start_idx_gr_block = i - end - end - end - end - #Finishing the UnitRange's - if start_idx_rf_block > 0 - push!(ranges, start_idx_rf_block:N) - push!(ranges_bool, true) - end - if start_idx_gr_block > 0 - push!(ranges, start_idx_gr_block:N) - push!(ranges_bool, false) - end - #Output - return ranges, ranges_bool + ranges = UnitRange{Int}[] + ranges_bool = Bool[] + start_idx_rf_block = 0 + start_idx_gr_block = 0 + #Split 1:N into Nblocks like kfoldperm + N = length(seqd.Δt) + k = min(N, Nblocks) + n, r = divrem(N, k) #N >= k, N < k + breaks = collect(1:n:(N + 1)) + for i in eachindex(breaks) + breaks[i] += i > r ? r : i - 1 + end + breaks = breaks[2:(end - 1)] #Remove borders, + #Iterate over B1 values to decide the simulation UnitRanges + for i in eachindex(seqd.Δt) + if abs(seqd.B1[i]) > 1e-9 #TODO: This is needed as the function ⏢ in get_rfs is not very accurate + if start_idx_rf_block == 0 #End RF block + start_idx_rf_block = i + end + if start_idx_gr_block > 0 #End of GR block + push!(ranges, start_idx_gr_block:(i - 1)) + push!(ranges_bool, false) + start_idx_gr_block = 0 + end + else + if start_idx_gr_block == 0 #Start GR block + start_idx_gr_block = i + end + if start_idx_rf_block > 0 #End of RF block + push!(ranges, start_idx_rf_block:(i - 1)) + push!(ranges_bool, true) + start_idx_rf_block = 0 + end + end + #More subdivisions + if i in breaks + if start_idx_rf_block > 0 #End of RF block + if length(start_idx_rf_block:(i - 1)) > 1 + push!(ranges, start_idx_rf_block:(i - 1)) + push!(ranges_bool, true) + start_idx_rf_block = i + end + end + if start_idx_gr_block > 0 #End of RF block + if length(start_idx_gr_block:(i - 1)) > 1 + push!(ranges, start_idx_gr_block:(i - 1)) + push!(ranges_bool, false) + start_idx_gr_block = i + end + end + end + end + #Finishing the UnitRange's + if start_idx_rf_block > 0 + push!(ranges, start_idx_rf_block:N) + push!(ranges_bool, true) + end + if start_idx_gr_block > 0 + push!(ranges, start_idx_gr_block:N) + push!(ranges_bool, false) + end + #Output + return ranges, ranges_bool end - """ out = simulate(obj::Phantom, seq::Sequence, sys::Scanner; sim_params, w) @@ -286,15 +320,15 @@ julia> plot_signal(raw) ``` """ function simulate( - obj::Phantom, seq::Sequence, sys::Scanner; - sim_params=Dict{String,Any}(), w=nothing + obj::Phantom, seq::Sequence, sys::Scanner; sim_params=Dict{String,Any}(), w=nothing ) #Simulation parameter unpacking, and setting defaults if key is not defined sim_params = default_sim_params(sim_params) # Simulation init seqd = discretize(seq; sampling_params=sim_params) # Sampling of Sequence waveforms parts, excitation_bool = get_sim_ranges(seqd; Nblocks=sim_params["Nblocks"]) # Generating simulation blocks - t_sim_parts = [seqd.t[p[1]] for p ∈ parts]; append!(t_sim_parts, seqd.t[end]) + t_sim_parts = [seqd.t[p[1]] for p in parts] + append!(t_sim_parts, seqd.t[end]) # Spins' state init (Magnetization, EPG, etc.), could include modifications to obj (e.g. T2*) Xt, obj = initialize_spins_state(obj, sim_params["sim_method"]) # Signal init @@ -303,31 +337,49 @@ function simulate( # Objects to GPU if sim_params["gpu"] #Default device!(sim_params["gpu_device"]) - gpu_name = name.(devices())[sim_params["gpu_device"]+1] - obj = obj |> gpu #Phantom + gpu_name = name.(devices())[sim_params["gpu_device"] + 1] + obj = obj |> gpu #Phantom seqd = seqd |> gpu #DiscreteSequence - Xt = Xt |> gpu #SpinStateRepresentation - sig = sig |> gpu #Signal + Xt = Xt |> gpu #SpinStateRepresentation + sig = sig |> gpu #Signal end if sim_params["precision"] == "f32" #Default - obj = obj |> f32 #Phantom + obj = obj |> f32 #Phantom seqd = seqd |> f32 #DiscreteSequence - Xt = Xt |> f32 #SpinStateRepresentation - sig = sig |> f32 #Signal + Xt = Xt |> f32 #SpinStateRepresentation + sig = sig |> f32 #Signal elseif sim_params["precision"] == "f64" - obj = obj |> f64 #Phantom + obj = obj |> f64 #Phantom seqd = seqd |> f64 #DiscreteSequence - Xt = Xt |> f64 #SpinStateRepresentation - sig = sig |> f64 #Signal + Xt = Xt |> f64 #SpinStateRepresentation + sig = sig |> f64 #Signal end + # Simulation - @info "Running simulation in the $(sim_params["gpu"] ? "GPU ($gpu_name)" : "CPU with $(sim_params["Nthreads"]) thread(s)")" koma_version=__VERSION__ sim_method = sim_params["sim_method"] spins = length(obj) time_points = length(seqd.t) adc_points=Ndims[1] - @time timed_tuple = @timed run_sim_time_iter!(obj, seqd, sig, Xt, sim_params["sim_method"]; Nblocks=length(parts), Nthreads=sim_params["Nthreads"], parts, excitation_bool, w) + @info "Running simulation in the $(sim_params["gpu"] ? "GPU ($gpu_name)" : "CPU with $(sim_params["Nthreads"]) thread(s)")" koma_version = + __VERSION__ sim_method = sim_params["sim_method"] spins = length(obj) time_points = length( + seqd.t + ) adc_points = Ndims[1] + @time timed_tuple = @timed run_sim_time_iter!( + obj, + seqd, + sig, + Xt, + sim_params["sim_method"]; + Nblocks=length(parts), + Nthreads=sim_params["Nthreads"], + parts, + excitation_bool, + w, + ) # Result to CPU, if already in the CPU it does nothing - sig = sum(sig; dims=length(Ndims)+1) |> cpu #Sum over threads + sig = sum(sig; dims=length(Ndims) + 1) |> cpu #Sum over threads sig .*= get_adc_phase_compensation(seq) Xt = Xt |> cpu - if sim_params["gpu"] GC.gc(true); CUDA.reclaim() end + if sim_params["gpu"] + GC.gc(true) + CUDA.reclaim() + end # Output if sim_params["return_type"] == "state" out = Xt @@ -344,7 +396,10 @@ function simulate( sim_params_raw["Nblocks"] = length(parts) sim_params_raw["sim_time_sec"] = timed_tuple.time sim_params_raw["allocations_bytes"] = timed_tuple.bytes - out = signal_to_raw_data(sig, seq; phantom_name=obj.name, sys=sys, sim_params=sim_params_raw) + + out = signal_to_raw_data( + sig, seq; phantom_name=obj.name, sys=sys, sim_params=sim_params_raw + ) end return out end @@ -366,12 +421,11 @@ Returns magnetization of spins distributed along `z` after running the Sequence - `mag`: (`::SpinStateRepresentation`) final state of the magnetization vector """ function simulate_slice_profile( - seq::Sequence; - z=range(-2.e-2, 2.e-2, 200), sim_params=Dict{String,Any}("Δt_rf" => 1e-6) + seq::Sequence; z=range(-2.e-2, 2.e-2, 200), sim_params=Dict{String,Any}("Δt_rf" => 1e-6) ) sim_params["return_type"] = "state" sys = Scanner() - obj = Phantom{Float64}(x=zeros(size(z)), z=Array(z)) + obj = Phantom{Float64}(; x=zeros(size(z)), z=Array(z)) mag = simulate(obj, seq, sys; sim_params) return mag end diff --git a/KomaMRICore/test/runtests.jl b/KomaMRICore/test/runtests.jl index cbd456429..9306533e0 100644 --- a/KomaMRICore/test/runtests.jl +++ b/KomaMRICore/test/runtests.jl @@ -166,7 +166,7 @@ end using Suppressor include(joinpath(@__DIR__, "test_files", "utils.jl")) - sig_jemris = signal_jemris() + sig_jemris = signal_sphere_jemris() seq = seq_epi_100x100_TE100_FOV230() obj = phantom_sphere() sys = Scanner() @@ -189,7 +189,7 @@ end using Suppressor include(joinpath(@__DIR__, "test_files", "utils.jl")) - sig_jemris = signal_jemris() + sig_jemris = signal_sphere_jemris() seq = seq_epi_100x100_TE100_FOV230() obj = phantom_sphere() sys = Scanner() @@ -208,11 +208,11 @@ end end -@testitem "Bloch_GPU" tags=[:important, :skipci, :core] begin +@testitem "Bloch_GPU" tags=[:important, :skipci, :core, :gpu] begin using Suppressor include(joinpath(@__DIR__, "test_files", "utils.jl")) - sig_jemris = signal_jemris() + sig_jemris = signal_sphere_jemris() seq = seq_epi_100x100_TE100_FOV230() obj = phantom_sphere() sys = Scanner() @@ -220,7 +220,8 @@ end sim_params = Dict{String, Any}( "gpu"=>true, "sim_method"=>KomaMRICore.Bloch(), - "return_type"=>"mat" + "return_type"=>"mat", + "precision"=>"f64" ) sig = @suppress simulate(obj, seq, sys; sim_params) sig = sig / prod(size(obj)) @@ -242,7 +243,7 @@ end N = 6 sys = Scanner() - obj = Phantom{Float64}(x=[0],T1=[T1],T2=[T2],Δw=[Δw]) + obj = Phantom{Float64}(x=[0.],T1=[T1],T2=[T2],Δw=[Δw]) rf_phase = [0, π/2] seq = Sequence() @@ -288,7 +289,7 @@ end N = 6 sys = Scanner() - obj = Phantom{Float64}(x=[0],T1=[T1],T2=[T2],Δw=[Δw]) + obj = Phantom{Float64}(x=[0.],T1=[T1],T2=[T2],Δw=[Δw]) rf_phase = [0, π/2] seq = Sequence() @@ -323,7 +324,7 @@ end @test error0 + error1 + error2 < 0.1 #NMRSE < 0.1% end -@testitem "Bloch_GPU_RF_accuracy" tags=[:important, :core] begin +@testitem "Bloch_GPU_RF_accuracy" tags=[:important, :core, :skipci, :gpu] begin using Suppressor Tadc = 1e-3 @@ -335,7 +336,7 @@ end N = 6 sys = Scanner() - obj = Phantom{Float64}(x=[0],T1=[T1],T2=[T2],Δw=[Δw]) + obj = Phantom{Float64}(x=[0.],T1=[T1],T2=[T2],Δw=[Δw]) rf_phase = [0, π/2] seq = Sequence() @@ -381,7 +382,7 @@ end N = 6 sys = Scanner() - obj = Phantom{Float64}(x=[0],T1=[T1],T2=[T2],Δw=[Δw]) + obj = Phantom{Float64}(x=[0.],T1=[T1],T2=[T2],Δw=[Δw]) rf_phase = 2π*rand() seq1 = Sequence() @@ -400,6 +401,156 @@ end end +@testitem "Bloch CPU_single_thread SimpleMotion" tags=[:important, :core] begin + using Suppressor + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain() + obj.motion = SimpleMotion([Translation(t_end=10.0, dx=0.0, dy=1.0, dz=0.0)]) + sim_params = Dict{String, Any}( + "gpu"=>false, + "Nthreads"=>1, + "sim_method"=>KomaMRICore.Bloch(), + "return_type"=>"mat" + ) + sig = @suppress simulate(obj, seq, sys; sim_params) + sig = sig / prod(size(obj)) + NMRSE(x, x_true) = sqrt.( sum(abs.(x .- x_true).^2) ./ sum(abs.(x_true).^2) ) * 100. + @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% +end + +@testitem "Bloch CPU_single_thread ArbitraryMotion" tags=[:important, :core] begin + using Suppressor + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain() + Ns = length(obj) + period_durations=[20.0] + dx = dz = zeros(Ns, 1) + dy = 1.0 .* ones(Ns, 1) + obj.motion = @suppress ArbitraryMotion( + period_durations, + dx, + dy, + dz) + sim_params = Dict{String, Any}( + "gpu"=>false, + "Nthreads"=>1, + "sim_method"=>KomaMRICore.Bloch(), + "return_type"=>"mat" + ) + sig = @suppress simulate(obj, seq, sys; sim_params) + sig = sig / prod(size(obj)) + NMRSE(x, x_true) = sqrt.( sum(abs.(x .- x_true).^2) ./ sum(abs.(x_true).^2) ) * 100. + @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% +end + + +@testitem "Bloch CPU_multi_thread SimpleMotion" tags=[:important, :core] begin + using Suppressor + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain() + obj.motion = SimpleMotion([Translation(t_end=10.0, dx=0.0, dy=1.0, dz=0.0)]) + sim_params = Dict{String, Any}( + "gpu"=>false, + "sim_method"=>KomaMRICore.Bloch(), + "return_type"=>"mat" + ) + sig = @suppress simulate(obj, seq, sys; sim_params) + sig = sig / prod(size(obj)) + NMRSE(x, x_true) = sqrt.( sum(abs.(x .- x_true).^2) ./ sum(abs.(x_true).^2) ) * 100. + @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% +end + +@testitem "Bloch CPU_multi_thread ArbitraryMotion" tags=[:important, :core] begin + using Suppressor + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain() + Ns = length(obj) + period_durations=[20.0] + dx = dz = zeros(Ns, 1) + dy = 1.0 .* ones(Ns, 1) + obj.motion = @suppress ArbitraryMotion( + period_durations, + dx, + dy, + dz) + sim_params = Dict{String, Any}( + "gpu"=>false, + "sim_method"=>KomaMRICore.Bloch(), + "return_type"=>"mat" + ) + sig = @suppress simulate(obj, seq, sys; sim_params) + sig = sig / prod(size(obj)) + NMRSE(x, x_true) = sqrt.( sum(abs.(x .- x_true).^2) ./ sum(abs.(x_true).^2) ) * 100. + @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% +end + +@testitem "Bloch GPU SimpleMotion" tags=[:important, :core, :skipci, :gpu] begin + using Suppressor + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain() + obj.motion = SimpleMotion([Translation(t_end=10.0, dx=0.0, dy=1.0, dz=0.0)]) + sim_params = Dict{String, Any}( + "gpu"=>true, + "sim_method"=>KomaMRICore.Bloch(), + "return_type"=>"mat", + "precision"=>"f64" + ) + sig = @suppress simulate(obj, seq, sys; sim_params) + sig = sig / prod(size(obj)) + NMRSE(x, x_true) = sqrt.( sum(abs.(x .- x_true).^2) ./ sum(abs.(x_true).^2) ) * 100. + @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% +end + +@testitem "Bloch GPU ArbitraryMotion" tags=[:important, :core, :skipci, :gpu] begin + using Suppressor + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain() + Ns = length(obj) + period_durations=[20.0] + dx = dz = zeros(Ns, 1) + dy = 1.0 .* ones(Ns, 1) + obj.motion = @suppress ArbitraryMotion( + period_durations, + dx, + dy, + dz) + sim_params = Dict{String, Any}( + "gpu"=>true, + "sim_method"=>KomaMRICore.Bloch(), + "return_type"=>"mat", + "precision"=>"f64" + ) + sig = @suppress simulate(obj, seq, sys; sim_params) + sig = sig / prod(size(obj)) + NMRSE(x, x_true) = sqrt.( sum(abs.(x .- x_true).^2) ./ sum(abs.(x_true).^2) ) * 100. + @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% +end + + @testitem "BlochDict_CPU_single_thread" tags=[:important, :core] begin using Suppressor include(joinpath(@__DIR__, "test_files", "utils.jl")) diff --git a/KomaMRICore/test/test_files/jemris_signals_epi_brain_motion.h5 b/KomaMRICore/test/test_files/jemris_signals_epi_brain_motion.h5 new file mode 100644 index 000000000..89955503e Binary files /dev/null and b/KomaMRICore/test/test_files/jemris_signals_epi_brain_motion.h5 differ diff --git a/KomaMRICore/test/test_files/utils.jl b/KomaMRICore/test/test_files/utils.jl index 553bd8706..d10066810 100644 --- a/KomaMRICore/test/test_files/utils.jl +++ b/KomaMRICore/test/test_files/utils.jl @@ -1,6 +1,6 @@ using HDF5 -function signal_jemris() +function signal_sphere_jemris() path = @__DIR__ sig = h5open(joinpath(path, "jemris_signals_epi_sphere_cs.h5"))["/signal/channels/00"] sig = sig[1,:] + 1im*sig[2,:] @@ -8,6 +8,14 @@ function signal_jemris() return sig end +function signal_brain_motion_jemris() + path = @__DIR__ + sig = h5open(joinpath(path, "jemris_signals_epi_brain_motion.h5"))["/signal/channels/00"] + sig = sig[1,:] + 1im*sig[2,:] + sig = sig[:] + return sig +end + function phantom_sphere() path = @__DIR__ fid = h5open(joinpath(path, "phantom_sphere.h5"), "r") @@ -16,6 +24,47 @@ function phantom_sphere() return Phantom(; name="sphere", x, y, z, ρ, T1, T2, T2s, Δw) end +function phantom_brain() + path = @__DIR__ + fid = h5open(joinpath(path, "../../../examples/2.phantoms/brain.h5"), "r") + data = read(fid["sample/data"]) + Δx = read(fid["sample/resolution"]) * 1e-3 #[m] + offset = read(fid["sample/offset"]) * 1e-3 #[m] + mask = data[1, :, :, :] .!= 0 + #Maps + ρ = data[1, :, :, :] + T1 = 1e-3 ./ data[2, :, :, :] + T2 = 1e-3 ./ data[3, :, :, :] + T2s = 1e-3 ./ data[4, :, :, :] + Δw = data[5, :, :, :] + #Positions + X, Y, Z = size(ρ) + FOVx = (X - 1) * Δx[1] #[m] + FOVy = (Y - 1) * Δx[2] #[m] + FOVz = (Z - 1) * Δx[3] #[m] + xx = reshape(((-FOVx / 2):Δx[1]:(FOVx / 2)), :, 1, 1) #[(-FOVx/2:Δx[1]:FOVx/2)...;] + yy = reshape(((-FOVy / 2):Δx[2]:(FOVy / 2)), 1, :, 1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;] + zz = reshape(((-FOVz / 2):Δx[3]:(FOVz / 2)), 1, 1, :) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;] + x = xx * 1 .+ yy * 0 .+ zz * 0 .+ offset[1] #spin x coordinates + y = xx * 0 .+ yy * 1 .+ zz * 0 .+ offset[2] #spin y coordinates + z = xx * 0 .+ yy * 0 .+ zz * 1 .+ offset[3] #spin z coordinates + v = 0 # m/s + + obj = Phantom(; + name="brain", + x=x[mask], + y=y[mask], + z=z[mask], + ρ=ρ[mask], + T1=T1[mask], + T2=T2[mask], + T2s=T2s[mask], + Δw=Δw[mask], + ) + + return obj +end + function seq_epi_100x100_TE100_FOV230() # Gyromagnetic constant diff --git a/KomaMRIFiles/Project.toml b/KomaMRIFiles/Project.toml index 65d422118..8957ee4c2 100644 --- a/KomaMRIFiles/Project.toml +++ b/KomaMRIFiles/Project.toml @@ -6,6 +6,7 @@ version = "0.8.2" [deps] FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" KomaMRIBase = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c" MAT = "23992714-dd62-5051-b70f-ba57cb901cac" MRIFiles = "5a6f062f-bf45-497d-b654-ad17aae2a530" diff --git a/KomaMRIFiles/src/KomaMRIFiles.jl b/KomaMRIFiles/src/KomaMRIFiles.jl index b4304bb9e..259639a62 100644 --- a/KomaMRIFiles/src/KomaMRIFiles.jl +++ b/KomaMRIFiles/src/KomaMRIFiles.jl @@ -1,7 +1,7 @@ module KomaMRIFiles using KomaMRIBase -using Scanf, FileIO, HDF5, MAT # IO related +using Scanf, FileIO, HDF5, MAT, InteractiveUtils # IO related using Reexport using MRIFiles @@ -12,9 +12,10 @@ import MRIFiles: insertNode include("Sequence/Pulseq.jl") include("Phantom/JEMRIS.jl") include("Phantom/MRiLab.jl") +include("Phantom/Phantom.jl") -export read_seq # Pulseq -export read_phantom_jemris, read_phantom_MRiLab # Phantom +export read_seq # Pulseq +export read_phantom_jemris, read_phantom_MRiLab, read_phantom, write_phantom # Phantom # Package version: KomaMRIFiles.__VERSION__ using Pkg diff --git a/KomaMRIFiles/src/Phantom/JEMRIS.jl b/KomaMRIFiles/src/Phantom/JEMRIS.jl index 239160cbb..05039aa0d 100644 --- a/KomaMRIFiles/src/Phantom/JEMRIS.jl +++ b/KomaMRIFiles/src/Phantom/JEMRIS.jl @@ -19,49 +19,49 @@ julia> plot_phantom_map(obj, :ρ) ``` """ function read_phantom_jemris(filename) - # A(:,:,:,1)=Sample.M0; - # I=find(Sample.T1); R1 =zeros(size(Sample.T1)); R1(I) =1./Sample.T1(I); - # I=find(Sample.T2); R2 =zeros(size(Sample.T2)); R2(I) =1./Sample.T2(I); - # I=find(Sample.T2S);R2S=zeros(size(Sample.T2S)); R2S(I)=1./Sample.T2S(I); - # A(:,:,:,2)=R1; #1/T1 - # A(:,:,:,3)=R2; #1/T2 - # A(:,:,:,4)=R2S; #1/T2s - # A(:,:,:,5)=Sample.DB; - fid = HDF5.h5open(filename) - data = read(fid["sample/data"]) - Δx = read(fid["sample/resolution"]) * 1e-3 #[m] - offset = read(fid["sample/offset"]) * 1e-3 #[m] - mask = data[1,:,:,:] .!= 0 - #Maps - ρ = data[1,:,:,:] - T1 = 1e-3 ./ data[2,:,:,:] - T2 = 1e-3 ./ data[3,:,:,:] - T2s = 1e-3 ./ data[4,:,:,:] - Δw = data[5,:,:,:] - #Positions - X, Y, Z = size(ρ) - FOVx = (X-1)*Δx[1] #[m] - FOVy = (Y-1)*Δx[2] #[m] - FOVz = (Z-1)*Δx[3] #[m] - xx = reshape((-FOVx/2:Δx[1]:FOVx/2),:,1,1) #[(-FOVx/2:Δx[1]:FOVx/2)...;] - yy = reshape((-FOVy/2:Δx[2]:FOVy/2),1,:,1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;] - zz = reshape((-FOVz/2:Δx[3]:FOVz/2),1,1,:) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;] - x = xx*1 .+ yy*0 .+ zz*0 .+ offset[1] #spin x coordinates - y = xx*0 .+ yy*1 .+ zz*0 .+ offset[2] #spin y coordinates - z = xx*0 .+ yy*0 .+ zz*1 .+ offset[3] #spin z coordinates - v = 0 # m/s + # A(:,:,:,1)=Sample.M0; + # I=find(Sample.T1); R1 =zeros(size(Sample.T1)); R1(I) =1./Sample.T1(I); + # I=find(Sample.T2); R2 =zeros(size(Sample.T2)); R2(I) =1./Sample.T2(I); + # I=find(Sample.T2S);R2S=zeros(size(Sample.T2S)); R2S(I)=1./Sample.T2S(I); + # A(:,:,:,2)=R1; #1/T1 + # A(:,:,:,3)=R2; #1/T2 + # A(:,:,:,4)=R2S; #1/T2s + # A(:,:,:,5)=Sample.DB; + fid = HDF5.h5open(filename) + data = read(fid["sample/data"]) + Δx = read(fid["sample/resolution"]) * 1e-3 #[m] + offset = read(fid["sample/offset"]) * 1e-3 #[m] + mask = data[1, :, :, :] .!= 0 + #Maps + ρ = data[1, :, :, :] + T1 = 1e-3 ./ data[2, :, :, :] + T2 = 1e-3 ./ data[3, :, :, :] + T2s = 1e-3 ./ data[4, :, :, :] + Δw = data[5, :, :, :] + #Positions + X, Y, Z = size(ρ) + FOVx = (X - 1) * Δx[1] #[m] + FOVy = (Y - 1) * Δx[2] #[m] + FOVz = (Z - 1) * Δx[3] #[m] + xx = reshape(((-FOVx / 2):Δx[1]:(FOVx / 2)), :, 1, 1) #[(-FOVx/2:Δx[1]:FOVx/2)...;] + yy = reshape(((-FOVy / 2):Δx[2]:(FOVy / 2)), 1, :, 1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;] + zz = reshape(((-FOVz / 2):Δx[3]:(FOVz / 2)), 1, 1, :) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;] + x = xx * 1 .+ yy * 0 .+ zz * 0 .+ offset[1] #spin x coordinates + y = xx * 0 .+ yy * 1 .+ zz * 0 .+ offset[2] #spin y coordinates + z = xx * 0 .+ yy * 0 .+ zz * 1 .+ offset[3] #spin z coordinates + v = 0 # m/s - obj = Phantom{Float64}( - name = basename(filename), - x = x[mask], - y = y[mask], - z = z[mask], - ρ = ρ[mask], - T1 = T1[mask], - T2 = T2[mask], - T2s = T2s[mask], - Δw = Δw[mask], - ux = (x,y,z,t)->v*t, - ) - return obj + obj = Phantom(; + name=basename(filename), + x=x[mask], + y=y[mask], + z=z[mask], + ρ=ρ[mask], + T1=T1[mask], + T2=T2[mask], + T2s=T2s[mask], + Δw=Δw[mask], + ) + + return obj end diff --git a/KomaMRIFiles/src/Phantom/MRiLab.jl b/KomaMRIFiles/src/Phantom/MRiLab.jl index 26a6696a4..f6f5ee3d4 100644 --- a/KomaMRIFiles/src/Phantom/MRiLab.jl +++ b/KomaMRIFiles/src/Phantom/MRiLab.jl @@ -18,55 +18,55 @@ julia> obj = read_phantom_MRiLab(obj_file) julia> plot_phantom_map(obj, :ρ) ``` """ -function read_phantom_MRiLab(filename; B0=1.5, offset=[0,0,0], FRange_filename="") +function read_phantom_MRiLab(filename; B0=1.5, offset=[0, 0, 0], FRange_filename="") data = MAT.matread(filename)["VObj"] - Δx = [data["XDimRes"], data["YDimRes"], data["ZDimRes"]] #[m] - mask = data["Rho"] .!= 0 - if FRange_filename!="" - data_mask = MAT.matread(FRange_filename)["VMag"] - mask .*= data_mask["FRange"] - end - #Maps - ρ = data["Rho"] - T1 = data["T1"] - T2 = data["T2"] - T2s = data["T2Star"] + Δx = [data["XDimRes"], data["YDimRes"], data["ZDimRes"]] #[m] + mask = data["Rho"] .!= 0 + if FRange_filename != "" + data_mask = MAT.matread(FRange_filename)["VMag"] + mask .*= data_mask["FRange"] + end + #Maps + ρ = data["Rho"] + T1 = data["T1"] + T2 = data["T2"] + T2s = data["T2Star"] CS = data["ChemShift"] #[Hz/T] - Δw = 2π .* CS .* B0 .* ones(size(ρ)) - #Positions - X, Y, Z = size(ρ) - FOVx = (X-1)*Δx[1] #[m] - FOVy = (Y-1)*Δx[2] #[m] - FOVz = (Z-1)*Δx[3] #[m] + Δw = 2π .* CS .* B0 .* ones(size(ρ)) + #Positions + X, Y, Z = size(ρ) + FOVx = (X - 1) * Δx[1] #[m] + FOVy = (Y - 1) * Δx[2] #[m] + FOVz = (Z - 1) * Δx[3] #[m] # [row,col,layer]=size(VOex.Mz); # VVar.ObjLoc = - # [((col+1)/2)*VOex.XDimRes; - # ((row+1)/2)*VOex.YDimRes ; - # ((layer+1)/2)*VOex.ZDimRes]; % Set matrix center as Object position for motion simulation + # [((col+1)/2)*VOex.XDimRes; + # ((row+1)/2)*VOex.YDimRes ; + # ((layer+1)/2)*VOex.ZDimRes]; % Set matrix center as Object position for motion simulation # VVar.ObjTurnLoc = - # [((col+1)/2)*VOex.XDimRes; - # ((row+1)/2)*VOex.YDimRes ; - # ((layer+1)/2)*VOex.ZDimRes]; % Set matrix center as Object origin for motion simulation - xx = reshape((-FOVx/2:Δx[1]:FOVx/2),:,1,1) #[(-FOVx/2:Δx[1]:FOVx/2)...;] - yy = reshape((-FOVy/2:Δx[2]:FOVy/2),1,:,1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;] - zz = reshape((-FOVz/2:Δx[3]:FOVz/2),1,1,:) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;] - x = xx*1 .+ yy*0 .+ zz*0 .+ offset[1] #spin x coordinates - y = xx*0 .+ yy*1 .+ zz*0 .+ offset[2] #spin y coordinates - z = xx*0 .+ yy*0 .+ zz*1 .+ offset[3] #spin z coordinates - v = 0 # m/s + # [((col+1)/2)*VOex.XDimRes; + # ((row+1)/2)*VOex.YDimRes ; + # ((layer+1)/2)*VOex.ZDimRes]; % Set matrix center as Object origin for motion simulation + xx = reshape(((-FOVx / 2):Δx[1]:(FOVx / 2)), :, 1, 1) #[(-FOVx/2:Δx[1]:FOVx/2)...;] + yy = reshape(((-FOVy / 2):Δx[2]:(FOVy / 2)), 1, :, 1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;] + zz = reshape(((-FOVz / 2):Δx[3]:(FOVz / 2)), 1, 1, :) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;] + x = xx * 1 .+ yy * 0 .+ zz * 0 .+ offset[1] #spin x coordinates + y = xx * 0 .+ yy * 1 .+ zz * 0 .+ offset[2] #spin y coordinates + z = xx * 0 .+ yy * 0 .+ zz * 1 .+ offset[3] #spin z coordinates + v = 0 # m/s - obj = Phantom( - name = basename(filename), - x = y[mask], - y = -x[mask], - z = -z[mask], - ρ = ρ[mask], - T1 = T1[mask], - T2 = T2[mask], - T2s = T2s[mask], - Δw = Δw[mask], - ux = (x,y,z,t)->v*t, - ) - return obj + obj = Phantom(; + name=basename(filename), + x=x[mask], + y=y[mask], + z=z[mask], + ρ=ρ[mask], + T1=T1[mask], + T2=T2[mask], + T2s=T2s[mask], + Δw=Δw[mask], + ) + + return obj end diff --git a/KomaMRIFiles/src/Phantom/Phantom.jl b/KomaMRIFiles/src/Phantom/Phantom.jl new file mode 100644 index 000000000..cd2bf15d8 --- /dev/null +++ b/KomaMRIFiles/src/Phantom/Phantom.jl @@ -0,0 +1,177 @@ +""" + phantom = read_phantom(filename) + +Reads a (.phantom) file and creates a Phantom structure from it +""" +function read_phantom(filename::String) + fid = HDF5.h5open(filename, "r") + fields = [] + version = read_attribute(fid, "Version") + dims = read_attribute(fid, "Dims") + Ns = read_attribute(fid, "Ns") + # Name + name = read_attribute(fid, "Name") + push!(fields, (:name, name)) + # Position and contrast + for key in ["position", "contrast"] + group = fid[key] + for label in HDF5.keys(group) + param = group[label] + values = read_param(param) + if values != "Default" + push!(fields, (Symbol(label), values)) + end + end + end + # Motion + motion_group = fid["motion"] + model = read_attribute(motion_group, "model") + import_motion!(fields, Ns, Symbol(model), motion_group) + + obj = Phantom(; fields...) + close(fid) + return obj +end + +function read_param(param::HDF5.Group) + if "type" in HDF5.keys(attrs(param)) + type = attrs(param)["type"] + if type == "Explicit" + values = read(param["values"]) + elseif type == "Indexed" + index = read(param["values"]) + @assert Ns == length(index) "Error: $(label) vector dimensions mismatch" + table = read(param["table"]) + N = read_attribute(param, "N") + @assert N == length(table) "Error: $(label) table dimensions mismatch" + values = table[index] + elseif type == "Default" + values = "Default" + end + else + values = read(param["values"]) + end + return values +end + +function import_motion!(fields::Array, Ns::Int, model::Symbol, motion_group::HDF5.Group) + return import_motion!(fields, Ns, Val(model), motion_group) +end +function import_motion!( + fields::Array, Ns::Int, model::Val{:NoMotion}, motion_group::HDF5.Group +) + return nothing +end +function import_motion!( + fields::Array, Ns::Int, model::Val{:SimpleMotion}, motion_group::HDF5.Group +) + types_group = motion_group["types"] + types = SimpleMotionType[] + for key in keys(types_group) + type_group = types_group[key] + type_str = split(key, "_")[2] + @assert type_str in last.(split.(string.(subtypes(SimpleMotionType)), ".")) "Simple Motion Type: $(type_str) has not been implemented in KomaMRIBase $(KomaMRIBase.__VERSION__)" + for SMT in subtypes(SimpleMotionType) + args = [] + if type_str == last(split(string(SMT), ".")) + for key in fieldnames(SMT) + push!(args, read_attribute(type_group, string(key))) + end + push!(types, SMT(args...)) + end + end + end + return push!(fields, (:motion, SimpleMotion(vcat(types...)))) +end +function import_motion!( + fields::Array, Ns::Int, model::Val{:ArbitraryMotion}, motion_group::HDF5.Group +) + dur = read(motion_group["period_durations"]) + K = read_attribute(motion_group, "K") + dx = zeros(Ns, K - 1) + dy = zeros(Ns, K - 1) + dz = zeros(Ns, K - 1) + for key in HDF5.keys(motion_group) + if key != "period_durations" + values = read_param(motion_group[key]) + if key == "dx" + dx = values + elseif key == "dy" + dy = values + elseif key == "dz" + dz = values + end + end + end + return push!(fields, (:motion, ArbitraryMotion(dur, dx, dy, dz))) +end + +""" + phantom = write_phantom(ph,filename) + +Writes a (.phantom) file from a Phantom struct. +""" +function write_phantom( + # By the moment, only "Explicit" type + # is considered when writing .phantom files + obj::Phantom, + filename::String; + store_coords=[:x, :y, :z], + store_contrasts=[:ρ, :T1, :T2, :T2s, :Δw], + store_motion=true, +) + # Create HDF5 phantom file + fid = h5open(filename, "w") + # Root attributes + HDF5.attributes(fid)["Version"] = string(KomaMRIFiles.__VERSION__) + HDF5.attributes(fid)["Name"] = obj.name + HDF5.attributes(fid)["Ns"] = length(obj.x) + dims = KomaMRIBase.get_dims(obj) + HDF5.attributes(fid)["Dims"] = sum(dims) + # Positions + pos = create_group(fid, "position") + for x in store_coords + create_group(pos, String(x))["values"] = getfield(obj, x) + end + # Contrast (Rho, T1, T2, T2s Deltaw) + contrast = create_group(fid, "contrast") + for x in store_contrasts + param = create_group(contrast, String(x)) + HDF5.attributes(param)["type"] = "Explicit" #TODO: consider "Indexed" type + param["values"] = getfield(obj, x) + end + # Motion + if store_motion + motion_group = create_group(fid, "motion") + export_motion!(motion_group, obj.motion) + end + return close(fid) +end + +function export_motion!(motion_group::HDF5.Group, motion::NoMotion) + return HDF5.attributes(motion_group)["model"] = "NoMotion" +end +function export_motion!(motion_group::HDF5.Group, motion::SimpleMotion) + HDF5.attributes(motion_group)["model"] = "SimpleMotion" + types_group = create_group(motion_group, "types") + counter = 1 + for (counter, sm_type) in enumerate(motion.types) + simple_motion_type = typeof(sm_type).name.name + type_group = create_group(types_group, "$(counter)_$simple_motion_type") + fields = fieldnames(typeof(sm_type)) + for field in fields + HDF5.attributes(type_group)[string(field)] = getfield(sm_type, field) + end + end +end +function export_motion!(motion_group::HDF5.Group, motion::ArbitraryMotion) + HDF5.attributes(motion_group)["model"] = "ArbitraryMotion" + HDF5.attributes(motion_group)["K"] = size(motion.dx)[2] + 1 + motion_group["period_durations"] = motion.period_durations + + for key in ["dx", "dy", "dz"] + d_group = create_group(motion_group, key) + HDF5.attributes(d_group)["type"] = "Explicit" #TODO: consider "Indexed" type + d_group["values"] = getfield(motion, Symbol(key)) + end +end \ No newline at end of file diff --git a/KomaMRIFiles/test/runtests.jl b/KomaMRIFiles/test/runtests.jl index a658fbd3e..deccbf82d 100644 --- a/KomaMRIFiles/test/runtests.jl +++ b/KomaMRIFiles/test/runtests.jl @@ -1,6 +1,6 @@ using TestItems, TestItemRunner -@run_package_tests filter=ti->!(:skipci in ti.tags)&&(:files in ti.tags) #verbose=true +@run_package_tests filter=t_start->!(:skipci in t_start.tags)&&(:files in t_start.tags) #verbose=true @testitem "Files" tags=[:files] begin using Suppressor @@ -40,7 +40,7 @@ using TestItems, TestItemRunner obj = read_phantom_jemris(path*"/test_files/column1d.h5") @test obj.name == "column1d.h5" end - # Test JEMRIS + # Test MRiLab @testset "MRiLab" begin path = @__DIR__ filename = path * "/test_files/brain_mrilab.mat" @@ -48,6 +48,49 @@ using TestItems, TestItemRunner obj = read_phantom_MRiLab(filename; FRange_filename) @test obj.name == "brain_mrilab.mat" end + # Test Phantom (.phantom) + @testset "Phantom" begin + using KomaMRIBase + path = @__DIR__ + # NoMotion + filename = path * "/test_files/brain_nomotion.phantom" + obj1 = brain_phantom2D() + write_phantom(obj1, filename) + obj2 = read_phantom(filename) + @test obj1 == obj2 + # SimpleMotion + filename = path * "/test_files/brain_simplemotion.phantom" + obj1 = brain_phantom2D() + obj1.motion = SimpleMotion([ + PeriodicRotation( + period=1.0, + yaw=45.0, + pitch=0.0, + roll=0.0), + Translation( + t_start=0.0, + t_end=0.5, + dx=0.0, + dy=0.02, + dz=0.0 + )]) + write_phantom(obj1, filename) + obj2 = read_phantom(filename) + @test obj1 == obj2 + # ArbitraryMotion + filename = path * "/test_files/brain_arbitrarymotion.phantom" + obj1 = brain_phantom2D() + Ns = length(obj1) + K = 10 + obj1.motion = ArbitraryMotion( + [1.0], + 0.01.*rand(Ns, K-1), + 0.01.*rand(Ns, K-1), + 0.01.*rand(Ns, K-1)) + write_phantom(obj1, filename) + obj2 = read_phantom(filename) + @test obj1 == obj2 + end end @testitem "Pulseq compat" tags=[:files, :pulseq] begin diff --git a/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom b/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom new file mode 100644 index 000000000..e33788427 Binary files /dev/null and b/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom differ diff --git a/KomaMRIFiles/test/test_files/brain_nomotion.phantom b/KomaMRIFiles/test/test_files/brain_nomotion.phantom new file mode 100644 index 000000000..a487a64b8 Binary files /dev/null and b/KomaMRIFiles/test/test_files/brain_nomotion.phantom differ diff --git a/KomaMRIFiles/test/test_files/brain_simplemotion.phantom b/KomaMRIFiles/test/test_files/brain_simplemotion.phantom new file mode 100644 index 000000000..c72c99584 Binary files /dev/null and b/KomaMRIFiles/test/test_files/brain_simplemotion.phantom differ diff --git a/KomaMRIPlots/src/KomaMRIPlots.jl b/KomaMRIPlots/src/KomaMRIPlots.jl index 754b2aa29..8594482e9 100644 --- a/KomaMRIPlots/src/KomaMRIPlots.jl +++ b/KomaMRIPlots/src/KomaMRIPlots.jl @@ -9,11 +9,23 @@ include("ui/DisplayFunctions.jl") using Reexport @reexport using PlotlyJS: savefig -export plot_seq, plot_M0, plot_M1, plot_M2, plot_eddy_currents, plot_seqd, - plot_slew_rate, plot_kspace, plot_phantom_map, plot_signal, plot_image, plot_dict +export plot_seq, + plot_M0, + plot_M1, + plot_M2, + plot_eddy_currents, + plot_seqd, + plot_slew_rate, + plot_kspace, + plot_phantom_map, + plot_signal, + plot_image, + plot_dict #Package version, KomaMRIPlots.__VERSION__ using Pkg -__VERSION__ = VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]) +__VERSION__ = VersionNumber( + Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"] +) end diff --git a/KomaMRIPlots/src/ui/DisplayFunctions.jl b/KomaMRIPlots/src/ui/DisplayFunctions.jl index 0741b4a43..4691b4256 100644 --- a/KomaMRIPlots/src/ui/DisplayFunctions.jl +++ b/KomaMRIPlots/src/ui/DisplayFunctions.jl @@ -14,74 +14,94 @@ Define colors for dark or light mode. - `sep_color`: (`::String`) color of separator lines """ function theme_chooser(darkmode) - if darkmode - bgcolor = "rgba(0,0,0,0)"#"rgb(13,16,17)" - text_color = "gray" - plot_bgcolor = "rgb(22,26,29)" #"rgb(33,37,41)" - grid_color = "rgb(40,52,66)" #rgb(40,40,40) - sep_color = "white" - else - bgcolor = "rgba(0,0,0,0)"#"white" - text_color = "gray"#"rgb(49,70,101)" - plot_bgcolor = "rgb(229,236,246)" - grid_color = "white" - sep_color = "black" - end - bgcolor, text_color, plot_bgcolor, grid_color, sep_color + if darkmode + bgcolor = "rgba(0,0,0,0)"#"rgb(13,16,17)" + text_color = "gray" + plot_bgcolor = "rgb(22,26,29)" #"rgb(33,37,41)" + grid_color = "rgb(40,52,66)" #rgb(40,40,40) + sep_color = "white" + else + bgcolor = "rgba(0,0,0,0)"#"white" + text_color = "gray"#"rgb(49,70,101)" + plot_bgcolor = "rgb(229,236,246)" + grid_color = "white" + sep_color = "black" + end + return bgcolor, text_color, plot_bgcolor, grid_color, sep_color end -function generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0) - #LAYOUT - bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) - l = Layout(;title=title, hovermode="closest", - xaxis_title="", - modebar=attr(orientation="h", yanchor="bottom", xanchor="right", y=1, x=0, bgcolor=bgcolor, color=text_color, activecolor=plot_bgcolor), - legend=attr(orientation="h", yanchor="bottom", xanchor="left", y=1, x=0), - plot_bgcolor=plot_bgcolor, - paper_bgcolor=bgcolor, - xaxis_gridcolor=grid_color, - yaxis_gridcolor=grid_color, - xaxis_zerolinecolor=grid_color, - yaxis_zerolinecolor=grid_color, - font_color=text_color, - yaxis_fixedrange = false, - xaxis=attr( - ticksuffix=" ms", - domain=range[:], - range=range[:], - rangeslider=attr(visible=slider), - rangeselector=attr( - buttons=[ - attr(count=1, - label="1m", - step=10, - stepmode="backward"), - attr(step="all") - ] - ), - ), - margin=attr(t=0,l=0,r=0,b=0) - ) +function generate_seq_time_layout_config( + title, width, height, range, slider, show_seq_blocks, darkmode; T0 +) + #LAYOUT + bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) + l = Layout(; + title=title, + hovermode="closest", + xaxis_title="", + modebar=attr(; + orientation="h", + yanchor="bottom", + xanchor="right", + y=1, + x=0, + bgcolor=bgcolor, + color=text_color, + activecolor=plot_bgcolor, + ), + legend=attr(; orientation="h", yanchor="bottom", xanchor="left", y=1, x=0), + plot_bgcolor=plot_bgcolor, + paper_bgcolor=bgcolor, + xaxis_gridcolor=grid_color, + yaxis_gridcolor=grid_color, + xaxis_zerolinecolor=grid_color, + yaxis_zerolinecolor=grid_color, + font_color=text_color, + yaxis_fixedrange=false, + xaxis=attr(; + ticksuffix=" ms", + domain=range[:], + range=range[:], + rangeslider=attr(; visible=slider), + rangeselector=attr(; + buttons=[ + attr(; count=1, label="1m", step=10, stepmode="backward"), + attr(; step="all"), + ], + ), + ), + margin=attr(; t=0, l=0, r=0, b=0), + ) if show_seq_blocks l.xaxis["tickvals"] = T0 * 1e3 end - if height !== nothing - l.height = height - end - if width !== nothing - l.width = width - end - #CONFIG - config = PlotConfig( - displaylogo=false, - toImageButtonOptions=attr( - format="svg", # one of png, svg, jpeg, webp - ).fields, - modeBarButtonsToRemove=["zoom", "select2d", "lasso2d", "autoScale", "resetScale2d", "pan", - "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"] - ) - - l, config + if height !== nothing + l.height = height + end + if width !== nothing + l.width = width + end + #CONFIG + config = PlotConfig(; + displaylogo=false, + toImageButtonOptions=attr(; + format="svg", # one of png, svg, jpeg, webp + ).fields, + modeBarButtonsToRemove=[ + "zoom", + "select2d", + "lasso2d", + "autoScale", + "resetScale2d", + "pan", + "tableRotation", + "resetCameraLastSave", + "zoomIn", + "zoomOut", + ], + ) + + return l, config end """ @@ -102,15 +122,15 @@ Interpolates a color map. This is used for plotting the kspace (refer to values """ function interp_map(c_map, t_interp) - idx = [c[1] for c = c_map] - R = [parse.(Int, split(c[2][5:end-1],", "))[1] for c = c_map] - G = [parse.(Int, split(c[2][5:end-1],", "))[2] for c = c_map] - B = [parse.(Int, split(c[2][5:end-1],", "))[3] for c = c_map] - r = linear_interpolation(idx,R)(t_interp) - g = linear_interpolation(idx,G)(t_interp) - b = linear_interpolation(idx,B)(t_interp) - c_map_interp = ["hsv($r, $g, $b)" for (r,g,b)=zip(r,g,b)] - c_map_interp + idx = [c[1] for c in c_map] + R = [parse.(Int, split(c[2][5:(end - 1)], ", "))[1] for c in c_map] + G = [parse.(Int, split(c[2][5:(end - 1)], ", "))[2] for c in c_map] + B = [parse.(Int, split(c[2][5:(end - 1)], ", "))[3] for c in c_map] + r = linear_interpolation(idx, R)(t_interp) + g = linear_interpolation(idx, G)(t_interp) + b = linear_interpolation(idx, B)(t_interp) + c_map_interp = ["hsv($r, $g, $b)" for (r, g, b) in zip(r, g, b)] + return c_map_interp end """ @@ -147,23 +167,23 @@ julia> plot_seq(seq) ``` """ function plot_seq( - seq::Sequence; - width=nothing, - height=nothing, - slider=false, - show_seq_blocks=false, - darkmode=false, - range=[], - title="", - xaxis="x", - yaxis="y", - showlegend=true, - freq_in_phase=false, - # Performance related - gl=false, - max_rf_samples=100, - show_adc=false - ) + seq::Sequence; + width=nothing, + height=nothing, + slider=false, + show_seq_blocks=false, + darkmode=false, + range=[], + title="", + xaxis="x", + yaxis="y", + showlegend=true, + freq_in_phase=false, + # Performance related + gl=false, + max_rf_samples=100, + show_adc=false, +) # Aux functions scatter_fun = gl ? scattergl : scatter @@ -347,33 +367,63 @@ julia> plot_M0(seq) ``` """ function plot_M0( - seq::Sequence; - width=nothing, - height=nothing, - slider=true, - show_seq_blocks=false, - darkmode=false, - range=[], - title="", - skip_rf=zeros(Bool, sum(is_RF_on.(seq))) - ) - #Times - t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) - t = t[1:end-1] - T0 = get_block_start_times(seq) - #M0 - ts = t .+ Δt - rf_idx, rf_type = KomaMRIBase.get_RF_types(seq, t) - k, _ = KomaMRIBase.get_kspace(seq; Δt=1, skip_rf) - #plots M0 - p = [scatter() for j=1:4] - p[1] = scatter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", name="M0x", legendgroup="Gx", marker=attr(color="#636EFA")) - p[2] = scatter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", name="M0y", legendgroup="Gy", marker=attr(color="#EF553B")) - p[3] = scatter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", name="M0z", legendgroup="Gz", marker=attr(color="#00CC96")) - p[4] = scatter(x=t[rf_idx]*1e3,y=rf_type,name="RFs",marker=attr(symbol="cross",size=8,color="orange"), mode="markers", showlegend=false) - #Layout and config - l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0) - return plot_koma(p, l; config) + seq::Sequence; + width=nothing, + height=nothing, + slider=true, + show_seq_blocks=false, + darkmode=false, + range=[], + title="", + skip_rf=zeros(Bool, sum(is_RF_on.(seq))), +) + #Times + t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) + t = t[1:(end - 1)] + T0 = get_block_start_times(seq) + #M0 + ts = t .+ Δt + rf_idx, rf_type = KomaMRIBase.get_RF_types(seq, t) + k, _ = KomaMRIBase.get_kspace(seq; Δt=1, skip_rf) + #plots M0 + p = [scatter() for j in 1:4] + p[1] = scatter(; + x=ts * 1e3, + y=k[:, 1], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", + name="M0x", + legendgroup="Gx", + marker=attr(; color="#636EFA"), + ) + p[2] = scatter(; + x=ts * 1e3, + y=k[:, 2], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", + name="M0y", + legendgroup="Gy", + marker=attr(; color="#EF553B"), + ) + p[3] = scatter(; + x=ts * 1e3, + y=k[:, 3], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms)", + name="M0z", + legendgroup="Gz", + marker=attr(; color="#00CC96"), + ) + p[4] = scatter(; + x=t[rf_idx] * 1e3, + y=rf_type, + name="RFs", + marker=attr(; symbol="cross", size=8, color="orange"), + mode="markers", + showlegend=false, + ) + #Layout and config + l, config = generate_seq_time_layout_config( + title, width, height, range, slider, show_seq_blocks, darkmode; T0 + ) + return plot_koma(p, l; config) end """ @@ -406,36 +456,65 @@ julia> plot_M1(seq) ``` """ function plot_M1( - seq::Sequence; - width=nothing, - height=nothing, - slider=true, - show_seq_blocks=false, - darkmode=false, - range=[], - title="", - skip_rf=zeros(Bool, sum(is_RF_on.(seq))) - ) - #Times - t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) - t = t[1:end-1] - T0 = get_block_start_times(seq) - #M1 - ts = t .+ Δt - rf_idx, rf_type = KomaMRIBase.get_RF_types(seq, t) - k, _ = KomaMRIBase.get_M1(seq; Δt=1, skip_rf) - #plots M1 - p = [scatter() for j=1:4] - p[1] = scatter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms²)", name="M1x", legendgroup="Gx", marker=attr(color="#636EFA")) - p[2] = scatter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms²)", name="M1y", legendgroup="Gy", marker=attr(color="#EF553B")) - p[3] = scatter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms²)", name="M1z", legendgroup="Gz", marker=attr(color="#00CC96")) - p[4] = scatter(x=t[rf_idx]*1e3,y=rf_type,name="RFs",marker=attr(symbol="cross",size=8,color="orange"), mode="markers", showlegend=false) - #Layout and config - l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0) - return plot_koma(p, l; config) + seq::Sequence; + width=nothing, + height=nothing, + slider=true, + show_seq_blocks=false, + darkmode=false, + range=[], + title="", + skip_rf=zeros(Bool, sum(is_RF_on.(seq))), +) + #Times + t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) + t = t[1:(end - 1)] + T0 = get_block_start_times(seq) + #M1 + ts = t .+ Δt + rf_idx, rf_type = KomaMRIBase.get_RF_types(seq, t) + k, _ = KomaMRIBase.get_M1(seq; Δt=1, skip_rf) + #plots M1 + p = [scatter() for j in 1:4] + p[1] = scatter(; + x=ts * 1e3, + y=k[:, 1], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms²)", + name="M1x", + legendgroup="Gx", + marker=attr(; color="#636EFA"), + ) + p[2] = scatter(; + x=ts * 1e3, + y=k[:, 2], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms²)", + name="M1y", + legendgroup="Gy", + marker=attr(; color="#EF553B"), + ) + p[3] = scatter(; + x=ts * 1e3, + y=k[:, 3], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms²)", + name="M1z", + legendgroup="Gz", + marker=attr(; color="#00CC96"), + ) + p[4] = scatter(; + x=t[rf_idx] * 1e3, + y=rf_type, + name="RFs", + marker=attr(; symbol="cross", size=8, color="orange"), + mode="markers", + showlegend=false, + ) + #Layout and config + l, config = generate_seq_time_layout_config( + title, width, height, range, slider, show_seq_blocks, darkmode; T0 + ) + return plot_koma(p, l; config) end - """ p = plot_M2(seq::Sequence; kwargs...) @@ -466,35 +545,64 @@ julia> plot_M2(seq) ``` """ function plot_M2( - seq::Sequence; - width = nothing, - height = nothing, - slider = true, - show_seq_blocks = false, - darkmode = false, - range = [], - title = "" - ) - #Times - t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) - t = t[1:end-1] - T0 = get_block_start_times(seq) - #M2 - ts = t .+ Δt - rf_idx, rf_type = KomaMRIBase.get_RF_types(seq, t) - k, _ = KomaMRIBase.get_M2(seq; Δt=1) - #Plor M2 - p = [scatter() for j=1:4] - p[1] = scatter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2x", legendgroup="Gx", marker=attr(color="#636EFA")) - p[2] = scatter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2y", legendgroup="Gy", marker=attr(color="#EF553B")) - p[3] = scatter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2z", legendgroup="Gz", marker=attr(color="#00CC96")) - p[4] = scatter(x=t[rf_idx]*1e3,y=rf_type,name="RFs",marker=attr(symbol="cross",size=8,color="orange"), mode="markers", showlegend=false) - #Layout and config - l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0) - return plot_koma(p, l; config) + seq::Sequence; + width=nothing, + height=nothing, + slider=true, + show_seq_blocks=false, + darkmode=false, + range=[], + title="", +) + #Times + t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) + t = t[1:(end - 1)] + T0 = get_block_start_times(seq) + #M2 + ts = t .+ Δt + rf_idx, rf_type = KomaMRIBase.get_RF_types(seq, t) + k, _ = KomaMRIBase.get_M2(seq; Δt=1) + #Plor M2 + p = [scatter() for j in 1:4] + p[1] = scatter(; + x=ts * 1e3, + y=k[:, 1], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", + name="M2x", + legendgroup="Gx", + marker=attr(; color="#636EFA"), + ) + p[2] = scatter(; + x=ts * 1e3, + y=k[:, 2], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", + name="M2y", + legendgroup="Gy", + marker=attr(; color="#EF553B"), + ) + p[3] = scatter(; + x=ts * 1e3, + y=k[:, 3], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", + name="M2z", + legendgroup="Gz", + marker=attr(; color="#00CC96"), + ) + p[4] = scatter(; + x=t[rf_idx] * 1e3, + y=rf_type, + name="RFs", + marker=attr(; symbol="cross", size=8, color="orange"), + mode="markers", + showlegend=false, + ) + #Layout and config + l, config = generate_seq_time_layout_config( + title, width, height, range, slider, show_seq_blocks, darkmode; T0 + ) + return plot_koma(p, l; config) end - """ p = plot_eddy_currents(seq::Sequence, λ; kwargs...) @@ -527,35 +635,59 @@ julia> plot_eddy_currents(seq, 80e-3) ``` """ function plot_eddy_currents( - seq::Sequence, λ; - α = ones(size(λ)), - width = nothing, - height = nothing, - slider = true, - show_seq_blocks = false, - darkmode = false, - range = [], - title = "" - ) - #Times - t, Δt = KomaMRIBase.get_variable_times(seq + ADC(100, 100e-3); Δt=1) - t = t[2:end] - T0 = get_block_start_times(seq) - Gx, Gy, Gz = KomaMRIBase.get_grads(seq, t) - #Eddy currents per lambda - Gec = zeros(length(t), 3) - for (i, l) in enumerate(λ) - aux, _ = KomaMRIBase.get_eddy_currents(seq + ADC(100, 100e-3); Δt=1, λ=l) - Gec .+= α[i] .* aux - end - #Plot eddy currents - p = [scatter() for j=1:4] - p[1] = scatter(x=t*1e3, y=(Gx*0 .+ Gec[:,1])*1e3, hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", name="ECx", legendgroup="Gx", marker=attr(color="#636EFA")) - p[2] = scatter(x=t*1e3, y=(Gy*0 .+ Gec[:,2])*1e3, hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", name="ECy", legendgroup="Gy", marker=attr(color="#EF553B")) - p[3] = scatter(x=t*1e3, y=(Gz*0 .+ Gec[:,3])*1e3, hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", name="ECz", legendgroup="Gz", marker=attr(color="#00CC96")) - #Layout and config - l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0) - return plot_koma(p, l; config) + seq::Sequence, + λ; + α=ones(size(λ)), + width=nothing, + height=nothing, + slider=true, + show_seq_blocks=false, + darkmode=false, + range=[], + title="", +) + #Times + t, Δt = KomaMRIBase.get_variable_times(seq + ADC(100, 100e-3); Δt=1) + t = t[2:end] + T0 = get_block_start_times(seq) + Gx, Gy, Gz = KomaMRIBase.get_grads(seq, t) + #Eddy currents per lambda + Gec = zeros(length(t), 3) + for (i, l) in enumerate(λ) + aux, _ = KomaMRIBase.get_eddy_currents(seq + ADC(100, 100e-3); Δt=1, λ=l) + Gec .+= α[i] .* aux + end + #Plot eddy currents + p = [scatter() for j in 1:4] + p[1] = scatter(; + x=t * 1e3, + y=(Gx * 0 .+ Gec[:, 1]) * 1e3, + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", + name="ECx", + legendgroup="Gx", + marker=attr(; color="#636EFA"), + ) + p[2] = scatter(; + x=t * 1e3, + y=(Gy * 0 .+ Gec[:, 2]) * 1e3, + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", + name="ECy", + legendgroup="Gy", + marker=attr(; color="#EF553B"), + ) + p[3] = scatter(; + x=t * 1e3, + y=(Gz * 0 .+ Gec[:, 3]) * 1e3, + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m)", + name="ECz", + legendgroup="Gz", + marker=attr(; color="#00CC96"), + ) + #Layout and config + l, config = generate_seq_time_layout_config( + title, width, height, range, slider, show_seq_blocks, darkmode; T0 + ) + return plot_koma(p, l; config) end """ @@ -588,33 +720,55 @@ julia> plot_slew_rate(seq) ``` """ function plot_slew_rate( - seq::Sequence; - width = nothing, - height = nothing, - slider = true, - show_seq_blocks = false, - darkmode = false, - range = [], - title = "" - ) - #Times - t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) - t = t[1:end-1] - T0 = get_block_start_times(seq) - ts = t .+ Δt - #Eddy currents per lambda - k, _ = KomaMRIBase.get_slew_rate(seq; Δt=1) - #Plot eddy currents - p = [scatter() for j=1:4] - p[1] = scatter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", name="SRx", legendgroup="Gx", marker=attr(color="#636EFA")) - p[2] = scatter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", name="SRy", legendgroup="Gy", marker=attr(color="#EF553B")) - p[3] = scatter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", name="SRz", legendgroup="Gz", marker=attr(color="#00CC96")) - #Layout and config - l, config = generate_seq_time_layout_config(title, width, height, range, slider, show_seq_blocks, darkmode; T0) - return plot_koma(p, l; config) + seq::Sequence; + width=nothing, + height=nothing, + slider=true, + show_seq_blocks=false, + darkmode=false, + range=[], + title="", +) + #Times + t, Δt = KomaMRIBase.get_variable_times(seq; Δt=1) + t = t[1:(end - 1)] + T0 = get_block_start_times(seq) + ts = t .+ Δt + #Eddy currents per lambda + k, _ = KomaMRIBase.get_slew_rate(seq; Δt=1) + #Plot eddy currents + p = [scatter() for j in 1:4] + p[1] = scatter(; + x=ts * 1e3, + y=k[:, 1], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", + name="SRx", + legendgroup="Gx", + marker=attr(; color="#636EFA"), + ) + p[2] = scatter(; + x=ts * 1e3, + y=k[:, 2], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", + name="SRy", + legendgroup="Gy", + marker=attr(; color="#EF553B"), + ) + p[3] = scatter(; + x=ts * 1e3, + y=k[:, 3], + hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m/ms)", + name="SRz", + legendgroup="Gz", + marker=attr(; color="#00CC96"), + ) + #Layout and config + l, config = generate_seq_time_layout_config( + title, width, height, range, slider, show_seq_blocks, darkmode; T0 + ) + return plot_koma(p, l; config) end - """ p = plot_image(image; height, width, zmin, zmax, darkmode, title) @@ -635,43 +789,60 @@ Plots an image matrix. - `p`: (`::PlotlyJS.SyncPlot`) plot of the image matrix """ function plot_image( - image; - height = 600, - width = nothing, - zmin = minimum(abs.(image[:])), - zmax = maximum(abs.(image[:])), - darkmode = false, - title = "" - ) - #Layout - bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) - l = Layout(;title=title,yaxis_title="y", - xaxis_title="x",margin=attr(t=50,l=0,r=0,b=0), - yaxis=attr(scaleanchor="x"), - font_color=text_color, - modebar=attr(orientation="v",bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),xaxis=attr(constrain="domain"),hovermode="closest", - paper_bgcolor=bgcolor, - plot_bgcolor=plot_bgcolor, - xaxis_gridcolor=grid_color, - yaxis_gridcolor=grid_color, - xaxis_zerolinecolor=grid_color, - yaxis_zerolinecolor=grid_color) + image; + height=600, + width=nothing, + zmin=minimum(abs.(image[:])), + zmax=maximum(abs.(image[:])), + darkmode=false, + title="", +) + #Layout + bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) + l = Layout(; + title=title, + yaxis_title="y", + xaxis_title="x", + margin=attr(; t=50, l=0, r=0, b=0), + yaxis=attr(; scaleanchor="x"), + font_color=text_color, + modebar=attr(; + orientation="v", bgcolor=bgcolor, color=text_color, activecolor=plot_bgcolor + ), + xaxis=attr(; constrain="domain"), + hovermode="closest", + paper_bgcolor=bgcolor, + plot_bgcolor=plot_bgcolor, + xaxis_gridcolor=grid_color, + yaxis_gridcolor=grid_color, + xaxis_zerolinecolor=grid_color, + yaxis_zerolinecolor=grid_color, + ) if height !== nothing l.height = height end if width !== nothing l.width = width end - #Plot - p = heatmap(z=image,transpose=false,zmin=zmin,zmax=zmax,colorscale="Greys") - config = PlotConfig( - displaylogo=false, - toImageButtonOptions=attr( - format="svg", # one of png, svg, jpeg, webp - ).fields, - modeBarButtonsToRemove=["zoom", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"] - ) - return plot_koma(p, l; config) + #Plot + p = heatmap(; z=image, transpose=false, zmin=zmin, zmax=zmax, colorscale="Greys") + config = PlotConfig(; + displaylogo=false, + toImageButtonOptions=attr(; + format="svg", # one of png, svg, jpeg, webp + ).fields, + modeBarButtonsToRemove=[ + "zoom", + "autoScale", + "resetScale2d", + "pan", + "tableRotation", + "resetCameraLastSave", + "zoomIn", + "zoomOut", + ], + ) + return plot_koma(p, l; config) end """ @@ -699,72 +870,121 @@ julia> seq = read_seq(seq_file) julia> plot_kspace(seq) ``` """ -function plot_kspace( - seq::Sequence; - width=nothing, - height=nothing, - darkmode=false - ) - bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) - #Calculations of theoretical k-space - kspace, kspace_adc = get_kspace(seq; Δt=1) #sim_params["Δt"]) - t_adc = KomaMRIBase.get_adc_sampling_times(seq) - #Colormap - c_map = [[t, "hsv($(floor(Int,(1-t)*255)), 100, 50)"] for t=range(0,1;length=10)] # range(s,b,N) only works in Julia 1.7.3 - c = "gray" - c2_idx = [] - counter = 0 - for s in seq - if is_ADC_on(s) - N = s.ADC.N[1] - append!(c2_idx, counter:N+counter-1) - counter += N - end - end - c2 = interp_map(c_map, c2_idx ./ maximum(c2_idx)) - #Layout - mink = minimum(kspace_adc,dims=1) - maxk = maximum(kspace_adc,dims=1) - dW = maximum(maxk .- mink, dims=2) * .3 - mink .-= dW - maxk .+= dW - #Layout - l = Layout(; - paper_bgcolor=bgcolor, - scene=attr(xaxis=attr(title="kx [m⁻¹]",range=[mink[1],maxk[1]],backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color), - yaxis=attr(title="ky [m⁻¹]",range=[mink[2],maxk[2]],backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color), - zaxis=attr(title="kz [m⁻¹]",range=[mink[3],maxk[3]],backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color)), - modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor), - legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0), - font_color=text_color, - scene_camera_eye=attr(x=0, y=0, z=1.7), - scene_camera_up=attr(x=0, y=1., z=0), - scene_aspectmode="cube", - margin=attr(t=0,l=0,r=0)) +function plot_kspace(seq::Sequence; width=nothing, height=nothing, darkmode=false) + bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) + #Calculations of theoretical k-space + kspace, kspace_adc = get_kspace(seq; Δt=1) #sim_params["Δt"]) + t_adc = KomaMRIBase.get_adc_sampling_times(seq) + #Colormap + c_map = [[t, "hsv($(floor(Int,(1-t)*255)), 100, 50)"] for t in range(0, 1; length=10)] # range(s,b,N) only works in Julia 1.7.3 + c = "gray" + c2_idx = [] + counter = 0 + for s in seq + if is_ADC_on(s) + N = s.ADC.N[1] + append!(c2_idx, counter:(N + counter - 1)) + counter += N + end + end + c2 = interp_map(c_map, c2_idx ./ maximum(c2_idx)) + #Layout + mink = minimum(kspace_adc; dims=1) + maxk = maximum(kspace_adc; dims=1) + dW = maximum(maxk .- mink; dims=2) * 0.3 + mink .-= dW + maxk .+= dW + #Layout + l = Layout(; + paper_bgcolor=bgcolor, + scene=attr(; + xaxis=attr(; + title="kx [m⁻¹]", + range=[mink[1], maxk[1]], + backgroundcolor=plot_bgcolor, + gridcolor=grid_color, + zerolinecolor=grid_color, + ), + yaxis=attr(; + title="ky [m⁻¹]", + range=[mink[2], maxk[2]], + backgroundcolor=plot_bgcolor, + gridcolor=grid_color, + zerolinecolor=grid_color, + ), + zaxis=attr(; + title="kz [m⁻¹]", + range=[mink[3], maxk[3]], + backgroundcolor=plot_bgcolor, + gridcolor=grid_color, + zerolinecolor=grid_color, + ), + ), + modebar=attr(; + orientation="h", + yanchor="bottom", + xanchor="right", + y=1, + x=0, + bgcolor=bgcolor, + color=text_color, + activecolor=plot_bgcolor, + ), + legend=attr(; orientation="h", yanchor="bottom", xanchor="left", y=1, x=0), + font_color=text_color, + scene_camera_eye=attr(; x=0, y=0, z=1.7), + scene_camera_up=attr(; x=0, y=1.0, z=0), + scene_aspectmode="cube", + margin=attr(; t=0, l=0, r=0), + ) if height !== nothing l.height = height end if width !== nothing l.width = width end - #Plot - p = [scatter() for j=1:3] - p[1] = scatter3d(x=kspace[:,1],y=kspace[:,2],z=kspace[:,3],mode="lines", - line=attr(color=c),name="Trajectory",hoverinfo="skip") - p[2] = scatter3d(x=kspace_adc[:,1],y=kspace_adc[:,2],z=kspace_adc[:,3],text=round.(t_adc*1e3,digits=3),mode="markers", - line=attr(color=c2),marker=attr(size=2),name="ADC",hovertemplate="kx: %{x:.1f} m⁻¹
ky: %{y:.1f} m⁻¹
kz: %{z:.1f} m⁻¹
t_acq: %{text} ms") - p[3] = scatter3d(x=[0],y=[0],z=[0],name="k=0",marker=attr(symbol="cross",size=10,color="red")) - config = PlotConfig( - displaylogo=false, - toImageButtonOptions=attr( - format="svg", # one of png, svg, jpeg, webp - ).fields, - modeBarButtonsToRemove=["zoom", "pan", "tableRotation", "resetCameraLastSave3d", "orbitRotation", "resetCameraDefault3d"] - ) - return plot_koma(p, l; config) + #Plot + p = [scatter() for j in 1:3] + p[1] = scatter3d(; + x=kspace[:, 1], + y=kspace[:, 2], + z=kspace[:, 3], + mode="lines", + line=attr(; color=c), + name="Trajectory", + hoverinfo="skip", + ) + p[2] = scatter3d(; + x=kspace_adc[:, 1], + y=kspace_adc[:, 2], + z=kspace_adc[:, 3], + text=round.(t_adc * 1e3, digits=3), + mode="markers", + line=attr(; color=c2), + marker=attr(; size=2), + name="ADC", + hovertemplate="kx: %{x:.1f} m⁻¹
ky: %{y:.1f} m⁻¹
kz: %{z:.1f} m⁻¹
t_acq: %{text} ms", + ) + p[3] = scatter3d(; + x=[0], y=[0], z=[0], name="k=0", marker=attr(; symbol="cross", size=10, color="red") + ) + config = PlotConfig(; + displaylogo=false, + toImageButtonOptions=attr(; + format="svg", # one of png, svg, jpeg, webp + ).fields, + modeBarButtonsToRemove=[ + "zoom", + "pan", + "tableRotation", + "resetCameraLastSave3d", + "orbitRotation", + "resetCameraDefault3d", + ], + ) + return plot_koma(p, l; config) end - """ p = plot_phantom_map(obj::Phantom, key::Symbol; kwargs...) @@ -776,7 +996,6 @@ Plots a phantom map for a specific spin parameter given by `key`. displaying different parameters of the phantom spins # Keywords -- `t0`: (`::Real`, `=0`, `[ms]`) time to see displacement of the phantom - `height`: (`::Integer`, `=600`) plot height - `width`: (`::Integer`, `=nothing`) plot width - `darkmode`: (`::Bool`, `=false`) boolean to indicate whether to display darkmode style @@ -800,123 +1019,336 @@ julia> plot_phantom_map(obj3D, :ρ) ``` """ function plot_phantom_map( - ph::Phantom, - key::Symbol; - t0=0, - height=600, - width=nothing, - darkmode=false, - view_2d=false, - colorbar=true, - kwargs... - ) - path = @__DIR__ - cmin_key = minimum(getproperty(ph,key)) - cmax_key = maximum(getproperty(ph,key)) - if key == :T1 || key == :T2 || key == :T2s - cmin_key = 0 - factor = 1e3 - unit = " ms" - if key == :T1 - cmax_key = 2500/factor - colors = MAT.matread(path*"/assets/T1cm.mat")["T1colormap"] - N, _ = size(colors) - idx = range(0,1;length=N) #range(0,T,N) works in Julia 1.7 - colormap = [[idx[n], "rgb($(floor(Int,colors[n,1]*255)),$(floor(Int,colors[n,2]*255)),$(floor(Int,colors[n,3]*255)))"] for n=1:N] - elseif key == :T2 || key == :T2s - if key == :T2 - cmax_key = 250/factor - end - colors = MAT.matread(path*"/assets/T2cm.mat")["T2colormap"] - N, _ = size(colors) - idx = range(0,1;length=N) #range(0,T,N) works in Julia 1.7 - colormap = [[idx[n], "rgb($(floor(Int,colors[n,1]*255)),$(floor(Int,colors[n,2]*255)),$(floor(Int,colors[n,3]*255)))"] for n=1:N] - end - elseif key == :x || key == :y || key == :z - factor = 1e2 - unit = " cm" - colormap="Greys" - elseif key == :Δw - factor = 1/(2π) - unit = " Hz" - colormap="Greys" - else - factor = 1 - cmin_key = 0 - unit="" - colormap="Greys" - end - cmin_key = get(kwargs, :cmin, factor * cmin_key) - cmax_key = get(kwargs, :cmax, factor * cmax_key) - x0 = -maximum(abs.([ph.x ph.y ph.z]))*1e2 - xf = maximum(abs.([ph.x ph.y ph.z]))*1e2 - #Layout - bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) - l = Layout(;title=ph.name*": "*string(key), - xaxis_title="x", - yaxis_title="y", - plot_bgcolor=plot_bgcolor, - paper_bgcolor=bgcolor, - xaxis_gridcolor=grid_color, - yaxis_gridcolor=grid_color, - xaxis_zerolinecolor=grid_color, - yaxis_zerolinecolor=grid_color, - font_color=text_color, - scene=attr( - xaxis=attr(title="x",range=[x0,xf],ticksuffix=" cm",backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color), - yaxis=attr(title="y",range=[x0,xf],ticksuffix=" cm",backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color), - zaxis=attr(title="z",range=[x0,xf],ticksuffix=" cm",backgroundcolor=plot_bgcolor,gridcolor=grid_color,zerolinecolor=grid_color), - aspectmode="manual", - aspectratio=attr(x=1,y=1,z=1)), - margin=attr(t=50,l=0,r=0), - modebar=attr(orientation="h",bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor), - xaxis=attr(constrain="domain"), - yaxis=attr(scaleanchor="x"), - hovermode="closest") + ph::Phantom, + key::Symbol; + height=700, + width=nothing, + darkmode=false, + view_2d=sum(KomaMRIBase.get_dims(ph)) < 3, + colorbar=true, + intermediate_time_samples=0, + max_time_samples=100, + max_spins=50000, + dt_frame=250, + kwargs..., +) + function process_times(motion::SimpleMotion) + t = times(motion) + # Interpolate time points (as many as indicated by intermediate_time_samples) + itp = interpolate( + ( + 1:(intermediate_time_samples + 1):(length(t) + intermediate_time_samples * (length(t) - 1)), + ), + t, + Gridded(Linear()), + ) + return itp.(1:(length(t) + intermediate_time_samples * (length(t) - 1))) + end + function process_times(motion::ArbitraryMotion) + t = times(motion) + # Interpolate time points (as many as indicated by intermediate_time_samples) + itp = interpolate( + ( + 1:(intermediate_time_samples + 1):(length(t) + intermediate_time_samples * (length(t) - 1)), + ), + t, + Gridded(Linear()), + ) + t = itp.(1:(length(t) + intermediate_time_samples * (length(t) - 1))) + # Decimate time points so their number is smaller than max_time_samples + ss = length(t) > max_time_samples ? length(t) ÷ max_time_samples : 1 + return t[1:ss:end] + end + process_times(motion::MotionModel) = times(motion) + + # IDEA: subsample phantoms which are too large + function decimate_uniform_phantom(ph::Phantom, num_points::Int) + dimx, dimy, dimz = KomaMRIBase.get_dims(ph) + ss = Int(ceil((length(ph) / num_points)^(1 / sum(KomaMRIBase.get_dims(ph))))) + ssx = dimx ? ss : 1 + ssy = dimy ? ss : 1 + ssz = dimz ? ss : 1 + ix = sortperm(ph.x)[1:ssx:end] + iy = sortperm(ph.y)[1:ssy:end] + iz = sortperm(ph.z)[1:ssz:end] + idx = intersect(ix, iy, iz) + return ph[idx] + end + + n_spins_before_decimate = length(ph) + ph = decimate_uniform_phantom(ph, max_spins) + n_spins_after_decimate = length(ph) + if n_spins_before_decimate > n_spins_after_decimate + @warn "Only $(n_spins_after_decimate) (approximately evenly spaced) of $(n_spins_before_decimate) spins are displayed in order to avoid CPU overload. + This affects display functions but not simulation functions." + end + + path = @__DIR__ + cmin_key = minimum(getproperty(ph, key)) + cmax_key = maximum(getproperty(ph, key)) + if key == :T1 || key == :T2 || key == :T2s + cmin_key = 0 + factor = 1e3 + unit = " ms" + if key == :T1 + cmax_key = 2500 / factor + colors = MAT.matread(path * "/assets/T1cm.mat")["T1colormap"] + N, _ = size(colors) + idx = range(0, 1; length=N) #range(0,T,N) works in Julia 1.7 + colormap = [ + [ + idx[n], + "rgb($(floor(Int,colors[n,1]*255)),$(floor(Int,colors[n,2]*255)),$(floor(Int,colors[n,3]*255)))", + ] for n in 1:N + ] + elseif key == :T2 || key == :T2s + if key == :T2 + cmax_key = 250 / factor + end + colors = MAT.matread(path * "/assets/T2cm.mat")["T2colormap"] + N, _ = size(colors) + idx = range(0, 1; length=N) #range(0,T,N) works in Julia 1.7 + colormap = [ + [ + idx[n], + "rgb($(floor(Int,colors[n,1]*255)),$(floor(Int,colors[n,2]*255)),$(floor(Int,colors[n,3]*255)))", + ] for n in 1:N + ] + end + elseif key == :x || key == :y || key == :z + factor = 1e2 + unit = " cm" + colormap = "Greys" + elseif key == :Δw + factor = 1 / (2π) + unit = " Hz" + colormap = "Greys" + else + factor = 1 + cmin_key = 0 + unit = "" + colormap = "Greys" + end + cmin_key = get(kwargs, :cmin, factor * cmin_key) + cmax_key = get(kwargs, :cmax, factor * cmax_key) + + sort_motions!(ph.motion) + t = process_times(ph.motion) + x, y, z = get_spin_coords(ph.motion, ph.x, ph.y, ph.z, t') + + x0 = -maximum(abs.([x y z])) * 1e2 + xf = maximum(abs.([x y z])) * 1e2 + + if view_2d + trace = [ + scatter(; + x=(x[:, 1]) * 1e2, + y=(y[:, 1]) * 1e2, + mode="markers", + marker=attr(; + color=getproperty(ph, key) * factor, + showscale=colorbar, + colorscale=colormap, + colorbar=attr(; ticksuffix=unit, title=string(key)), + cmin=cmin_key, + cmax=cmax_key, + size=4, + ), + showlegend=false, + text=round.(getproperty(ph, key) * factor, digits=4), + hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
$(string(key)): %{text}$unit", + ), + ] + frames = PlotlyFrame[ + frame(; + data=[ + attr(; + x=(x[:, i]) * 1e2, + y=(y[:, i]) * 1e2, + z=(z[:, i]) * 1e2, + zmin=0, + zmax=1, + ), + ], + name="frame_$i", + traces=[0], + ) for i in 1:length(t) + ] + else + trace = [ + scatter3d(; + x=(x[:, 1]) * 1e2, + y=(y[:, 1]) * 1e2, + z=(z[:, 1]) * 1e2, + mode="markers", + marker=attr(; + color=getproperty(ph, key) * factor, + showscale=colorbar, + colorscale=colormap, + colorbar=attr(; ticksuffix=unit, title=string(key)), + cmin=cmin_key, + cmax=cmax_key, + size=2, + ), + showlegend=false, + text=round.(getproperty(ph, key) * factor, digits=4), + hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
z: %{z:.1f} cm
$(string(key)): %{text}$unit", + ), + ] + frames = PlotlyFrame[ + frame(; + data=[ + attr(; + x=(x[:, i]) * 1e2, + y=(y[:, i]) * 1e2, + z=(z[:, i]) * 1e2, + zmin=0, + zmax=1, + ), + ], + name="frame_$i", + traces=[0], + ) for i in 1:length(t) + ] + end + + sliders_attr = [ + attr(; + visible=length(t) > 1, + pad=attr(; b=10, t=60), + len=0.85, + x=0.15, + y=0.1, + t=0, + steps=[ + attr(; + label=round(t0; digits=2), + method="animate", + args=[ + ["frame_$i"], + attr(; + visible=[fill(false, i - 1); true; fill(false, length(t) - i)], + mode="immediate", + transition=attr(; duration=0), + frame=attr(; duration=5, redraw=true), + ), + ], + ) for (i, t0) in enumerate(t) + ], + currentvalue=attr(; + prefix="t = ", + suffix=" s", + xanchor="right", + font=attr(; color="#888", size=20), + ), + ), + ] + + buttons_attr = [ + attr(; + label="▶", # play symbol + method="animate", + args=[ + nothing, + attr(; + fromcurrent=true, + transition=(duration=dt_frame,), + frame=attr(; duration=dt_frame, redraw=true), + ), + ], + ), + attr(; + label="◼", # pause symbol + method="animate", + args=[ + [nothing], + attr(; + mode="immediate", + fromcurrent=true, + transition=attr(; duration=dt_frame), + frame=attr(; duration=dt_frame, redraw=true), + ), + ], + ), + ] + + #Layout + bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) + l = Layout(; + title=ph.name * ": " * string(key), + xaxis_title="x", + yaxis_title="y", + xaxis_range=[x0, xf], + yaxis_range=[x0, xf], + xaxis_ticksuffix=" cm", + yaxis_ticksuffix=" cm", + plot_bgcolor=plot_bgcolor, + paper_bgcolor=bgcolor, + xaxis_gridcolor=grid_color, + yaxis_gridcolor=grid_color, + xaxis_zerolinecolor=grid_color, + yaxis_zerolinecolor=grid_color, + font_color=text_color, + updatemenus=[ + attr(; + visible=length(t) > 1, + x=0.1, + y=0.05, + yanchor="top", + xanchor="center", + showactive=true, + direction="left", + type="buttons", + pad=attr(; r=10, t=70), + buttons=buttons_attr, + ), + ], + sliders=sliders_attr, + scene=attr(; + xaxis=attr(; + title="x", + range=[x0, xf], + ticksuffix=" cm", + backgroundcolor=plot_bgcolor, + gridcolor=grid_color, + zerolinecolor=grid_color, + ), + yaxis=attr(; + title="y", + range=[x0, xf], + ticksuffix=" cm", + backgroundcolor=plot_bgcolor, + gridcolor=grid_color, + zerolinecolor=grid_color, + ), + zaxis=attr(; + title="z", + range=[x0, xf], + ticksuffix=" cm", + backgroundcolor=plot_bgcolor, + gridcolor=grid_color, + zerolinecolor=grid_color, + ), + aspectmode="manual", + aspectratio=attr(; x=1, y=1, z=1), + ), + margin=attr(; t=50, l=0, r=0, b=0), + modebar=attr(; + orientation="h", bgcolor=bgcolor, color=text_color, activecolor=plot_bgcolor + ), + xaxis=attr(; constrain="domain"), + yaxis=attr(; scaleanchor="x"), + hovermode="closest", + ) + if height !== nothing l.height = height end if width !== nothing l.width = width end - if view_2d - h = scatter( x=(ph.x .+ ph.ux(ph.x,ph.y,ph.z,t0*1e-3))*1e2, - y=(ph.y .+ ph.uy(ph.x,ph.y,ph.z,t0*1e-3))*1e2, - mode="markers", - marker=attr(color=getproperty(ph,key)*factor, - showscale=colorbar, - colorscale=colormap, - colorbar=attr(ticksuffix=unit, title=string(key)), - cmin=cmin_key, - cmax=cmax_key, - size=4 - ), - text=round.(getproperty(ph,key)*factor,digits=4), - hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
$(string(key)): %{text}$unit") - else - h = scatter3d( x=(ph.x .+ ph.ux(ph.x,ph.y,ph.z,t0*1e-3))*1e2, - y=(ph.y .+ ph.uy(ph.x,ph.y,ph.z,t0*1e-3))*1e2, - z=(ph.z .+ ph.uz(ph.x,ph.y,ph.z,t0*1e-3))*1e2, - mode="markers", - marker=attr(color=getproperty(ph,key)*factor, - showscale=colorbar, - colorscale=colormap, - colorbar=attr(ticksuffix=unit, title=string(key)), - cmin=cmin_key, - cmax=cmax_key, - size=2 - ), - text=round.(getproperty(ph,key)*factor,digits=4), - hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
z: %{z:.1f} cm
$(string(key)): %{text}$unit") - end - config = PlotConfig( - displaylogo=false, - toImageButtonOptions=attr( - format="svg", # one of png, svg, jpeg, webp - ).fields, - modeBarButtonsToRemove=["zoom", "pan", "tableRotation", "resetCameraLastSave3d", "orbitRotation", "resetCameraDefault3d"] - ) - return plot_koma(h, l; config) + + return Plot(trace, l, frames) end """ @@ -950,96 +1382,114 @@ julia> plot_signal(raw) ``` """ function plot_signal( - raw::RawAcquisitionData; - width=nothing, - height=nothing, - slider=true, - show_sim_blocks=false, - darkmode=false, - range=[], - gl=false - ) - not_Koma = raw.params["systemVendor"] != "KomaMRI.jl" - t = [] - signal = [] - current_t0 = 0 - for p in raw.profiles - dt = p.head.sample_time_us != 0 ? p.head.sample_time_us * 1e-3 : 1 - t0 = p.head.acquisition_time_stamp * 1e-3 #This parameter is used in Koma to store the time offset - N = p.head.number_of_samples != 0 ? p.head.number_of_samples : 1 - if not_Koma - t0 = current_t0 * dt - current_t0 += N - end - if N != 1 - append!(t, t0.+(0:dt:dt*(N-1))) - else - append!(t, t0) - end - append!(signal, p.data[:,1]) #Just one coil - #To generate gap - append!(t, t[end]) - append!(signal, [Inf + Inf*1im]) - end - #Show simulation blocks - shapes = [] - annotations = [] - type_names = ["precession", "excitation"] - if !not_Koma && show_sim_blocks - t_sim_parts = raw.params["userParameters"]["t_sim_parts"] - type_sim_parts = raw.params["userParameters"]["type_sim_parts"] - - current_type = -1 - for i = eachindex(t_sim_parts[1:end-1]) - aux = rect( - xref="x", yref="paper", - x0=t_sim_parts[i]*1e3, y0=0, - x1=t_sim_parts[i+1]*1e3, y1=1, - fillcolor=type_sim_parts[i] ? "Purple" : "Blue", - opacity=.1, - layer="below", line_width=2, - ) - push!(shapes, aux) - - if type_sim_parts[i] != current_type - aux = attr(xref="x", yref="paper", x=t_sim_parts[i]*1e3, y=1, showarrow=false, text=type_names[type_sim_parts[i]+1]) - push!(annotations, aux) - current_type = type_sim_parts[i] - end - end - end - #PLOT - bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) - l = Layout(; hovermode="closest", - xaxis_title="", - modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor), - legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0), - plot_bgcolor=plot_bgcolor, - paper_bgcolor=bgcolor, - xaxis_gridcolor=grid_color, - yaxis_gridcolor=grid_color, - xaxis_zerolinecolor=grid_color, - yaxis_zerolinecolor=grid_color, - font_color=text_color, - yaxis_fixedrange = false, - xaxis=attr( - ticksuffix=" ms", - range=range[:], - rangeslider=attr(visible=slider), - rangeselector=attr( - buttons=[ - attr(count=1, - label="1m", - step=10, - stepmode="backward"), - attr(step="all") - ] - ), - ), - shapes = shapes, - annotations = annotations, - margin=attr(t=0,l=0,r=0,b=0) - ) + raw::RawAcquisitionData; + width=nothing, + height=nothing, + slider=true, + show_sim_blocks=false, + darkmode=false, + range=[], + gl=false, +) + not_Koma = raw.params["systemVendor"] != "KomaMRI.jl" + t = [] + signal = [] + current_t0 = 0 + for p in raw.profiles + dt = p.head.sample_time_us != 0 ? p.head.sample_time_us * 1e-3 : 1 + t0 = p.head.acquisition_time_stamp * 1e-3 #This parameter is used in Koma to store the time offset + N = p.head.number_of_samples != 0 ? p.head.number_of_samples : 1 + if not_Koma + t0 = current_t0 * dt + current_t0 += N + end + if N != 1 + append!(t, t0 .+ (0:dt:(dt * (N - 1)))) + else + append!(t, t0) + end + append!(signal, p.data[:, 1]) #Just one coil + #To generate gap + append!(t, t[end]) + append!(signal, [Inf + Inf * 1im]) + end + #Show simulation blocks + shapes = [] + annotations = [] + type_names = ["precession", "excitation"] + if !not_Koma && show_sim_blocks + t_sim_parts = raw.params["userParameters"]["t_sim_parts"] + type_sim_parts = raw.params["userParameters"]["type_sim_parts"] + + current_type = -1 + for i in eachindex(t_sim_parts[1:(end - 1)]) + aux = rect(; + xref="x", + yref="paper", + x0=t_sim_parts[i] * 1e3, + y0=0, + x1=t_sim_parts[i + 1] * 1e3, + y1=1, + fillcolor=type_sim_parts[i] ? "Purple" : "Blue", + opacity=0.1, + layer="below", + line_width=2, + ) + push!(shapes, aux) + + if type_sim_parts[i] != current_type + aux = attr(; + xref="x", + yref="paper", + x=t_sim_parts[i] * 1e3, + y=1, + showarrow=false, + text=type_names[type_sim_parts[i] + 1], + ) + push!(annotations, aux) + current_type = type_sim_parts[i] + end + end + end + #PLOT + bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode) + l = Layout(; + hovermode="closest", + xaxis_title="", + modebar=attr(; + orientation="h", + yanchor="bottom", + xanchor="right", + y=1, + x=0, + bgcolor=bgcolor, + color=text_color, + activecolor=plot_bgcolor, + ), + legend=attr(; orientation="h", yanchor="bottom", xanchor="left", y=1, x=0), + plot_bgcolor=plot_bgcolor, + paper_bgcolor=bgcolor, + xaxis_gridcolor=grid_color, + yaxis_gridcolor=grid_color, + xaxis_zerolinecolor=grid_color, + yaxis_zerolinecolor=grid_color, + font_color=text_color, + yaxis_fixedrange=false, + xaxis=attr(; + ticksuffix=" ms", + range=range[:], + rangeslider=attr(; visible=slider), + rangeselector=attr(; + buttons=[ + attr(; count=1, label="1m", step=10, stepmode="backward"), + attr(; step="all"), + ], + ), + ), + shapes=shapes, + annotations=annotations, + margin=attr(; t=0, l=0, r=0, b=0), + ) if height !== nothing l.height = height end @@ -1047,19 +1497,34 @@ function plot_signal( l.width = width end scatter_fun = gl ? scattergl : scatter - p = [scatter_fun() for j=1:3] - p[1] = scatter(x=t,y=abs.(signal), name="|S(t)|",hovertemplate="(%{x:.4f} ms, %{y:.3f} a.u.)") - p[2] = scatter_fun(x=t,y=real.(signal),name="Re{S(t)}",hovertemplate="(%{x:.4f} ms, %{y:.3f} a.u.)") - p[3] = scatter_fun(x=t,y=imag.(signal),name="Im{S(t)}",hovertemplate="(%{x:.4f} ms, %{y:.3f} a.u.)") - config = PlotConfig( - displaylogo=false, - toImageButtonOptions=attr( - format="svg", # one of png, svg, jpeg, webp - ).fields, - modeBarButtonsToRemove=["zoom", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"] - # modeBarButtonsToRemove=["zoom", "select2d", "lasso2d", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"] - ) - return plot_koma(p, l; config) + p = [scatter_fun() for j in 1:3] + p[1] = scatter(; + x=t, y=abs.(signal), name="|S(t)|", hovertemplate="(%{x:.4f} ms, %{y:.3f} a.u.)" + ) + p[2] = scatter_fun(; + x=t, y=real.(signal), name="Re{S(t)}", hovertemplate="(%{x:.4f} ms, %{y:.3f} a.u.)" + ) + p[3] = scatter_fun(; + x=t, y=imag.(signal), name="Im{S(t)}", hovertemplate="(%{x:.4f} ms, %{y:.3f} a.u.)" + ) + config = PlotConfig(; + displaylogo=false, + toImageButtonOptions=attr(; + format="svg", # one of png, svg, jpeg, webp + ).fields, + modeBarButtonsToRemove=[ + "zoom", + "autoScale", + "resetScale2d", + "pan", + "tableRotation", + "resetCameraLastSave", + "zoomIn", + "zoomOut", + ], + # modeBarButtonsToRemove=["zoom", "select2d", "lasso2d", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"] + ) + return plot_koma(p, l; config) end """ @@ -1074,29 +1539,29 @@ Generates an HTML table based on the dictionary `dict`. - `str`: (`::String`) dictionary as an HTML table """ function plot_dict(dict::Dict) - html = """ - - - - - - - - - - """ - i = 1 - for (key,val) = dict - html *= """ - - - - - - """ - i += 1 - end - html *= "
#NameValue
$i$(string(key))$(string(val))
" + html = """ + + + + + + + + + + """ + i = 1 + for (key, val) in dict + html *= """ + + + + + + """ + i += 1 + end + return html *= "
#NameValue
$i$(string(key))$(string(val))
" end """ @@ -1124,13 +1589,56 @@ julia> plot_seqd(seq) ``` """ function plot_seqd(seq::Sequence; sampling_params=KomaMRIBase.default_sampling_params()) - seqd = KomaMRIBase.discretize(seq; sampling_params) - Gx = scattergl(x=seqd.t*1e3, y=seqd.Gx*1e3, name="Gx", mode="markers+lines", marker_symbol=:circle) - Gy = scattergl(x=seqd.t*1e3, y=seqd.Gy*1e3, name="Gy", mode="markers+lines", marker_symbol=:circle) - Gz = scattergl(x=seqd.t*1e3, y=seqd.Gz*1e3, name="Gz", mode="markers+lines", marker_symbol=:circle) - B1_abs = scattergl(x=seqd.t*1e3, y=abs.(seqd.B1*1e6), name="|B1|", mode="markers+lines", marker_symbol=:circle) - B1_angle = scattergl(x=seqd.t*1e3, y=angle.(seqd.B1), name="∠B1", mode="markers+lines", marker_symbol=:circle) - ADC = scattergl(x=seqd.t[seqd.ADC]*1e3, y=ones(sum(seqd.ADC)), name="ADC", mode="markers", marker_symbol=:x) - B1_Δf = scattergl(x=seqd.t*1e3, y=abs.(seqd.Δf*1e-3), name="B1_Δf", mode="markers+lines", marker_symbol=:circle, visible="legendonly") - plot_koma([Gx,Gy,Gz,B1_abs,B1_angle,ADC,B1_Δf]) + seqd = KomaMRIBase.discretize(seq; sampling_params) + Gx = scattergl(; + x=seqd.t * 1e3, + y=seqd.Gx * 1e3, + name="Gx", + mode="markers+lines", + marker_symbol=:circle, + ) + Gy = scattergl(; + x=seqd.t * 1e3, + y=seqd.Gy * 1e3, + name="Gy", + mode="markers+lines", + marker_symbol=:circle, + ) + Gz = scattergl(; + x=seqd.t * 1e3, + y=seqd.Gz * 1e3, + name="Gz", + mode="markers+lines", + marker_symbol=:circle, + ) + B1_abs = scattergl(; + x=seqd.t * 1e3, + y=abs.(seqd.B1 * 1e6), + name="|B1|", + mode="markers+lines", + marker_symbol=:circle, + ) + B1_angle = scattergl(; + x=seqd.t * 1e3, + y=angle.(seqd.B1), + name="∠B1", + mode="markers+lines", + marker_symbol=:circle, + ) + ADC = scattergl(; + x=seqd.t[seqd.ADC] * 1e3, + y=zeros(sum(seqd.ADC)), + name="ADC", + mode="markers", + marker_symbol=:x, + ) + B1_Δf = scattergl(; + x=seqd.t * 1e3, + y=abs.(seqd.Δf * 1e-3), + name="B1_Δf", + mode="markers+lines", + marker_symbol=:circle, + visible="legendonly", + ) + return plot_koma([Gx, Gy, Gz, B1_abs, B1_angle, ADC, B1_Δf]) end diff --git a/docs/src/reference/2-koma-base.md b/docs/src/reference/2-koma-base.md index fa27fa4c4..9778de4c9 100644 --- a/docs/src/reference/2-koma-base.md +++ b/docs/src/reference/2-koma-base.md @@ -16,6 +16,30 @@ Phantom brain_phantom2D brain_phantom3D pelvis_phantom2D +heart_phantom +``` + +### `MotionModel` +```@docs +get_spin_coords +``` + +### `SimpleMotion <: MotionModel` + +```@docs +SimpleMotion +Translation +Rotation +HeartBeat +PeriodicTranslation +PeriodicRotation +PeriodicHeartBeat +``` + +### `ArbitraryMotion <: MotionModel` + +```@docs +ArbitraryMotion ``` ## `Sequence`-related functions diff --git a/docs/src/reference/4-koma-files.md b/docs/src/reference/4-koma-files.md index eec615a82..2640f7205 100644 --- a/docs/src/reference/4-koma-files.md +++ b/docs/src/reference/4-koma-files.md @@ -15,4 +15,6 @@ read_seq ```@docs read_phantom_jemris read_phantom_MRiLab +read_phantom +write_phantom ``` \ No newline at end of file diff --git a/examples/2.phantoms/brain_motion.phantom b/examples/2.phantoms/brain_motion.phantom index 99adac4bb..c3412d7f0 100644 Binary files a/examples/2.phantoms/brain_motion.phantom and b/examples/2.phantoms/brain_motion.phantom differ diff --git a/examples/2.phantoms/contracting_ring.phantom b/examples/2.phantoms/contracting_ring.phantom new file mode 100644 index 000000000..e5c7fb7c1 Binary files /dev/null and b/examples/2.phantoms/contracting_ring.phantom differ diff --git a/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl b/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl index 7e3c9933c..e2f63f0d9 100644 --- a/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl +++ b/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl @@ -1,7 +1,7 @@ using KomaMRI, Suppressor, MAT obj = @suppress read_phantom_jemris("../../../2.phantoms/brain.h5") -obj.uy = (x,y,z,t)-> 0.1f0 * t +obj.uy = (x,y,z,t)-> 0.1f0 * t # Hacer que el fichero brain_motion.phantom tenga este movimiento (con SimpleMotion y ArbitraryMotion) seq = @suppress read_seq("../sequences/EPI/epi_100x100_TE100_FOV230.seq") sys = Scanner() #Time diff --git a/src/ui/ExportUIFunctions.jl b/src/ui/ExportUIFunctions.jl index 22aa68fb1..665b38609 100644 --- a/src/ui/ExportUIFunctions.jl +++ b/src/ui/ExportUIFunctions.jl @@ -5,47 +5,50 @@ function setup_blink_window(; darkmode=true, frame=true, dev_tools=false, show_w komamri_src_ui = @__DIR__ komamri_root = dirname(dirname(komamri_src_ui)) # Asset folders - assets = AssetRegistry.register(komamri_root*"/assets/") - scripts = AssetRegistry.register(komamri_src_ui*"/scripts/") - css = AssetRegistry.register(komamri_src_ui*"/css/") + assets = AssetRegistry.register(komamri_root * "/assets/") + scripts = AssetRegistry.register(komamri_src_ui * "/scripts/") + css = AssetRegistry.register(komamri_src_ui * "/css/") # Images logo = joinpath(assets, "logo-dark.svg") home_image = joinpath(assets, "home-image.svg") - app_icon = komamri_root*"/assets/app-icon.png" + app_icon = komamri_root * "/assets/app-icon.png" # JS bsjs = joinpath(scripts, "bootstrap.bundle.min.js") #this already has Popper - bscss = joinpath(css,"bootstrap.min.css") - bsiconcss = joinpath(css,"bootstrap-icons.css") - jquery = joinpath(scripts,"jquery-3.4.1.slim.min.js") + bscss = joinpath(css, "bootstrap.min.css") + bsiconcss = joinpath(css, "bootstrap-icons.css") + jquery = joinpath(scripts, "jquery-3.4.1.slim.min.js") # mathjaxsetup = joinpath(scripts, "mathjaxsetup.js") # KaTeX katexrender = joinpath(scripts, "auto-render.min.js") - katexjs = joinpath(scripts,"katex.min.js") - katexcss = joinpath(css,"katex.min.css") + katexjs = joinpath(scripts, "katex.min.js") + katexcss = joinpath(css, "katex.min.css") # User defined JS and CSS - customcss = joinpath(css,"custom.css") - customjs = joinpath(scripts,"custom.js") - sidebarcss = joinpath(css,"sidebars.css") + customcss = joinpath(css, "custom.css") + customjs = joinpath(scripts, "custom.js") + sidebarcss = joinpath(css, "sidebars.css") # Custom icons - icons = joinpath(css,"icons.css") + icons = joinpath(css, "icons.css") ## WINDOW - w = Blink.Window(Dict( - "title"=>"KomaUI", - "autoHideMenuBar"=>true, - "frame"=>frame, #removes title bar - "node-integration"=>true, - :icon=>app_icon, - "width"=>1200, - "height"=>800, - :show=>show_window - ), async=false); + w = Blink.Window( + Dict( + "title" => "KomaUI", + "autoHideMenuBar" => true, + "frame" => frame, #removes title bar + "node-integration" => true, + :icon => app_icon, + "width" => 1200, + "height" => 800, + :show => show_window, + ); + async=false, + ) ## NAV BAR - sidebar = open(f->read(f, String), komamri_src_ui*"/html/sidebar.html") - sidebar = replace(sidebar, "LOGO"=>logo) + sidebar = open(f -> read(f, String), komamri_src_ui * "/html/sidebar.html") + sidebar = replace(sidebar, "LOGO" => logo) ## CONTENT - index = open(f->read(f, String), komamri_src_ui*"/html/index.html") - index = replace(index, "ICON"=>home_image) + index = open(f -> read(f, String), komamri_src_ui * "/html/index.html") + index = replace(index, "ICON" => home_image) @sync begin ## CSS @@ -66,9 +69,9 @@ function setup_blink_window(; darkmode=true, frame=true, dev_tools=false, show_w end #Update GUI's home - body!(w, *(sidebar, index), async=false) #async false is important to set the background correctly + body!(w, *(sidebar, index); async=false) #async false is important to set the background correctly if darkmode - @js_ w document.getElementById("main").style="background-color:rgb(13,16,17);" + @js_ w document.getElementById("main").style = "background-color:rgb(13,16,17);" end # Return the Blink window if dev_tools @@ -129,30 +132,32 @@ function setup_raw() "encodedSize" => [2, 2, 1], "reconSize" => [2, 2, 1], "number_of_samples" => 4, - "encodedFOV" => [100., 100., 1], - "trajectory" => "other" + "encodedFOV" => [100.0, 100.0, 1], + "trajectory" => "other", ), [ - KomaMRICore.Profile(AcquisitionHeader( - trajectory_dimensions=2, sample_time_us=1), - [0. 0. 1 1; 0 1 1 1]./2, - reshape([0.; 0im; 0; 0], 4, 1) - ) - ] + KomaMRICore.Profile( + AcquisitionHeader(; trajectory_dimensions=2, sample_time_us=1), + [0.0 0.0 1 1; 0 1 1 1] ./ 2, + reshape([0.0; 0im; 0; 0], 4, 1), + ), + ], ) # Return the default RawAcquisitionData struct return raw end - """ Updates the blink window with a loading content """ function display_loading!(w::Window, msg::String) komamri_src_ui = @__DIR__ - loading = replace(open((f) -> read(f, String), komamri_src_ui*"/html/loading.html"), "LOADDES"=>msg) - content!(w, "div#content", loading) + loading = replace( + open((f) -> read(f, String), komamri_src_ui * "/html/loading.html"), + "LOADDES" => msg, + ) + return content!(w, "div#content", loading) end """ @@ -161,25 +166,32 @@ Fileficker callback when sequence is loaded function callback_filepicker(filename::String, w::Window, seq::Sequence) if filename != "" filename_extension = splitext(filename)[end] - if filename_extension ==".seqk" # Koma - seq = JLD2.load(FileIO.File{FileIO.DataFormat{:JLD2}}(filename),"seq") + if filename_extension == ".seqk" # Koma + seq = JLD2.load(FileIO.File{FileIO.DataFormat{:JLD2}}(filename), "seq") elseif filename_extension == ".seq" # Pulseq seq = read_seq(filename) # Pulseq read end - @js_ w (@var name = $(basename(filename)); - document.getElementById("seqname").innerHTML = "" + name + ""; - Toasty("0", "Loaded " + name + " successfully", """ - - """)) + @js_ w ( + @var name = $(basename(filename)); + document.getElementById("seqname").innerHTML = + "" + name + ""; + Toasty( + "0", + "Loaded " + name + " successfully", + """ + +""", + ) + ) return seq end return seq @@ -192,24 +204,31 @@ function callback_filepicker(filename::String, w::Window, obj::Phantom) if filename != "" filename_extension = splitext(filename)[end] if filename_extension == ".phantom" # Koma - obj = JLD2.load(FileIO.File{FileIO.DataFormat{:JLD2}}(filename),"phantom") + obj = read_phantom(filename) elseif filename_extension == ".h5" # JEMRIS obj = read_phantom_jemris(filename) end - @js_ w (@var name = $(basename(filename)); - document.getElementById("phaname").innerHTML="" + name + "";; - Toasty("0", "Loaded " + name + " successfully", """ - - """)) + @js_ w ( + @var name = $(basename(filename)); + document.getElementById("phaname").innerHTML = + "" + name + ""; + Toasty( + "0", + "Loaded " + name + " successfully", + """ + +""", + ) + ) return obj end return obj @@ -224,20 +243,27 @@ function callback_filepicker(filename::String, w::Window, raw::RawAcquisitionDat if raw.params["systemVendor"] != "KomaMRI.jl" @warn "ISMRMRD files generated externally could cause problems during the reconstruction. We are currently improving compatibility." end - @js_ w (@var name = $(basename(filename)); - document.getElementById("rawname").innerHTML="" + name + "";; - Toasty("0", "Loaded " + name + " successfully", """ - - """)) + @js_ w ( + @var name = $(basename(filename)); + document.getElementById("rawname").innerHTML = + "" + name + ""; + Toasty( + "0", + "Loaded " + name + " successfully", + """ + +""", + ) + ) return raw end return raw @@ -281,20 +307,40 @@ end """ Updates the UI with phantom plots """ -function view_ui_phantom!(obj::Phantom, w::Window, seq::Sequence, buttons_obj::Vector{Widget{:button, Int64}}; key=:ρ, darkmode=true) +function view_ui_phantom!( + obj::Phantom, + w::Window, + seq::Sequence, + buttons_obj::Vector{Widget{:button,Int64}}; + key=:ρ, + darkmode=true, +) display_loading!(w, "Plotting phantom ...") - widget_plot = @manipulate for t0_ms in 1e3*range(0, dur(seq); length=5) - plot_phantom_map(obj, key; t0=t0_ms, darkmode) - end + widget_plot = plot_phantom_map(obj, key; intermediate_time_samples=5, darkmode) div_content = dom"div"(hbox(buttons_obj...), widget_plot) content!(w, "div#content", div_content) @js_ w document.getElementById("content").dataset.content = "phantom" end -function view_ui!(obj::Phantom, w::Window, seq::Sequence, buttons_obj::Vector{Widget{:button, Int64}}; key=:ρ, darkmode=true) - view_ui_phantom!(obj, w, seq, buttons_obj; key, darkmode) +function view_ui!( + obj::Phantom, + w::Window, + seq::Sequence, + buttons_obj::Vector{Widget{:button,Int64}}; + key=:ρ, + darkmode=true, +) + return view_ui_phantom!(obj, w, seq, buttons_obj; key, darkmode) end -function view_ui!(cnt::Integer, w::Window, obj::Phantom, seq::Sequence, buttons_obj::Vector{Widget{:button, Int64}}; key=:ρ, darkmode=true) - view_ui_phantom!(obj, w, seq, buttons_obj; key, darkmode) +function view_ui!( + cnt::Integer, + w::Window, + obj::Phantom, + seq::Sequence, + buttons_obj::Vector{Widget{:button,Int64}}; + key=:ρ, + darkmode=true, +) + return view_ui_phantom!(obj, w, seq, buttons_obj; key, darkmode) end """ @@ -302,35 +348,36 @@ Updates the UI with scanner information """ function view_ui!(sys::Scanner, w::Window) display_loading!(w, "Displaying scanner parameters ...") - sys_dict = Dict("B0" => sys.B0, - "B1" => sys.B1, - "Gmax" => sys.Gmax, - "Smax" => sys.Smax, - "ADC_dt" => sys.ADC_Δt, - "seq_dt" => sys.seq_Δt, - "GR_dt" => sys.GR_Δt, - "RF_dt" => sys.RF_Δt, - "RF_ring_down_T" => sys.RF_ring_down_T, - "RF_dead_time_T" => sys.RF_dead_time_T, - "ADC_dead_time_T" => sys.ADC_dead_time_T) + sys_dict = Dict( + "B0" => sys.B0, + "B1" => sys.B1, + "Gmax" => sys.Gmax, + "Smax" => sys.Smax, + "ADC_dt" => sys.ADC_Δt, + "seq_dt" => sys.seq_Δt, + "GR_dt" => sys.GR_Δt, + "RF_dt" => sys.RF_Δt, + "RF_ring_down_T" => sys.RF_ring_down_T, + "RF_dead_time_T" => sys.RF_dead_time_T, + "ADC_dead_time_T" => sys.ADC_dead_time_T, + ) plt = plot_dict(sys_dict) title = """

Scanner parameters

""" - content!(w, "div#content", title*plt) + content!(w, "div#content", title * plt) @js_ w document.getElementById("content").dataset.content = "scanneparams" end """ Updates the UI with simulation parameters information """ -function view_ui!(sim_params::Dict{String, Any}, w::Window) +function view_ui!(sim_params::Dict{String,Any}, w::Window) display_loading!(w, "Displaying simulation parameters ...") plt = plot_dict(sim_params) title = """

Simulation parameters

""" - content!(w, "div#content", title*plt) + content!(w, "div#content", title * plt) @js_ w document.getElementById("content").dataset.content = "simparams" end - """ Updates the UI with raw signal plot """ @@ -344,11 +391,11 @@ end """ Updates the UI with reconstruction parameters information """ -function view_ui!(rec_params::Dict{Symbol, Any}, w::Window) +function view_ui!(rec_params::Dict{Symbol,Any}, w::Window) display_loading!(w, "Displaying reconstruction parameters ...") plt = plot_dict(rec_params) title = """

Reconstruction parameters

""" - content!(w, "div#content", title*plt) + content!(w, "div#content", title * plt) @js_ w document.getElementById("content").dataset.content = "recparams" end @@ -361,7 +408,13 @@ function view_ui!(img::Array, w::Window; type="absi", darkmode=true) display_loading!(w, "Plotting image magnitude ...") widget_plot = @manipulate for slice in 1:size(img, 3) aux = abs.(img) * prod(size(img)[1:2]) - plot_image(aux[:, :, slice], zmin=minimum(aux[:]), zmax=maximum(aux[:]); darkmode, title="Reconstruction ($slice/$(size(img, 3)))") + plot_image( + aux[:, :, slice]; + zmin=minimum(aux[:]), + zmax=maximum(aux[:]), + darkmode, + title="Reconstruction ($slice/$(size(img, 3)))", + ) end content!(w, "div#content", dom"div"(widget_plot)) @js_ w document.getElementById("content").dataset.content = "absi" @@ -369,7 +422,13 @@ function view_ui!(img::Array, w::Window; type="absi", darkmode=true) display_loading!(w, "Plotting image phase ...") widget_plot = @manipulate for slice in 1:size(img, 3) aux = angle.(img[:, :, slice]) - plot_image(aux, zmin=-π, zmax=π; darkmode, title="Reconstruction ($slice/$(size(img, 3)))") + plot_image( + aux; + zmin=-π, + zmax=π, + darkmode, + title="Reconstruction ($slice/$(size(img, 3)))", + ) end content!(w, "div#content", dom"div"(widget_plot)) @js_ w document.getElementById("content").dataset.content = "angi" @@ -378,7 +437,13 @@ function view_ui!(img::Array, w::Window; type="absi", darkmode=true) widget_plot = @manipulate for slice in 1:size(img, 3) kspace = fftc(img) aux = log.(abs.(kspace[:, :, slice]) .+ 1) - plot_image(aux, zmin=0, zmax=.1*maximum(aux[:]); darkmode, title="Reconstruction ($slice/$(size(img, 3)))") + plot_image( + aux; + zmin=0, + zmax=0.1 * maximum(aux[:]), + darkmode, + title="Reconstruction ($slice/$(size(img, 3)))", + ) end content!(w, "div#content", dom"div"(widget_plot)) @js_ w document.getElementById("content").dataset.content = "absk" @@ -388,30 +453,42 @@ end """ Saves to .mat in the UI (displays a Toast message en the UI) """ -function save_ui!(w::Window, seq::Sequence, obj::Phantom, sys::Scanner, raw, img, rec_params, mat_folder; type="all") - if type=="all" - str_toast = export_2_mat(seq, obj, sys, raw, rec_params, img, mat_folder; type, matfilename="") - @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);) +function save_ui!( + w::Window, + seq::Sequence, + obj::Phantom, + sys::Scanner, + raw, + img, + rec_params, + mat_folder; + type="all", +) + if type == "all" + str_toast = export_2_mat( + seq, obj, sys, raw, rec_params, img, mat_folder; type, matfilename="" + ) + @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files", msg)) @js_ w document.getElementById("content").dataset.content = "matfolder" - elseif type=="sequence" + elseif type == "sequence" str_toast = export_2_mat(seq, obj, sys, raw, rec_params, img, mat_folder; type) - @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);) + @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files", msg)) @js_ w document.getElementById("content").dataset.content = "matfolderseq" - elseif type=="phantom" + elseif type == "phantom" str_toast = export_2_mat(seq, obj, sys, raw, rec_params, img, mat_folder; type) - @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);) + @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files", msg)) @js_ w document.getElementById("content").dataset.content = "matfolderpha" - elseif type=="scanner" + elseif type == "scanner" str_toast = export_2_mat(seq, obj, sys, raw, rec_params, img, mat_folder; type) - @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);) + @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files", msg)) @js_ w document.getElementById("content").dataset.content = "matfoldersca" - elseif type=="raw" + elseif type == "raw" str_toast = export_2_mat(seq, obj, sys, raw, rec_params, img, mat_folder; type) - @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);) + @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files", msg)) @js_ w document.getElementById("content").dataset.content = "matfolderraw" - elseif type=="image" + elseif type == "image" str_toast = export_2_mat(seq, obj, sys, raw, rec_params, img, mat_folder; type) - @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files" , msg);) + @js_ w (@var msg = $str_toast; Toasty("1", "Saved .mat files", msg)) @js_ w document.getElementById("content").dataset.content = "matfolderima" - end + end end