From afa7e9fe9897cb83c0a765d0515d85131cd8aed7 Mon Sep 17 00:00:00 2001 From: Pablo Villacorta Aylagas Date: Wed, 27 May 2026 10:21:20 +0200 Subject: [PATCH 1/2] (WIP) First commit. Tests with variable rho (Pending tests for KomaFiles). --- KomaMRIBase/src/KomaMRIBase.jl | 8 +- KomaMRIBase/src/datatypes/Phantom.jl | 23 ++-- .../phantom/TimeDependentProperty.jl | 39 +++++++ KomaMRIBase/src/motion/MotionList.jl | 2 - KomaMRIBase/src/motion/SpinSpan.jl | 15 +++ .../src/{motion => timing}/Interpolation.jl | 0 .../src/{motion => timing}/TimeCurve.jl | 0 KomaMRIBase/test/runtests.jl | 29 +++++ KomaMRICore/src/simulation/Flow.jl | 58 +++++----- KomaMRICore/src/simulation/Functors.jl | 1 + .../SimMethods/Bloch/cpu/BlochCPU.jl | 11 +- .../SimMethods/Bloch/gpu/BlochGPU.jl | 10 +- .../SimMethods/BlochDict/BlochDict.jl | 10 +- .../BlochMagnus/cpu/BlochMagnusCPU.jl | 6 +- .../SimMethods/BlochSimple/BlochSimple.jl | 12 +- .../simulation/SimMethods/SimulationMethod.jl | 2 +- KomaMRIPlots/src/ui/DisplayFunctions.jl | 109 ++++++++++++------ 17 files changed, 236 insertions(+), 99 deletions(-) create mode 100644 KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl rename KomaMRIBase/src/{motion => timing}/Interpolation.jl (100%) rename KomaMRIBase/src/{motion => timing}/TimeCurve.jl (100%) diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl index 0f3624e11..a13dcaf35 100644 --- a/KomaMRIBase/src/KomaMRIBase.jl +++ b/KomaMRIBase/src/KomaMRIBase.jl @@ -70,6 +70,9 @@ include("datatypes/sequence/TimingChecks.jl") include("datatypes/sequence/HardwareChecks.jl") include("datatypes/sequence/AddBlockMacro.jl") include("datatypes/sequence/DelayDuration.jl") +# Timing +include("timing/TimeCurve.jl") +include("timing/Interpolation.jl") # Motion include("motion/MotionList.jl") include("motion/NoMotion.jl") @@ -102,8 +105,9 @@ export Translate, TranslateX, TranslateY, TranslateZ export Rotate, RotateX, RotateY, RotateZ, CenterOfMass export HeartBeat, Path, FlowPath export TimeRange, Periodic, TimeCurve -export SpinRange, AllSpins -export get_spin_coords +export SpinRange, AllSpins, spin_indicator +export get_spin_coords, get_spin_property, get_phantom_property, get_spin_property_at_end +export TimeDependentProperty # Secondary export get_kspace, rotx, roty, rotz # Additionals diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl index b09e97b7c..42758e1b2 100644 --- a/KomaMRIBase/src/datatypes/Phantom.jl +++ b/KomaMRIBase/src/datatypes/Phantom.jl @@ -1,3 +1,5 @@ +include("phantom/TimeDependentProperty.jl") + """ obj = Phantom(name, x, y, z, ρ, T1, T2, T2s, Δw, Dλ1, Dλ2, Dθ, motion) @@ -9,7 +11,7 @@ a property value representing a spin. This struct serves as an input for the sim - `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 -- `ρ`: (`::AbstractVector{T<:Real}`) spin proton density vector +- `ρ`: (`::Union{AbstractVector{T}, TimeDependentProperty{T}}`) spin proton density - `T1`: (`::AbstractVector{T<:Real}`, `[s]`) spin T1 parameter vector - `T2`: (`::AbstractVector{T<:Real}`, `[s]`) spin T2 parameter vector - `T2s`: (`::AbstractVector{T<:Real}`, `[s]`) spin T2s parameter vector @@ -34,7 +36,7 @@ julia> obj.ρ x::AbstractVector{T} = @isdefined(T) ? T[] : Float64[] y::AbstractVector{T} = zeros(eltype(x), size(x)) z::AbstractVector{T} = zeros(eltype(x), size(x)) - ρ::AbstractVector{T} = ones(eltype(x), size(x)) + ρ::Union{AbstractVector{T}, TimeDependentProperty{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 @@ -53,9 +55,12 @@ end const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String, fieldnames(Phantom)) const VECTOR_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) <: AbstractVector, fieldnames(Phantom)) +get_phantom_property(obj::Phantom, key::Symbol) = get_property_value(getfield(obj, key)) +get_property_value(p) = p + """Size and length of a phantom""" -size(x::Phantom) = size(x.ρ) -Base.length(x::Phantom) = length(x.ρ) +size(x::Phantom) = size(x.x) +Base.length(x::Phantom) = length(x.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 @@ -99,10 +104,12 @@ end for field in VECTOR_PHANTOM_FIELDS push!(fields, (field, [getfield(obj1, field); getfield(obj2, field)])) end - return Phantom(; - name = name, - fields..., - motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2))) + return Phantom(; + name=name, + fields..., + ρ=vcat_property(getfield(obj1, :ρ), getfield(obj2, :ρ)), + motion=vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)), + ) end """Scalar multiplication of a phantom""" diff --git a/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl b/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl new file mode 100644 index 000000000..704129b25 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl @@ -0,0 +1,39 @@ +struct TimeDependentProperty{T<:Real} + value::AbstractMatrix{T} + time::TimeCurve{T} +end + +get_property_value(p::TimeDependentProperty) = p.value[:, 1] + +get_spin_property(p::AbstractVector, t) = p +function get_spin_property(p::TimeDependentProperty, t) + Ns = size(p.value, 1) + t_unit = unit_time(t, p.time) + return resample(interpolate(p.value, Gridded(Linear()), Val(Ns), t_unit), t_unit) +end + +get_spin_property_at_end(p::AbstractMatrix) = selectdim(p, 2, size(p, 2)) +get_spin_property_at_end(p::AbstractVector) = p + +Base.:(==)(p1::TimeDependentProperty, p2::TimeDependentProperty) = p1.value == p2.value && p1.time == p2.time +Base.:(≈)(p1::TimeDependentProperty, p2::TimeDependentProperty) = p1.value ≈ p2.value && p1.time ≈ p2.time +Base.:(==)(::TimeDependentProperty, ::AbstractVector) = false +Base.:(==)(::AbstractVector, ::TimeDependentProperty) = false +Base.:(≈)(::TimeDependentProperty, ::AbstractVector) = false +Base.:(≈)(::AbstractVector, ::TimeDependentProperty) = false + +Base.getindex(p::TimeDependentProperty, i) = TimeDependentProperty(p.value[i, :], p.time) +Base.view(p::TimeDependentProperty, i) = TimeDependentProperty(@view(p.value[i, :]), p.time) + +times(::AbstractVector) = Float64[] +times(p::TimeDependentProperty) = times(p.time) + +*(α::Real, p::TimeDependentProperty) = TimeDependentProperty(α .* p.value, p.time) + +vcat_property(p1::AbstractVector, p2::AbstractVector) = [p1; p2] +vcat_property(p1::TimeDependentProperty, p2::AbstractVector) = TimeDependentProperty([p1.value; repeat(p2, 1, size(p1.value, 2))], p1.time) +vcat_property(p1::AbstractVector, p2::TimeDependentProperty) = TimeDependentProperty([repeat(p1, 1, size(p2.value, 2)); p2.value], p2.time) +vcat_property(p1::TimeDependentProperty, p2::TimeDependentProperty) = begin + @assert p1.time == p2.time "Time ranges must be the same" + return TimeDependentProperty([p1.value; p2.value], p1.time) +end \ No newline at end of file diff --git a/KomaMRIBase/src/motion/MotionList.jl b/KomaMRIBase/src/motion/MotionList.jl index 99d6ecbad..207a173a7 100644 --- a/KomaMRIBase/src/motion/MotionList.jl +++ b/KomaMRIBase/src/motion/MotionList.jl @@ -1,6 +1,4 @@ -include("Interpolation.jl") include("SpinSpan.jl") -include("TimeCurve.jl") include("Action.jl") include("Motion.jl") diff --git a/KomaMRIBase/src/motion/SpinSpan.jl b/KomaMRIBase/src/motion/SpinSpan.jl index fa6060aa3..cd0d3f003 100644 --- a/KomaMRIBase/src/motion/SpinSpan.jl +++ b/KomaMRIBase/src/motion/SpinSpan.jl @@ -72,3 +72,18 @@ expand(sr::SpinRange, Ns::Int) = sr intersect_idx(a, b) = findall(x -> x in a, b) intersect_idx(a, b::BitVector) = findall(x -> x in a, findall(x->x==true, b)) intersect_idx(a::BitVector, b) = findall(x -> x in findall(x->x==true, a), b) + +"""Per-spin weights for broadcast masking (GPU-friendly; no gather on `replace_by`).""" +function spin_indicator(spins::AllSpins, w::AbstractVector{T}) where {T} + fill!(w, one(T)) + return w +end + +function spin_indicator(spins::SpinRange, w::AbstractVector{T}) where {T} + fill!(w, zero(T)) + r = spins.range + @inbounds for i in r + i <= length(w) && (w[i] = one(T)) + end + return w +end diff --git a/KomaMRIBase/src/motion/Interpolation.jl b/KomaMRIBase/src/timing/Interpolation.jl similarity index 100% rename from KomaMRIBase/src/motion/Interpolation.jl rename to KomaMRIBase/src/timing/Interpolation.jl diff --git a/KomaMRIBase/src/motion/TimeCurve.jl b/KomaMRIBase/src/timing/TimeCurve.jl similarity index 100% rename from KomaMRIBase/src/motion/TimeCurve.jl rename to KomaMRIBase/src/timing/TimeCurve.jl diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl index 17ac2c5fe..1677233e8 100644 --- a/KomaMRIBase/test/runtests.jl +++ b/KomaMRIBase/test/runtests.jl @@ -1572,6 +1572,25 @@ end obj2.motion = tr oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba + + tc = TimeRange(0.0, 1.0) + pd = [1.0 0.0; 0.5 1.0] + ρt = TimeDependentProperty(pd, tc) + # TimeDependentProperty + AbstractVector + obj1.ρ = ρt + obj2.ρ = ρ + oba.ρ = KomaMRIBase.vcat_property(obj1.ρ, obj2.ρ) + @test obj1 + obj2 == oba + # AbstractVector + TimeDependentProperty + obj1.ρ = ρ + obj2.ρ = ρt + oba.ρ = KomaMRIBase.vcat_property(obj1.ρ, obj2.ρ) + @test obj1 + obj2 == oba + # TimeDependentProperty + TimeDependentProperty + obj1.ρ = ρt + obj2.ρ = ρt + oba.ρ = KomaMRIBase.vcat_property(obj1.ρ, obj2.ρ) + @test obj1 + obj2 == oba end @testset "Scalar multiplication" begin 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θ, motion=MotionList(tr, rt)) @@ -1579,6 +1598,16 @@ end 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θ, motion=MotionList(tr, rt)) @test c * obj1 == obc end + @testset "Time-dependent ρ" begin + tc = TimeRange(0.0, 1.0) + pd = [1.0 0.0; 0.5 1.0] + ρt = TimeDependentProperty(pd, tc) + obj = Phantom(ρ=ρt, x=x[1:2], y=y[1:2], z=z[1:2]) + @test KomaMRIBase.get_spin_property(obj.ρ, [0.5]) ≈ [0.5, 0.75] + ρt2 = KomaMRIBase.get_spin_property(obj.ρ, [0.0, 1.0]) + @test ρt2 ≈ [1.0 0.0; 0.5 1.0] + @test KomaMRIBase.get_spin_property_at_end(ρt2) ≈ [1.0, 1.0] + end @testset "Brain Phantom 2D" begin ph = brain_phantom2D() @test ph.name == "brain2D_axial" diff --git a/KomaMRICore/src/simulation/Flow.jl b/KomaMRICore/src/simulation/Flow.jl index 99a782ba8..495cb8dd5 100644 --- a/KomaMRICore/src/simulation/Flow.jl +++ b/KomaMRICore/src/simulation/Flow.jl @@ -3,17 +3,29 @@ function outflow_spin_reset!(args...; kwargs...) end function outflow_spin_reset!(spin_state_matrix, t, ml::MotionList; replace_by=0, seq_t=0, add_t0=false) - for m in ml.motions + for m in ml.motions outflow_spin_reset!(spin_state_matrix, t, m; replace_by=replace_by, seq_t=seq_t, add_t0=add_t0) end return nothing end -function outflow_spin_reset!(spin_state_matrix, t, m::Motion; replace_by=0, seq_t=0, add_t0=false) +function outflow_spin_reset!(spin_state_matrix, t, m::Motion; replace_by=0, seq_t=0, add_t0=false) outflow_spin_reset!(spin_state_matrix, t, m.action, m.time, m.spins; replace_by=replace_by, seq_t=seq_t, add_t0=add_t0) return nothing end +function _flow_spin_weights(spin_span, prototype::AbstractVector{T}) where {T} + w = similar(prototype) + return spin_indicator(spin_span, w) +end + +_flow_reset_mask(mask::AbstractVector) = @view mask[lastindex(mask):lastindex(mask)] +_flow_reset_mask(mask::AbstractMatrix) = selectdim(mask, 2, size(mask, 2)) + +flow_replace_mag(replace_by::AbstractVector) = replace_by +flow_replace_mag(replace_by::AbstractMatrix) = selectdim(replace_by, 2, size(replace_by, 2)) +flow_replace_mag(replace_by::Number) = replace_by + function outflow_spin_reset!( spin_state_matrix::AbstractArray, t, @@ -24,17 +36,14 @@ function outflow_spin_reset!( seq_t=0, add_t0=false, ) - # Initialize time: add t0 and normalize ts = KomaMRIBase.unit_time(init_time(t, seq_t, add_t0), time_curve) - # Get spin state range affected by the spin span - idx = KomaMRIBase.get_indexing_range(spin_span) - spin_state_matrix = @view(spin_state_matrix[idx, :]) - replace_by = replace_view(replace_by, idx) - # Obtain mask + prototype = ndims(spin_state_matrix) == 1 ? spin_state_matrix : selectdim(spin_state_matrix, 2, 1) + w = _flow_spin_weights(spin_span, prototype) mask = get_mask(action.spin_reset, ts) - # Modify spin state: reset and replace by initial value - spin_state_matrix .*= (1 .- mask) - spin_state_matrix .+= replace_by .* mask + indicator = reshape(w, :, 1) + weighted_mask = mask .* indicator + spin_state_matrix .*= (1 .- weighted_mask) + spin_state_matrix .+= replace_by .* weighted_mask return nothing end @@ -48,20 +57,14 @@ function outflow_spin_reset!( seq_t=0, add_t0=false, ) - # Initialize time: add t0 and normalize ts = KomaMRIBase.unit_time(init_time(t, seq_t, add_t0), time_curve) - # Get spin state range affected by the spin span - idx = KomaMRIBase.get_indexing_range(spin_span) - M = @view(M[idx]) - replace_by = replace_view(replace_by, idx) - # Obtain mask - mask = get_mask(action.spin_reset, ts) - mask = @view(mask[:, end]) - # Modify spin state: reset and replace by initial value - M.xy .*= (1 .- mask) - M.z .*= (1 .- mask) - M.xy .+= 0 .* mask - M.z .+= replace_by .* mask + w = _flow_spin_weights(spin_span, M.z) + mask = _flow_reset_mask(get_mask(action.spin_reset, ts)) + ρ = flow_replace_mag(replace_by) + M.xy .*= (1 .- mask .* w) + M.z .*= (1 .- mask .* w) + M.xy .+= zero(eltype(M.z)) .* mask .* w + M.z .+= ρ .* mask .* w return nothing end @@ -73,13 +76,6 @@ function init_time(t, seq_t, add_t0) return t end -function replace_view(replace_by::AbstractArray, idx) - return @view(replace_by[idx]) -end -function replace_view(replace_by, idx) - return replace_by -end - function get_mask(spin_reset, t) itp = KomaMRIBase.interpolate(spin_reset, KomaMRIBase.Gridded(KomaMRIBase.Constant{KomaMRIBase.Next}()), Val(size(spin_reset, 1)), t) return KomaMRIBase.resample(itp, t) diff --git a/KomaMRICore/src/simulation/Functors.jl b/KomaMRICore/src/simulation/Functors.jl index 52dcc4fb2..4aac9de7d 100644 --- a/KomaMRICore/src/simulation/Functors.jl +++ b/KomaMRICore/src/simulation/Functors.jl @@ -112,6 +112,7 @@ adapt_storage(T::Type{<:Real}, xs::MotionList) = MotionList(paramtype.(T, xs.mot @functor Path @functor FlowPath @functor TimeCurve +@functor TimeDependentProperty # Spinor @functor Spinor # DiscreteSequence diff --git a/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl b/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl index bf3f6004b..910c8f5ab 100644 --- a/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl +++ b/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl @@ -65,6 +65,7 @@ function run_spin_precession!( fill!(ϕ, zero(T)) block_time = zero(T) sample = 1 + ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[1]) @. Bz_old = x * seq.Gx[1] + y * seq.Gy[1] + z * seq.Gz[1] + ΔBz #Simulation @@ -91,9 +92,9 @@ function run_spin_precession!( end #Final Spin-State @. M.xy = M.xy * exp(-block_time / p.T2) * cis(ϕ) - @. M.z = M.z * exp(-block_time / p.T1) + p.ρ * (T(1) - exp(-block_time / p.T1)) + @. M.z = M.z * exp(-block_time / p.T1) + ρ_end * (T(1) - exp(-block_time / p.T1)) #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end) return nothing end @@ -133,6 +134,8 @@ function run_spin_excitation!( for i in eachindex(seq.Δt) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i]) + ρ = get_spin_property(p.ρ, seq.t[i, :]) + ρ_end = get_spin_property_at_end(ρ) #Effective field @. Bz = (seq.Gx[i] * x + seq.Gy[i] * y + seq.Gz[i] * z) + ΔBz - seq.Δf[i] / T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z) @. B = sqrt(abs(seq.B1[i])^2 + abs(Bz)^2) @@ -144,9 +147,9 @@ function run_spin_excitation!( mul!(Spinor(α, β), M, Maux_xy, Maux_z) #Relaxation @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2) - @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + p.ρ * (T(1) - exp(-seq.Δt[i] / p.T1)) + @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + ρ * (T(1) - exp(-seq.Δt[i] / p.T1)) #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=ρ_end) #Acquire signal if seq.ADC[i + 1] # ADC at the end of the time step sig[sample] = sum(M.xy) diff --git a/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl b/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl index dba40b09f..d6e7f5259 100644 --- a/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl +++ b/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl @@ -37,13 +37,14 @@ function run_spin_precession!( ) where {T<:Real, SM<:BlochLikeSimMethods} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') + ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) has_adc = length(sig) > 0 #Precession precession_kernel!(backend, groupsize)( pre.sig_output, M.xy, M.z, - x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)), + x, y, z, pre.ΔBz, p.T1, p.T2, ρ_end, UInt32(length(M.xy)), seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.ADC, UInt32(length(seq.t)), Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), sim_method, @@ -57,7 +58,7 @@ function run_spin_precession!( end #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end) return nothing end @@ -74,13 +75,14 @@ function run_spin_excitation!( ) where {T<:Real, SM<:BlochLikeSimMethods} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') + ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) has_adc = length(sig) > 0 #Excitation excitation_kernel!(backend, groupsize)( pre.sig_output, M.xy, M.z, - x, y, z, pre.ΔBz, p.T1, p.T2, p.ρ, UInt32(length(M.xy)), + x, y, z, pre.ΔBz, p.T1, p.T2, ρ_end, UInt32(length(M.xy)), seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.Δf, seq.B1, seq.ψ, seq.ADC, UInt32(length(seq.t)), Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), sim_method, @@ -94,7 +96,7 @@ function run_spin_excitation!( end #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t', p.motion; replace_by=p.ρ) # TODO: reset state inside kernel + outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end) # TODO: reset state inside kernel return nothing end diff --git a/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl b/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl index 7db8d6f0f..1209d31f5 100644 --- a/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl +++ b/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl @@ -46,6 +46,8 @@ function run_spin_precession!( ) where {T<:Real} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') + ρ = get_spin_property(p.ρ, seq.t[2:end]') + ρ_end = get_spin_property_at_end(ρ) #Effective field Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw ./ T(2π .* γ) #Rotation @@ -65,16 +67,16 @@ function run_spin_precession!( sig[:, :, 1] .= @views transpose(Mxy[:, findall(seq.ADC[2:end])]) if sim_method.save_Mz - Mz = M.z .* exp.(-tp' ./ p.T1) .+ p.ρ .* (1 .- exp.(-tp' ./ p.T1)) #Calculate intermediate points + Mz = M.z .* exp.(-tp' ./ p.T1) .+ ρ .* (1 .- exp.(-tp' ./ p.T1)) #Calculate intermediate points #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(Mz, seq.t[2:end]', p.motion; replace_by=p.ρ) + outflow_spin_reset!(Mz, seq.t[2:end]', p.motion; replace_by=ρ_end) sig[:, :, 2] .= @views transpose(Mz[:, findall(seq.ADC[2:end])]) #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 + M.z .= M.z .* exp.(-dur ./ p.T1) .+ ρ_end .* (1 .- exp.(-dur ./ p.T1)) #Jump to the last point end #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=ρ_end) return nothing end diff --git a/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl b/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl index 4e65121ae..8edde6930 100644 --- a/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl +++ b/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl @@ -99,6 +99,8 @@ function run_spin_excitation!( for i in eachindex(seq.Δt) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i + 1]) + ρ = get_spin_property(p.ρ, seq.t[i, :]) + ρ_end = get_spin_property_at_end(ρ) #Effective field @. ωxy_new = seq.B1[i + 1] * B_to_ω @. ωz_new = (seq.Gx[i + 1] * x + seq.Gy[i + 1] * y + seq.Gz[i + 1] * z + ΔBz) * B_to_ω + seq.Δf[i + 1] * T(2π) @@ -111,9 +113,9 @@ function run_spin_excitation!( mul!(Spinor(α, β), M, Maux_xy, Maux_z) #Relaxation @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2) - @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + p.ρ * (T(1) - exp(-seq.Δt[i] / p.T1)) + @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + ρ * (T(1) - exp(-seq.Δt[i] / p.T1)) #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=ρ_end) #Acquire signal if seq.ADC[i + 1] # ADC at the end of the time step sig[sample] = sum(M.xy) diff --git a/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl b/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl index 4c63d8b44..1e2c2ab67 100644 --- a/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl +++ b/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl @@ -33,6 +33,7 @@ function run_spin_precession!( ) where {T<:Real} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') + ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) #Effective field Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw ./ T(2π .* γ) #Rotation @@ -45,11 +46,12 @@ function run_spin_precession!( tp = cumsum(seq.Δt) # t' = t - t0 dur = sum(seq.Δt) # Total length, used for signal relaxation Mxy = M.xy .* exp.(-tp' ./ p.T2) .* cis.(ϕ) #This assumes Δw and T2 are constant in time + Mz = M.z .* exp.(-dur ./ p.T1) .+ ρ_end .* (1 .- exp.(-dur ./ p.T1)) M.xy .= Mxy[:, end] - M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1)) + M.z .= Mz[:, end] #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(Mxy, seq.t[2:end]', p.motion) - outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=ρ_end) #Acquired signal sig .= @views transpose(sum(Mxy[:, findall(seq.ADC[2:end])]; dims=1)) #<--- TODO: add coil sensitivities return nothing @@ -94,6 +96,8 @@ function run_spin_excitation!( ) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, s.t) + ρ = get_spin_property(p.ρ, s.t) + ρ_end = get_spin_property_at_end(ρ) #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(x,y,z) Bz = (s.Gx .* x .+ s.Gy .* y .+ s.Gz .* z) .+ ΔBz @@ -104,9 +108,9 @@ function run_spin_excitation!( 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) + ρ * (1 - exp(-s.Δt / p.T1)) #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, s.tnew, p.motion; replace_by=p.ρ) + outflow_spin_reset!(M, s.tnew, p.motion; replace_by=ρ_end) #Acquire signal # TODO: use sim_method and sys to modify sig if s.ADC # ADC at the end of the time step diff --git a/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl b/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl index 23b93c350..cdb4453cb 100644 --- a/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl +++ b/KomaMRICore/src/simulation/SimMethods/SimulationMethod.jl @@ -28,7 +28,7 @@ function initialize_spins_state( ) where {T<:Real} Nspins = length(obj) Mxy = zeros(T, Nspins) - Mz = obj.ρ + Mz = get_phantom_property(obj, :ρ) Xt = Mag{T}(Mxy, Mz) return Xt, obj end diff --git a/KomaMRIPlots/src/ui/DisplayFunctions.jl b/KomaMRIPlots/src/ui/DisplayFunctions.jl index c696674cf..f6b22e04f 100644 --- a/KomaMRIPlots/src/ui/DisplayFunctions.jl +++ b/KomaMRIPlots/src/ui/DisplayFunctions.jl @@ -1127,6 +1127,68 @@ function plot_kspace(seq::Sequence; width=nothing, height=nothing, darkmode=fals return plot_koma(p, l; config) end +motion_time_nodes(::NoMotion) = Float64[] +motion_time_nodes(motion) = times(motion) + +function append_time_dependent_property_nodes!(nodes, p::TimeDependentProperty) + append!(nodes, times(p)) + return nothing +end +append_time_dependent_property_nodes!(nodes, ::Any) = nothing + +function property_time_nodes(obj::Phantom) + nodes = Float64[] + for field in fieldnames(Phantom) + append_time_dependent_property_nodes!(nodes, getfield(obj, field)) + end + return nodes +end + +function phantom_time_nodes(obj::Phantom) + return unique(sort([motion_time_nodes(obj.motion); property_time_nodes(obj)])) +end + +phantom_time_node_count(obj::Phantom) = length(phantom_time_nodes(obj)) + +property_extrema(p::AbstractVector) = extrema(p) +property_extrema(p::TimeDependentProperty) = extrema(p.value) + +function spin_coords_for_plot(mv::NoMotion, x, y, z, t) + nt = length(t) + nt == 1 && return x, y, z + return repeat(x, 1, nt), repeat(y, 1, nt), repeat(z, 1, nt) +end +spin_coords_for_plot(motion, x, y, z, t) = get_spin_coords(motion, x, y, z, t') + +spin_property_frames(p::AbstractVector, t) = + repeat(reshape(KomaMRIBase.get_spin_property(p, t'), :, 1), 1, length(t)) +spin_property_frames(p::TimeDependentProperty, t) = KomaMRIBase.get_spin_property(p, t') + +function process_phantom_times(time_nodes, time_samples, ::Type{T}) where {T} + if isempty(time_nodes) + return [zero(T)] + elseif length(time_nodes) > 1 + time_min, time_max = minimum(time_nodes), maximum(time_nodes) + t = collect(range(time_min, time_max, length=time_samples)) + merged_times = sort(unique([t; time_nodes])) + if length(merged_times) > time_samples + to_keep = time_nodes + remaining = setdiff(merged_times, to_keep) + num_needed = time_samples - length(to_keep) + extra = remaining[round.(Int, range(1, length(remaining), length=num_needed))] + merged_times = sort([to_keep; extra]) + elseif length(merged_times) < time_samples + extra_points = setdiff(t, merged_times) + extra_needed = time_samples - length(merged_times) + additional = sort(extra_points)[1:extra_needed] + merged_times = sort([merged_times; additional]) + end + return merged_times + else + return time_nodes + end +end + """ p = plot_phantom_map(obj::Phantom, key::Symbol; kwargs...) @@ -1171,37 +1233,9 @@ function plot_phantom_map( view_2d=sum(KomaMRIBase.get_dims(obj)) < 3, colorbar=true, max_spins=20_000, - time_samples= obj.motion isa NoMotion ? 0 : length(times(obj.motion)), + time_samples=phantom_time_node_count(obj), kwargs..., ) - function process_times(::NoMotion) - return [zero(eltype(obj.x))] - end - - function process_times(motion) - time_nodes = times(motion) - if length(time_nodes) > 1 - time_min, time_max = minimum(time_nodes), maximum(time_nodes) - t = collect(range(time_min, time_max, length=time_samples)) - merged_times = sort(unique([t; time_nodes])) - if length(merged_times) > time_samples - to_keep = time_nodes - remaining = setdiff(merged_times, to_keep) - num_needed = time_samples - length(to_keep) - extra = remaining[round.(Int, range(1, length(remaining), length=num_needed))] - merged_times = sort([to_keep; extra]) - elseif length(merged_times) < time_samples - extra_points = setdiff(t, merged_times) - extra_needed = time_samples - length(merged_times) - additional = sort(extra_points)[1:extra_needed] - merged_times = sort([merged_times; additional]) - end - return merged_times - else - return time_nodes - end - end - function decimate_uniform_phantom(obj, num_points::Int) ss = Int(ceil(length(obj) / num_points)) return obj[1:ss:end] @@ -1213,8 +1247,8 @@ function plot_phantom_map( end path = @__DIR__ - cmin_key = minimum(getproperty(obj, key)) - cmax_key = maximum(getproperty(obj, key)) + prop = getproperty(obj, key) + cmin_key, cmax_key = property_extrema(prop) if key == :T1 || key == :T2 || key == :T2s cmin_key = 0 factor = 1e3 @@ -1253,8 +1287,9 @@ function plot_phantom_map( cmin_key = get(kwargs, :zmin, factor * cmin_key) cmax_key = get(kwargs, :zmax, factor * cmax_key) - t = process_times(obj.motion) - x, y, z = get_spin_coords(obj.motion, obj.x, obj.y, obj.z, t') + t = process_phantom_times(phantom_time_nodes(obj), time_samples, eltype(obj.x)) + x, y, z = spin_coords_for_plot(obj.motion, obj.x, obj.y, obj.z, t) + values = spin_property_frames(prop, t) x0 = -maximum(abs.([x y z])) * 1e2 xf = maximum(abs.([x y z])) * 1e2 @@ -1304,7 +1339,7 @@ function plot_phantom_map( x=dims[1] ? (x[:,i])*1e2 : (y[:,i])*1e2, y=dims[1] & dims[2] ? (y[:,i])*1e2 : (z[:,i])*1e2, mode="markers", - marker=attr(color=getproperty(obj,key)*factor, + marker=attr(color=values[:, i]*factor, showscale=colorbar, colorscale=colormap, colorbar=attr(ticksuffix=unit, title=string(key)), @@ -1314,7 +1349,7 @@ function plot_phantom_map( ), visible=i==1, showlegend=false, - text=round.(getproperty(obj,key)*factor,digits=4), + text=round.(values[:, i]*factor,digits=4), hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
$(string(key)): %{text}$unit")) end else # 3D @@ -1357,7 +1392,7 @@ function plot_phantom_map( y=(y[:,i])*1e2, z=(z[:,i])*1e2, mode="markers", - marker=attr(color=getproperty(obj,key)*factor, + marker=attr(color=values[:, i]*factor, showscale=colorbar, colorscale=colormap, colorbar=attr(ticksuffix=unit, title=string(key)), @@ -1367,7 +1402,7 @@ function plot_phantom_map( ), visible=i==1, showlegend=false, - text=round.(getproperty(obj,key)*factor,digits=4), + text=round.(values[:, i]*factor,digits=4), hovertemplate="x: %{x:.1f} cm
y: %{y:.1f} cm
z: %{z:.1f} cm
$(string(key)): %{text}$unit")) end end From b1f95f67e43684a487416fd04685efbccfb59f44 Mon Sep 17 00:00:00 2001 From: Pablo Villacorta Aylagas Date: Wed, 27 May 2026 13:45:07 +0200 Subject: [PATCH 2/2] Extend variable rho to remaining Phantom properties. Corresponding tests. --- KomaMRIBase/src/KomaMRIBase.jl | 7 ++- KomaMRIBase/src/datatypes/Phantom.jl | 58 +++++++++++-------- .../phantom/TimeDependentProperty.jl | 37 +++++++++++- KomaMRIBase/src/motion/MotionList.jl | 8 +-- KomaMRIBase/src/motion/NoMotion.jl | 10 ++-- KomaMRIBase/src/timing/Interpolation.jl | 4 ++ KomaMRIBase/test/runtests.jl | 28 +++++---- .../SimMethods/Bloch/cpu/BlochCPU.jl | 27 ++++++--- .../SimMethods/Bloch/gpu/BlochGPU.jl | 16 ++--- .../SimMethods/BlochDict/BlochDict.jl | 13 +++-- .../BlochMagnus/cpu/BlochMagnusCPU.jl | 10 ++-- .../SimMethods/BlochSimple/BlochSimple.jl | 18 +++--- KomaMRIFiles/src/Phantom/Phantom.jl | 58 ++++++++++++++++--- KomaMRIFiles/test/runtests.jl | 12 ++++ 14 files changed, 220 insertions(+), 86 deletions(-) diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl index a13dcaf35..e0808bfa1 100644 --- a/KomaMRIBase/src/KomaMRIBase.jl +++ b/KomaMRIBase/src/KomaMRIBase.jl @@ -106,8 +106,11 @@ export Rotate, RotateX, RotateY, RotateZ, CenterOfMass export HeartBeat, Path, FlowPath export TimeRange, Periodic, TimeCurve export SpinRange, AllSpins, spin_indicator -export get_spin_coords, get_spin_property, get_phantom_property, get_spin_property_at_end -export TimeDependentProperty +export get_spin_coords +export TimeDependentProperty, get_phantom_property +export get_spin_property, get_spin_property_at_end +export get_spin_properties, get_spin_properties_block_end + # Secondary export get_kspace, rotx, roty, rotz # Additionals diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl index 42758e1b2..8236cfa5e 100644 --- a/KomaMRIBase/src/datatypes/Phantom.jl +++ b/KomaMRIBase/src/datatypes/Phantom.jl @@ -11,14 +11,14 @@ a property value representing a spin. This struct serves as an input for the sim - `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 -- `ρ`: (`::Union{AbstractVector{T}, TimeDependentProperty{T}}`) spin proton density -- `T1`: (`::AbstractVector{T<:Real}`, `[s]`) spin T1 parameter vector -- `T2`: (`::AbstractVector{T<:Real}`, `[s]`) spin T2 parameter vector -- `T2s`: (`::AbstractVector{T<:Real}`, `[s]`) spin T2s parameter vector -- `Δw`: (`::AbstractVector{T<:Real}`, `[rad/s]`) spin off-resonance parameter vector -- `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 +- `ρ`: (`::PhantomProperty{T}`) spin proton density +- `T1`: (`::PhantomProperty{T}`, `[s]`) spin T1 parameter vector +- `T2`: (`::PhantomProperty{T}`, `[s]`) spin T2 parameter vector +- `T2s`: (`::PhantomProperty{T}`, `[s]`) spin T2s parameter vector +- `Δw`: (`::PhantomProperty{T}`, `[rad/s]`) spin off-resonance parameter vector +- `Dλ1`: (`::PhantomProperty{T}`) spin Dλ1 (diffusion) parameter vector +- `Dλ2`: (`::PhantomProperty{T}`) spin Dλ2 (diffusion) parameter vector +- `Dθ`: (`::PhantomProperty{T}`) spin Dθ (diffusion) parameter vector - `motion`: (`::Union{NoMotion, Motion{T<:Real} MotionList{T<:Real}}`) motion # Returns @@ -31,29 +31,33 @@ julia> obj = Phantom(x=[0.0]) julia> obj.ρ ``` """ +const PhantomProperty{T<:Real} = Union{AbstractVector{T}, TimeDependentProperty{T}} + @with_kw mutable struct Phantom{T<:Real} name::String = "spins" x::AbstractVector{T} = @isdefined(T) ? T[] : Float64[] y::AbstractVector{T} = zeros(eltype(x), size(x)) z::AbstractVector{T} = zeros(eltype(x), size(x)) - ρ::Union{AbstractVector{T}, TimeDependentProperty{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 + ρ::PhantomProperty{T} = ones(eltype(x), size(x)) + T1::PhantomProperty{T} = ones(eltype(x), size(x)) * 1_000_000 + T2::PhantomProperty{T} = ones(eltype(x), size(x)) * 1_000_000 + T2s::PhantomProperty{T} = ones(eltype(x), size(x)) * 1_000_000 #Off-resonance related - Δw::AbstractVector{T} = zeros(eltype(x), size(x)) + Δw::PhantomProperty{T} = zeros(eltype(x), size(x)) #χ::Vector{SusceptibilityModel} #Diffusion - 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)) + Dλ1::PhantomProperty{T} = zeros(eltype(x), size(x)) + Dλ2::PhantomProperty{T} = zeros(eltype(x), size(x)) + Dθ::PhantomProperty{T} = zeros(eltype(x), size(x)) #Diff::Vector{DiffusionModel} #Diffusion map #Motion motion::Union{NoMotion, Motion{T}, MotionList{T}} = NoMotion() end -const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String, fieldnames(Phantom)) -const VECTOR_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) <: AbstractVector, fieldnames(Phantom)) +const POSITION_PHANTOM_FIELDS = (:x, :y, :z) +const PROPERTY_PHANTOM_FIELDS = (:ρ, :T1, :T2, :T2s, :Δw, :Dλ1, :Dλ2, :Dθ) + +const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String, fieldnames(Phantom)) get_phantom_property(obj::Phantom, key::Symbol) = get_property_value(getfield(obj, key)) get_property_value(p) = p @@ -100,22 +104,28 @@ end """Addition of phantoms""" +(obj1::Phantom, obj2::Phantom) = begin name = first(obj1.name * "+" * obj2.name, 50) # The name is limited to 50 characters - fields = [] - for field in VECTOR_PHANTOM_FIELDS - push!(fields, (field, [getfield(obj1, field); getfield(obj2, field)])) + fields = Pair{Symbol,Any}[] + for field in POSITION_PHANTOM_FIELDS + push!(fields, field => [getfield(obj1, field); getfield(obj2, field)]) + end + for field in PROPERTY_PHANTOM_FIELDS + push!(fields, field => vcat_property(getfield(obj1, field), getfield(obj2, field))) end return Phantom(; name=name, fields..., - ρ=vcat_property(getfield(obj1, :ρ), getfield(obj2, :ρ)), - motion=vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)), + motion=vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)), ) end """Scalar multiplication of a phantom""" *(α::Real, obj::Phantom) = begin obj1 = deepcopy(obj) - obj1.ρ .*= α + if obj1.ρ isa AbstractVector + obj1.ρ .*= α + else + obj1.ρ = α * obj1.ρ + end return obj1 end diff --git a/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl b/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl index 704129b25..c9dbb84a4 100644 --- a/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl +++ b/KomaMRIBase/src/datatypes/phantom/TimeDependentProperty.jl @@ -5,6 +5,35 @@ end get_property_value(p::TimeDependentProperty) = p.value[:, 1] +"""Spin-resolved phantom properties at simulation times `t` (vector or matrix columns).""" +function get_spin_properties(p, t) + return ( + ρ=get_spin_property(p.ρ, t), + T1=get_spin_property(p.T1, t), + T2=get_spin_property(p.T2, t), + T2s=get_spin_property(p.T2s, t), + Δw=get_spin_property(p.Δw, t), + Dλ1=get_spin_property(p.Dλ1, t), + Dλ2=get_spin_property(p.Dλ2, t), + Dθ=get_spin_property(p.Dθ, t), + ) +end + +"""Block-end spin properties used for GPU kernels and block-final relaxation.""" +function get_spin_properties_block_end(p, t) + props = get_spin_properties(p, t) + return ( + ρ=get_spin_property_at_end(props.ρ), + T1=get_spin_property_at_end(props.T1), + T2=get_spin_property_at_end(props.T2), + T2s=get_spin_property_at_end(props.T2s), + Δw=get_spin_property_at_end(props.Δw), + Dλ1=get_spin_property_at_end(props.Dλ1), + Dλ2=get_spin_property_at_end(props.Dλ2), + Dθ=get_spin_property_at_end(props.Dθ), + ) +end + get_spin_property(p::AbstractVector, t) = p function get_spin_property(p::TimeDependentProperty, t) Ns = size(p.value, 1) @@ -12,7 +41,13 @@ function get_spin_property(p::TimeDependentProperty, t) return resample(interpolate(p.value, Gridded(Linear()), Val(Ns), t_unit), t_unit) end -get_spin_property_at_end(p::AbstractMatrix) = selectdim(p, 2, size(p, 2)) +""" + get_spin_property_at_end(p::AbstractMatrix) + +Returns spin-wise property values at the last queried time in `p`, with shape +`(Nspins, Ntimes)`. +""" +get_spin_property_at_end(p::AbstractMatrix) = vec(selectdim(p, 2, size(p, 2))) get_spin_property_at_end(p::AbstractVector) = p Base.:(==)(p1::TimeDependentProperty, p2::TimeDependentProperty) = p1.value == p2.value && p1.time == p2.time diff --git a/KomaMRIBase/src/motion/MotionList.jl b/KomaMRIBase/src/motion/MotionList.jl index 207a173a7..941bdad5f 100644 --- a/KomaMRIBase/src/motion/MotionList.jl +++ b/KomaMRIBase/src/motion/MotionList.jl @@ -52,7 +52,7 @@ end # https://github.com/JuliaHealth/KomaMRI.jl/issues/480 """ Addition of MotionLists """ # MotionList + MotionList -function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real} +function vcat_motion(m1::MotionList{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real} mv_aux = Motion{T}[] for m in m1.motions m_aux = deepcopy(m) @@ -68,7 +68,7 @@ function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Rea return MotionList(mv_aux...) end # Motion + Motion -function Base.vcat(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} +function vcat_motion(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} mv_aux = Motion{T}[] m_aux = deepcopy(m1) m_aux.spins = expand(m_aux.spins, Ns1) @@ -80,8 +80,8 @@ function Base.vcat(m1::Motion{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} return MotionList(mv_aux...) end # Motion + MotionList -Base.vcat(m1::MotionList{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} = vcat(m2, m1, Ns2, Ns1) -function Base.vcat(m1::Motion{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real} +vcat_motion(m1::MotionList{T}, m2::Motion{T}, Ns1, Ns2) where {T<:Real} = vcat_motion(m2, m1, Ns2, Ns1) +function vcat_motion(m1::Motion{T}, m2::MotionList{T}, Ns1, Ns2) where {T<:Real} mv_aux = Motion{T}[] m_aux = deepcopy(m1) m_aux.spins = expand(m_aux.spins, Ns1) diff --git a/KomaMRIBase/src/motion/NoMotion.jl b/KomaMRIBase/src/motion/NoMotion.jl index 3c5561c43..9a639bdbd 100644 --- a/KomaMRIBase/src/motion/NoMotion.jl +++ b/KomaMRIBase/src/motion/NoMotion.jl @@ -18,10 +18,10 @@ Base.view(mv::NoMotion, p) = mv """ Addition of NoMotions """ # NoMotion + NoMotion -Base.vcat(m1::NoMotion, m2::NoMotion, Ns1, Ns2) = m1 +vcat_motion(m1::NoMotion, m2::NoMotion, Ns1, Ns2) = m1 # NoMotion + MotionList -Base.vcat(m1::MotionList, m2::NoMotion, Ns1, Ns2) = vcat(m2, m1, 0, Ns1) -function Base.vcat(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T} +vcat_motion(m1::MotionList, m2::NoMotion, Ns1, Ns2) = vcat_motion(m2, m1, 0, Ns1) +function vcat_motion(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T} mv_aux = Motion{T}[] for m in m2.motions m_aux = deepcopy(m) @@ -32,8 +32,8 @@ function Base.vcat(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T} return MotionList(mv_aux) end # NoMotion + Motion -Base.vcat(m1::Motion, m2::NoMotion, Ns1, Ns2) = vcat(m2, m1, 0, Ns1) -function Base.vcat(m1::NoMotion, m2::Motion{T}, Ns1, Ns2) where {T} +vcat_motion(m1::Motion, m2::NoMotion, Ns1, Ns2) = vcat_motion(m2, m1, 0, Ns1) +function vcat_motion(m1::NoMotion, m2::Motion{T}, Ns1, Ns2) where {T} m_aux = deepcopy(m2) m_aux.spins = expand(m_aux.spins, Ns2) m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1) diff --git a/KomaMRIBase/src/timing/Interpolation.jl b/KomaMRIBase/src/timing/Interpolation.jl index 09b699867..1227ca68d 100644 --- a/KomaMRIBase/src/timing/Interpolation.jl +++ b/KomaMRIBase/src/timing/Interpolation.jl @@ -54,6 +54,10 @@ function resample(itp::Interpolator2D, t) Ns = size(itp.coefs, 1) id = _similar(t, Ns) copyto!(id, collect(range(oneunit(eltype(t)), eltype(t)(Ns), Ns))) + if t isa AbstractVector + # Produce an outer-product evaluation: Ns spins × length(t) times. + return itp.(reshape(id, :, 1), reshape(t, 1, :)) + end return itp.(id, t) end diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl index 1677233e8..e3a7c9067 100644 --- a/KomaMRIBase/test/runtests.jl +++ b/KomaMRIBase/test/runtests.jl @@ -1527,7 +1527,7 @@ end 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]], - vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) ) # NOTE: these vcat methods must be simplified once the Vector{<:Motion} approach is accomplished: # https://github.com/JuliaHealth/KomaMRI.jl/issues/480 @@ -1535,42 +1535,42 @@ end @test obj1 + obj2 == oba # NoMotion + MotionList obj2.motion = MotionList(tr, rt) - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # MotionList + NoMotion obj1.motion = MotionList(tr, rt) obj2.motion = NoMotion() - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # NoMotion + Motion obj1.motion = NoMotion() obj2.motion = tr - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # Motion + NoMotion obj1.motion = tr obj2.motion = NoMotion() - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # MotionList + MotionList obj1.motion = MotionList(tr, rt) obj2.motion = MotionList(tr, rt) - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # Motion + Motion obj1.motion = tr obj2.motion = rt - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # Motion + MotionList obj1.motion = tr obj2.motion = MotionList(tr, rt) - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba # MotionList + Motion obj1.motion = MotionList(tr, rt) obj2.motion = tr - oba.motion = vcat(obj1.motion, obj2.motion, length(obj1), length(obj2)) + oba.motion = KomaMRIBase.vcat_motion(obj1.motion, obj2.motion, length(obj1), length(obj2)) @test obj1 + obj2 == oba tc = TimeRange(0.0, 1.0) @@ -1606,7 +1606,15 @@ end @test KomaMRIBase.get_spin_property(obj.ρ, [0.5]) ≈ [0.5, 0.75] ρt2 = KomaMRIBase.get_spin_property(obj.ρ, [0.0, 1.0]) @test ρt2 ≈ [1.0 0.0; 0.5 1.0] - @test KomaMRIBase.get_spin_property_at_end(ρt2) ≈ [1.0, 1.0] + @test KomaMRIBase.get_spin_property_at_end(ρt2) ≈ [0.0, 1.0] + end + @testset "Time-dependent T1" begin + tc = TimeRange(0.0, 1.0) + T1t = TimeDependentProperty([0.5 1.0; 1.5 2.0], tc) + obj = Phantom(T1=T1t, x=x[1:2], y=y[1:2], z=z[1:2]) + props = KomaMRIBase.get_spin_properties(obj, [0.5]) + @test props.T1 ≈ [0.75, 1.75] + @test KomaMRIBase.get_spin_properties_block_end(obj, [0.0, 1.0]).T1 ≈ [1.0, 2.0] end @testset "Brain Phantom 2D" begin ph = brain_phantom2D() diff --git a/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl b/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl index 910c8f5ab..f70009e30 100644 --- a/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl +++ b/KomaMRICore/src/simulation/SimMethods/Bloch/cpu/BlochCPU.jl @@ -33,7 +33,7 @@ function prealloc(sim_method::Bloch, backend::KA.CPU, obj::Phantom{T}, M::Mag{T} similar(M.xy), similar(M.xy) ), - obj.Δw ./ T(2π .* γ) + zeros(T, size(obj.x)), ) end @@ -65,22 +65,28 @@ function run_spin_precession!( fill!(ϕ, zero(T)) block_time = zero(T) sample = 1 - ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) + props_end = get_spin_properties_block_end(p, seq.t') + ρ_end = props_end.ρ x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[1]) + Δw = get_spin_property(p.Δw, seq.t[1]) + @. ΔBz = Δw / T(2π * γ) @. Bz_old = x * seq.Gx[1] + y * seq.Gy[1] + z * seq.Gz[1] + ΔBz #Simulation for i in eachindex(seq.Δt) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i + 1]) #Effective Field + Δw = get_spin_property(p.Δw, seq.t[i + 1]) + @. ΔBz = Δw / T(2π * γ) @. Bz_new = x * seq.Gx[i + 1] + y * seq.Gy[i + 1] + z * seq.Gz[i + 1] + ΔBz #Rotation @. ϕ += (Bz_old + Bz_new) * T(-π * γ) * seq.Δt[i] block_time += seq.Δt[i] #Acquired Signal if seq.ADC[i + 1] + T2 = get_spin_property(p.T2, seq.t[i + 1]) #Update signal - @. Mxy = exp(-block_time / p.T2) * M.xy * cis(ϕ) + @. Mxy = exp(-block_time / T2) * M.xy * cis(ϕ) #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(Mxy, seq.t[i + 1], p.motion) #Acquired signal @@ -91,8 +97,9 @@ function run_spin_precession!( Bz_old .= Bz_new end #Final Spin-State - @. M.xy = M.xy * exp(-block_time / p.T2) * cis(ϕ) - @. M.z = M.z * exp(-block_time / p.T1) + ρ_end * (T(1) - exp(-block_time / p.T1)) + T2_end = props_end.T2 + @. M.xy = M.xy * exp(-block_time / T2_end) * cis(ϕ) + @. M.z = M.z * exp(-block_time / props_end.T1) + ρ_end * (T(1) - exp(-block_time / props_end.T1)) #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end) return nothing @@ -134,10 +141,12 @@ function run_spin_excitation!( for i in eachindex(seq.Δt) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i]) - ρ = get_spin_property(p.ρ, seq.t[i, :]) + props = get_spin_properties(p, seq.t[i, :]) + ρ = props.ρ ρ_end = get_spin_property_at_end(ρ) #Effective field - @. Bz = (seq.Gx[i] * x + seq.Gy[i] * y + seq.Gz[i] * z) + ΔBz - seq.Δf[i] / T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z) + @. ΔBz = props.Δw / T(2π * γ) - seq.Δf[i] / T(γ) # ΔB_0 = (B_0 - ω_rf/γ), Need to add a component here to model scanner's dB0(x,y,z) + @. Bz = seq.Gx[i] * x + seq.Gy[i] * y + seq.Gz[i] * z + ΔBz @. B = sqrt(abs(seq.B1[i])^2 + abs(Bz)^2) @. B[B == 0] = eps(T) #Spinor Rotation @@ -146,8 +155,8 @@ function run_spin_excitation!( @. β = -Complex{T}(im) * (seq.B1[i] / B) * sin(φ_half) mul!(Spinor(α, β), M, Maux_xy, Maux_z) #Relaxation - @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2) - @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + ρ * (T(1) - exp(-seq.Δt[i] / p.T1)) + @. M.xy = M.xy * exp(-seq.Δt[i] / props.T2) + @. M.z = M.z * exp(-seq.Δt[i] / props.T1) + ρ * (T(1) - exp(-seq.Δt[i] / props.T1)) #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=ρ_end) #Acquire signal diff --git a/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl b/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl index d6e7f5259..139f3443c 100644 --- a/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl +++ b/KomaMRICore/src/simulation/SimMethods/Bloch/gpu/BlochGPU.jl @@ -21,7 +21,7 @@ function prealloc( return BlochGPUPrealloc( KA.zeros(backend, Complex{T}, (cld(size(obj.x, 1), groupsize), max_block_length)), KA.zeros(backend, Complex{T}, 1, max_block_length), - obj.Δw ./ T(2π .* γ) + obj.Δw ./ T(2π .* γ), ) end @@ -37,14 +37,15 @@ function run_spin_precession!( ) where {T<:Real, SM<:BlochLikeSimMethods} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') - ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) + props_end = get_spin_properties_block_end(p, seq.t') + @. pre.ΔBz = props_end.Δw / T(2π * γ) has_adc = length(sig) > 0 #Precession precession_kernel!(backend, groupsize)( pre.sig_output, M.xy, M.z, - x, y, z, pre.ΔBz, p.T1, p.T2, ρ_end, UInt32(length(M.xy)), + x, y, z, pre.ΔBz, props_end.T1, props_end.T2, props_end.ρ, UInt32(length(M.xy)), seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.ADC, UInt32(length(seq.t)), Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), sim_method, @@ -58,7 +59,7 @@ function run_spin_precession!( end #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end) + outflow_spin_reset!(M, seq.t', p.motion; replace_by=props_end.ρ) return nothing end @@ -75,14 +76,15 @@ function run_spin_excitation!( ) where {T<:Real, SM<:BlochLikeSimMethods} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') - ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) + props_end = get_spin_properties_block_end(p, seq.t') + @. pre.ΔBz = props_end.Δw / T(2π * γ) has_adc = length(sig) > 0 #Excitation excitation_kernel!(backend, groupsize)( pre.sig_output, M.xy, M.z, - x, y, z, pre.ΔBz, p.T1, p.T2, ρ_end, UInt32(length(M.xy)), + x, y, z, pre.ΔBz, props_end.T1, props_end.T2, props_end.ρ, UInt32(length(M.xy)), seq.Gx, seq.Gy, seq.Gz, seq.Δt, seq.Δf, seq.B1, seq.ψ, seq.ADC, UInt32(length(seq.t)), Val(!(p.motion isa NoMotion)), Val(supports_warp_reduction(backend)), Val(has_adc), sim_method, @@ -96,7 +98,7 @@ function run_spin_excitation!( end #Reset Spin-State (Magnetization). Only for FlowPath - outflow_spin_reset!(M, seq.t', p.motion; replace_by=ρ_end) # TODO: reset state inside kernel + outflow_spin_reset!(M, seq.t', p.motion; replace_by=props_end.ρ) # TODO: reset state inside kernel return nothing end diff --git a/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl b/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl index 1209d31f5..52f7b9923 100644 --- a/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl +++ b/KomaMRICore/src/simulation/SimMethods/BlochDict/BlochDict.jl @@ -46,10 +46,12 @@ function run_spin_precession!( ) where {T<:Real} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') - ρ = get_spin_property(p.ρ, seq.t[2:end]') + t_acq = seq.t[2:end]' + props = get_spin_properties(p, t_acq) + ρ = props.ρ ρ_end = get_spin_property_at_end(ρ) #Effective field - Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw ./ T(2π .* γ) + Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ props.Δw ./ T(2π .* γ) #Rotation if is_ADC_on(seq) ϕ = T(-2π .* γ) .* cumtrapz(seq.Δt', Bz) @@ -59,7 +61,7 @@ function run_spin_precession!( #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 .* exp.(-tp' ./ p.T2) .* cis.(ϕ) #This assumes Δw and T2 are constant in time + Mxy = M.xy .* exp.(-tp' ./ props.T2) .* cis.(ϕ) M.xy .= Mxy[:, end] #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(Mxy, seq.t[2:end]', p.motion) @@ -67,13 +69,14 @@ function run_spin_precession!( sig[:, :, 1] .= @views transpose(Mxy[:, findall(seq.ADC[2:end])]) if sim_method.save_Mz - Mz = M.z .* exp.(-tp' ./ p.T1) .+ ρ .* (1 .- exp.(-tp' ./ p.T1)) #Calculate intermediate points + Mz = M.z .* exp.(-tp' ./ props.T1) .+ ρ .* (1 .- exp.(-tp' ./ props.T1)) #Calculate intermediate points #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(Mz, seq.t[2:end]', p.motion; replace_by=ρ_end) sig[:, :, 2] .= @views transpose(Mz[:, findall(seq.ADC[2:end])]) #Save state to signal M.z .= Mz[:, end] else - M.z .= M.z .* exp.(-dur ./ p.T1) .+ ρ_end .* (1 .- exp.(-dur ./ p.T1)) #Jump to the last point + T1_end = get_spin_property_at_end(props.T1) + M.z .= M.z .* exp.(-dur ./ T1_end) .+ ρ_end .* (1 .- exp.(-dur ./ T1_end)) #Jump to the last point end #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(M, seq.t[2:end]', p.motion; replace_by=ρ_end) diff --git a/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl b/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl index 8edde6930..97acd8e4c 100644 --- a/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl +++ b/KomaMRICore/src/simulation/SimMethods/BlochMagnus/cpu/BlochMagnusCPU.jl @@ -45,7 +45,7 @@ function prealloc(sim_method::BlochMagnus, backend::KA.CPU, obj::Phantom{T}, M:: similar(M.xy), similar(M.xy) ), - obj.Δw ./ T(2π .* γ) + zeros(T, size(obj.x)), ) end @@ -99,9 +99,11 @@ function run_spin_excitation!( for i in eachindex(seq.Δt) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t[i + 1]) - ρ = get_spin_property(p.ρ, seq.t[i, :]) + props = get_spin_properties(p, seq.t[i, :]) + ρ = props.ρ ρ_end = get_spin_property_at_end(ρ) #Effective field + @. ΔBz = props.Δw / T(2π * γ) @. ωxy_new = seq.B1[i + 1] * B_to_ω @. ωz_new = (seq.Gx[i + 1] * x + seq.Gy[i + 1] * y + seq.Gz[i + 1] * z + ΔBz) * B_to_ω + seq.Δf[i + 1] * T(2π) effective_rotation_vector!(θxy, θz, ωxy_old, ωz_old, ωxy_new, ωz_new, seq.Δt[i], sim_method) @@ -112,8 +114,8 @@ function run_spin_excitation!( @. β = -Complex{T}(im) * (θxy / θ) * sin(θ / T(2)) mul!(Spinor(α, β), M, Maux_xy, Maux_z) #Relaxation - @. M.xy = M.xy * exp(-seq.Δt[i] / p.T2) - @. M.z = M.z * exp(-seq.Δt[i] / p.T1) + ρ * (T(1) - exp(-seq.Δt[i] / p.T1)) + @. M.xy = M.xy * exp(-seq.Δt[i] / props.T2) + @. M.z = M.z * exp(-seq.Δt[i] / props.T1) + ρ * (T(1) - exp(-seq.Δt[i] / props.T1)) #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(M, seq.t[i + 1, :], p.motion; replace_by=ρ_end) #Acquire signal diff --git a/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl b/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl index 1e2c2ab67..9fc2899ec 100644 --- a/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl +++ b/KomaMRICore/src/simulation/SimMethods/BlochSimple/BlochSimple.jl @@ -33,9 +33,10 @@ function run_spin_precession!( ) where {T<:Real} #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, seq.t') - ρ_end = get_spin_property_at_end(get_spin_property(p.ρ, seq.t')) + props = get_spin_properties(p, seq.t') + ρ_end = get_spin_property_at_end(props.ρ) #Effective field - Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ p.Δw ./ T(2π .* γ) + Bz = x .* seq.Gx' .+ y .* seq.Gy' .+ z .* seq.Gz' .+ props.Δw ./ T(2π .* γ) #Rotation if is_ADC_on(seq) ϕ = T(-2π .* γ) .* cumtrapz(seq.Δt', Bz) @@ -45,8 +46,8 @@ function run_spin_precession!( #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 .* exp.(-tp' ./ p.T2) .* cis.(ϕ) #This assumes Δw and T2 are constant in time - Mz = M.z .* exp.(-dur ./ p.T1) .+ ρ_end .* (1 .- exp.(-dur ./ p.T1)) + Mxy = M.xy .* exp.(-tp' ./ props.T2) .* cis.(ϕ) + Mz = M.z .* exp.(-dur ./ get_spin_property_at_end(props.T1)) .+ ρ_end .* (1 .- exp.(-dur ./ get_spin_property_at_end(props.T1))) M.xy .= Mxy[:, end] M.z .= Mz[:, end] #Reset Spin-State (Magnetization). Only for FlowPath @@ -96,10 +97,11 @@ function run_spin_excitation!( ) #Motion x, y, z = get_spin_coords(p.motion, p.x, p.y, p.z, s.t) - ρ = get_spin_property(p.ρ, s.t) + props = get_spin_properties(p, s.t) + ρ = props.ρ ρ_end = get_spin_property_at_end(ρ) #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(x,y,z) + ΔBz = props.Δ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) @@ -107,8 +109,8 @@ function run_spin_excitation!( φ = T(-2π .* γ) .* (B .* s.Δt) 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) + ρ * (1 - exp(-s.Δt / p.T1)) + @. M.xy = M.xy * exp(-s.Δt / props.T2) + @. M.z = M.z * exp(-s.Δt / props.T1) + ρ * (1 - exp(-s.Δt / props.T1)) #Reset Spin-State (Magnetization). Only for FlowPath outflow_spin_reset!(M, s.tnew, p.motion; replace_by=ρ_end) #Acquire signal diff --git a/KomaMRIFiles/src/Phantom/Phantom.jl b/KomaMRIFiles/src/Phantom/Phantom.jl index 3589b7f8d..ba3c7eeee 100644 --- a/KomaMRIFiles/src/Phantom/Phantom.jl +++ b/KomaMRIFiles/src/Phantom/Phantom.jl @@ -19,12 +19,14 @@ function read_phantom(filename::String) name = read_attribute(fid, "Name") push!(phantom_fields, (:name, name)) # Position and contrast - for key in ["position", "contrast"] - group = fid[key] - for label in HDF5.keys(group) - values = read(group[label]) - push!(phantom_fields, (Symbol(label), values)) - end + pos = fid["position"] + for label in HDF5.keys(pos) + push!(phantom_fields, (Symbol(label), read(pos[label]))) + end + contrast = fid["contrast"] + T = eltype(phantom_fields[end][2]) + for label in HDF5.keys(contrast) + import_contrast_property!(phantom_fields, contrast, label, T) end # Motion if "motion" in keys(fid) @@ -88,6 +90,48 @@ function import_motion_subfield!(motion_subfields::Array, subfield_value::Union{ push!(motion_subfields, subfield_value) return nothing end +function read_timecurve_group(group::HDF5.Group) + t = read(group, "t") + t_unit = read(group, "t_unit") + periodic = haskey(HDF5.attributes(group), "periodic") ? read_attribute(group, "periodic") : false + periods = haskey(HDF5.attributes(group), "periods") ? read_attribute(group, "periods") : 1.0 + return TimeCurve(; t, t_unit, periodic, periods) +end + +function write_timecurve_group!(group::HDF5.Group, tc::TimeCurve) + HDF5.attributes(group)["type"] = "TimeCurve" + group["t"] = tc.t + group["t_unit"] = tc.t_unit + HDF5.attributes(group)["periodic"] = tc.periodic + HDF5.attributes(group)["periods"] = tc.periods +end + +function export_contrast_property!(contrast::HDF5.Group, label::Symbol, prop) + if prop isa TimeDependentProperty + grp = create_group(contrast, string(label)) + HDF5.attributes(grp)["type"] = "TimeDependentProperty" + grp["value"] = prop.value + time_grp = create_group(grp, "time") + write_timecurve_group!(time_grp, prop.time) + else + contrast[string(label)] = prop + end + return nothing +end + +function import_contrast_property!(phantom_fields::Array, contrast::HDF5.Group, label::String, ::Type{T}) where {T<:Real} + node = contrast[label] + if node isa HDF5.Dataset + push!(phantom_fields, (Symbol(label), read(node))) + else + @assert read_attribute(node, "type") == "TimeDependentProperty" + value = read(node["value"]) + time = read_timecurve_group(node["time"]) + push!(phantom_fields, (Symbol(label), TimeDependentProperty(value, time))) + end + return nothing +end + function import_motion_subfield!(motion_subfields::Array, subfield_value::String, key::String, T::Type{<:Real}) if subfield_value in ["true", "false"] return push!(motion_subfields, subfield_value == "true" ? true : false) @@ -132,7 +176,7 @@ function write_phantom( # Contrast (Rho, T1, T2, T2s Deltaw) contrast = create_group(fid, "contrast") for x in store_contrasts - contrast[String(x)] = getfield(obj, x) + export_contrast_property!(contrast, x, getfield(obj, x)) end # Motion if !(obj.motion isa NoMotion) & store_motion diff --git a/KomaMRIFiles/test/runtests.jl b/KomaMRIFiles/test/runtests.jl index 0a1a935c4..24d5fb5c9 100644 --- a/KomaMRIFiles/test/runtests.jl +++ b/KomaMRIFiles/test/runtests.jl @@ -41,6 +41,18 @@ end obj2 = read_phantom(filename) @test obj1 == obj2 end + @testset "TimeDependentProperty" begin + pth = @__DIR__ + filename = pth * "/test_files/phantom/brain_timedep_w.phantom" + obj1 = brain_phantom2D() + tc = TimeRange(0.0, 1.0) + Ns = length(obj1) + obj1.ρ = TimeDependentProperty([obj1.ρ obj1.ρ], tc) + obj1.T1 = TimeDependentProperty([obj1.T1 2 .* obj1.T1], tc) + write_phantom(obj1, filename) + obj2 = read_phantom(filename) + @test obj1 == obj2 + end @testset "ArbitraryAction" begin pth = @__DIR__ filename = pth * "/test_files/phantom/brain_arbitrarymotion_w.phantom"