diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl index 370007e1a..ca76a7462 100644 --- a/KomaMRIBase/src/KomaMRIBase.jl +++ b/KomaMRIBase/src/KomaMRIBase.jl @@ -46,10 +46,12 @@ export kfoldperm, trapz, cumtrapz # Phantom export brain_phantom2D, brain_phantom3D, pelvis_phantom2D, heart_phantom # Motion -export AbstractMotion, Motion, MotionVector, NoMotion +export MotionList, NoMotion export SimpleMotion, ArbitraryMotion -export Translation, Rotation, HeartBeat, Trajectory, FlowTrajectory -export TimeScale, TimeRange, Periodic +export Translation, TranslationX, TranslationY, TranslationZ +export Rotation, RotationX, RotationY, RotationZ +export HeartBeat, Trajectory, FlowTrajectory +export AbstractTimeSpan, TimeRange, Periodic export sort_motions!, get_spin_coords # Secondary export get_kspace, rotx, roty, rotz diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl index 4de25d2df..9f5cee974 100644 --- a/KomaMRIBase/src/datatypes/Phantom.jl +++ b/KomaMRIBase/src/datatypes/Phantom.jl @@ -1,7 +1,7 @@ # TimeScale: include("../timing/TimeScale.jl") # Motion: -abstract type AbstractMotion{T<:Real} end +abstract type AbstractMotionList{T<:Real} end include("phantom/Motion.jl") include("phantom/NoMotion.jl") @@ -54,7 +54,7 @@ julia> obj.ρ Dθ::AbstractVector{T} = zeros(eltype(x), size(x)) #Diff::Vector{DiffusionModel} #Diffusion map #Motion - motion::AbstractMotion{T} = NoMotion{eltype(x)}() + motion::AbstractMotionList{T} = NoMotion{eltype(x)}() end """Size and length of a phantom""" @@ -192,15 +192,15 @@ function heart_phantom( Dλ1=Dλ1[ρ .!= 0], Dλ2=Dλ2[ρ .!= 0], Dθ=Dθ[ρ .!= 0], - motion=MotionVector( + motion=MotionList( HeartBeat(; - times=Periodic(; period=period, asymmetry=asymmetry), + time=Periodic(; period=period, asymmetry=asymmetry), circumferential_strain=circumferential_strain, radial_strain=radial_strain, longitudinal_strain=0.0, ), Rotation(; - times=Periodic(; period=period, asymmetry=asymmetry), + time=Periodic(; period=period, asymmetry=asymmetry), yaw=rotation_angle, pitch=0.0, roll=0.0, diff --git a/KomaMRIBase/src/datatypes/phantom/Motion.jl b/KomaMRIBase/src/datatypes/phantom/Motion.jl index eafde6d02..aa5b46a3d 100644 --- a/KomaMRIBase/src/datatypes/phantom/Motion.jl +++ b/KomaMRIBase/src/datatypes/phantom/Motion.jl @@ -1,53 +1,53 @@ -abstract type Motion{T<:Real} end +abstract type AbstractMotion{T<:Real} end -is_composable(m::Motion) = false +is_composable(m::AbstractMotion) = false -struct MotionVector{T<:Real} <: AbstractMotion{T} - motions::Vector{<:Motion{T}} +struct MotionList{T<:Real} <: AbstractMotionList{T} + motions::Vector{<:AbstractMotion{T}} end -MotionVector(motions...) = length([motions]) > 0 ? MotionVector([motions...]) : @error "You must provide at least one motion as input argument. If you do not want to define motion, use `NoMotion{T}()`" +MotionList(motions...) = length([motions]) > 0 ? MotionList([motions...]) : @error "You must provide at least one motion as input argument. If you do not want to define motion, use `NoMotion{T}()`" include("motion/SimpleMotion.jl") include("motion/ArbitraryMotion.jl") -Base.getindex(mv::MotionVector, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = MotionVector(getindex.(mv.motions, Ref(p))) -Base.view(mv::MotionVector, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = MotionVector(view.(mv.motions, Ref(p))) +Base.getindex(mv::MotionList, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = MotionList(getindex.(mv.motions, Ref(p))) +Base.view(mv::MotionList, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = MotionList(view.(mv.motions, Ref(p))) -""" Addition of MotionVectors """ -function Base.vcat(m1::MotionVector{T}, m2::MotionVector{T}, Ns1::Int, Ns2::Int) where {T<:Real} +""" Addition of MotionLists """ +function Base.vcat(m1::MotionList{T}, m2::MotionList{T}, Ns1::Int, Ns2::Int) where {T<:Real} mv1 = m1.motions - mv1_aux = Motion{T}[] + mv1_aux = AbstractMotion{T}[] for i in 1:length(mv1) if typeof(mv1[i]) <: ArbitraryMotion zeros1 = similar(mv1[i].dx, Ns2, size(mv1[i].dx, 2)) zeros1 .= zero(T) - push!(mv1_aux, typeof(mv1[i])(mv1[i].times, [[getfield(mv1[i], d); zeros1] for d in filter(x -> x != :times, fieldnames(typeof(mv1[i])))]...)) + push!(mv1_aux, typeof(mv1[i])(mv1[i].time, [[getfield(mv1[i], d); zeros1] for d in filter(x -> x != :time, fieldnames(typeof(mv1[i])))]...)) else push!(mv1_aux, mv1[i]) end end mv2 = m2.motions - mv2_aux = Motion{T}[] + mv2_aux = AbstractMotion{T}[] for i in 1:length(mv2) if typeof(mv2[i]) <: ArbitraryMotion zeros2 = similar(mv2[i].dx, Ns1, size(mv2[i].dx, 2)) zeros2 .= zero(T) - push!(mv2_aux, typeof(mv2[i])(mv2[i].times, [[zeros2; getfield(mv2[i], d)] for d in filter(x -> x != :times, fieldnames(typeof(mv2[i])))]...)) + push!(mv2_aux, typeof(mv2[i])(mv2[i].time, [[zeros2; getfield(mv2[i], d)] for d in filter(x -> x != :time, fieldnames(typeof(mv2[i])))]...)) else push!(mv2_aux, mv2[i]) end end - return MotionVector([mv1_aux; mv2_aux]) + return MotionList([mv1_aux; mv2_aux]) end """ Compare two motion vectors """ -function Base.:(==)(mv1::MotionVector{T}, mv2::MotionVector{T}) where {T<:Real} +function Base.:(==)(mv1::MotionList{T}, mv2::MotionList{T}) where {T<:Real} sort_motions!(mv1) sort_motions!(mv2) return reduce(&, mv1.motions .== mv2.motions) end -function Base.:(≈)(mv1::MotionVector{T}, mv2::MotionVector{T}) where {T<:Real} +function Base.:(≈)(mv1::MotionList{T}, mv2::MotionList{T}) where {T<:Real} sort_motions!(mv1) sort_motions!(mv2) return reduce(&, mv1.motions .≈ mv2.motions) @@ -60,7 +60,7 @@ Calculates the position of each spin at a set of arbitrary time instants, i.e. t For each dimension (x, y, z), the output matrix has ``N_{\text{spins}}`` rows and `length(t)` columns. # Arguments -- `motion`: (`::Vector{<:Motion{T<:Real}}`) phantom motion +- `motion`: (`::Vector{<:AbstractMotion{T<:Real}}`) 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 @@ -70,7 +70,7 @@ For each dimension (x, y, z), the output matrix has ``N_{\text{spins}}`` rows an - `x, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time """ function get_spin_coords( - mv::MotionVector{T}, + mv::MotionList{T}, x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, @@ -81,7 +81,7 @@ function get_spin_coords( # Buffers for displacements: ux, uy, uz = similar(xt), similar(yt), similar(zt) - # Composable motions: they need to be run sequentially + # Composable motions: they need to be run sequentially. Note that they depend on xt, yt , and zt for m in Iterators.filter(is_composable, mv.motions) displacement_x!(ux, m, xt, yt, zt, t) displacement_y!(uy, m, xt, yt, zt, t) @@ -105,16 +105,17 @@ end """ times = times(motion) """ -times(m::Motion) = times(m.times) -times(mv::MotionVector{T}) where {T<:Real} = begin +times(m::AbstractMotion) = times(m.time) +times(mv::MotionList{T}) where {T<:Real} = begin nodes = reduce(vcat, [times(m) for m in mv.motions]; init=[zero(T)]) return unique(sort(nodes)) end """ - sort_motions! + sort_motions!(motion_list) +sort_motions motions in a list according to their starting time """ -function sort_motions!(mv::MotionVector{T}) where {T<:Real} +function sort_motions!(mv::MotionList{T}) where {T<:Real} sort!(mv.motions; by=m -> times(m)[1]) return nothing end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/NoMotion.jl b/KomaMRIBase/src/datatypes/phantom/NoMotion.jl index 09f7774c8..ed52d186d 100644 --- a/KomaMRIBase/src/datatypes/phantom/NoMotion.jl +++ b/KomaMRIBase/src/datatypes/phantom/NoMotion.jl @@ -1,11 +1,11 @@ -struct NoMotion{T<:Real} <: AbstractMotion{T} end +struct NoMotion{T<:Real} <: AbstractMotionList{T} end Base.getindex(mv::NoMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = mv Base.view(mv::NoMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = mv """ Addition of NoMotions """ -Base.vcat(m1::NoMotion{T}, m2::AbstractMotion{T}, Ns1::Int, Ns2::Int) where {T<:Real} = m2 -Base.vcat(m1::AbstractMotion{T}, m2::NoMotion{T}, Ns1::Int, Ns2::Int) where {T<:Real} = m1 +Base.vcat(m1::NoMotion{T}, m2::AbstractMotionList{T}, Ns1::Int, Ns2::Int) where {T<:Real} = m2 +Base.vcat(m1::AbstractMotionList{T}, m2::NoMotion{T}, Ns1::Int, Ns2::Int) where {T<:Real} = m1 Base.:(==)(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = true Base.:(≈)(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = true diff --git a/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl index 8b6a2bb53..5a41f425e 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl @@ -19,13 +19,13 @@ const Interpolator2D = Interpolations.GriddedInterpolation{ """ ArbitraryMotion """ -abstract type ArbitraryMotion{T<:Real} <: Motion{T} end +abstract type ArbitraryMotion{T<:Real} <: AbstractMotion{T} end function Base.getindex(motion::ArbitraryMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) - return typeof(motion)(motion.times, [getfield(motion, d)[p,:] for d in filter(x -> x != :times, fieldnames(typeof(motion)))]...) + return typeof(motion)(motion.time, [getfield(motion, d)[p,:] for d in filter(x -> x != :time, fieldnames(typeof(motion)))]...) end function Base.view(motion::ArbitraryMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) - return typeof(motion)(motion.times, [@view(getfield(motion, d)[p,:]) for d in filter(x -> x != :times, fieldnames(typeof(motion)))]...) + return typeof(motion)(motion.time, [@view(getfield(motion, d)[p,:]) for d in filter(x -> x != :time, fieldnames(typeof(motion)))]...) end Base.:(==)(m1::ArbitraryMotion, m2::ArbitraryMotion) = (typeof(m1) == typeof(m2)) & reduce(&, [getfield(m1, field) == getfield(m2, field) for field in fieldnames(typeof(m1))]) @@ -68,7 +68,7 @@ function displacement_x!( t::AbstractArray{T}, ) where {T<:Real} itp = interpolate(motion.dx, Gridded(Linear()), Val(size(x,1))) - ux .= resample(itp, unit_time(t, motion.times)) + ux .= resample(itp, unit_time(t, motion.time)) return nothing end @@ -81,7 +81,7 @@ function displacement_y!( t::AbstractArray{T}, ) where {T<:Real} itp = interpolate(motion.dy, Gridded(Linear()), Val(size(x,1))) - uy .= resample(itp, unit_time(t, motion.times)) + uy .= resample(itp, unit_time(t, motion.time)) return nothing end @@ -94,7 +94,7 @@ function displacement_z!( t::AbstractArray{T}, ) where {T<:Real} itp = interpolate(motion.dz, Gridded(Linear()), Val(size(x,1))) - uz .= resample(itp, unit_time(t, motion.times)) + uz .= resample(itp, unit_time(t, motion.time)) return nothing end diff --git a/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl index bebd4ad06..f431d19a0 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl @@ -21,7 +21,7 @@ julia> motion = SimpleMotion( ) ``` """ -abstract type SimpleMotion{T<:Real} <: Motion{T} end +abstract type SimpleMotion{T<:Real} <: AbstractMotion{T} end Base.getindex(motion::SimpleMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = motion Base.view(motion::SimpleMotion, p::Union{AbstractRange, AbstractVector, Colon, Integer}) = motion diff --git a/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/FlowTrajectory.jl b/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/FlowTrajectory.jl index 56d08fb3a..c459625a4 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/FlowTrajectory.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/FlowTrajectory.jl @@ -1,7 +1,7 @@ -struct FlowTrajectory{T<:Real, TS<:TimeScale{T}} <: ArbitraryMotion{T} - times::TS +struct FlowTrajectory{T<:Real, TS<:AbstractTimeSpan{T}} <: ArbitraryMotion{T} + time::TS dx::AbstractArray{T} dy::AbstractArray{T} dz::AbstractArray{T} - resetmag::AbstractArray{Bool} + spin_reset::AbstractArray{Bool} end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/Trajectory.jl b/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/Trajectory.jl index 0aae0b094..a059ac618 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/Trajectory.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/arbitrarymotion/Trajectory.jl @@ -1,5 +1,5 @@ -struct Trajectory{T<:Real, TS<:TimeScale{T}} <: ArbitraryMotion{T} - times::TS +struct Trajectory{T<:Real, TS<:AbstractTimeSpan{T}} <: ArbitraryMotion{T} + time::TS dx::AbstractArray{T} dy::AbstractArray{T} dz::AbstractArray{T} diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl index 09a8146c4..ea1fdc211 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl @@ -1,11 +1,11 @@ @doc raw""" - heartbeat = HeartBeat(times, circumferential_strain, radial_strain, longitudinal_strain) + heartbeat = HeartBeat(time, circumferential_strain, radial_strain, longitudinal_strain) HeartBeat struct. It produces a heartbeat-like motion, characterised by three types of strain: Circumferential, Radial and Longitudinal # Arguments -- `times`: (`::TimeScale{T<:Real}`, `[s]`) time scale +- `time`: (`::AbstractTimeSpan{T<:Real}`, `[s]`) time scale - `circumferential_strain`: (`::Real`) contraction parameter - `radial_strain`: (`::Real`) contraction parameter - `longitudinal_strain`: (`::Real`) contraction parameter @@ -15,11 +15,11 @@ Circumferential, Radial and Longitudinal # Examples ```julia-repl -julia> hb = HeartBeat(times=Periodic(period=1.0, asymmetry=0.3), circumferential_strain=-0.3, radial_strain=-0.2, longitudinal_strain=0.0) +julia> hb = HeartBeat(time=Periodic(period=1.0, asymmetry=0.3), circumferential_strain=-0.3, radial_strain=-0.2, longitudinal_strain=0.0) ``` """ -@with_kw struct HeartBeat{T<:Real, TS<:TimeScale{T}} <: SimpleMotion{T} - times :: TS +@with_kw struct HeartBeat{T<:Real, TS<:AbstractTimeSpan{T}} <: SimpleMotion{T} + time :: TS circumferential_strain :: T radial_strain :: T longitudinal_strain :: T = typeof(circumferential_strain)(0.0) @@ -35,7 +35,7 @@ function displacement_x!( z::AbstractArray{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) r = sqrt.(x .^ 2 + y .^ 2) θ = atan.(y, x) Δ_circunferential = motion.circumferential_strain * maximum(r) @@ -57,7 +57,7 @@ function displacement_y!( z::AbstractArray{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) r = sqrt.(x .^ 2 + y .^ 2) θ = atan.(y, x) Δ_circunferential = motion.circumferential_strain * maximum(r) @@ -79,7 +79,7 @@ function displacement_z!( z::AbstractArray{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) uz .= t_unit .* (z .* motion.longitudinal_strain) return nothing 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 index 48628b5b8..e511f78d0 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl @@ -1,5 +1,5 @@ @doc raw""" - rotation = Rotation(times, pitch, roll, yaw) + rotation = Rotation(time, pitch, roll, yaw) Rotation motion struct. It produces a rotation of the phantom in the three axes: x (pitch), y (roll), and z (yaw). @@ -38,7 +38,7 @@ R &= R_z(\alpha) R_y(\beta) R_x(\gamma) \\ ``` # Arguments -- `times`: (`::TimeScale{T<:Real}`, `[s]`) time scale +- `time`: (`::AbstractTimeSpan{T<:Real}`, `[s]`) time scale - `pitch`: (`::Real`, `[º]`) rotation in x - `roll`: (`::Real`, `[º]`) rotation in y - `yaw`: (`::Real`, `[º]`) rotation in z @@ -48,16 +48,20 @@ R &= R_z(\alpha) R_y(\beta) R_x(\gamma) \\ # Examples ```julia-repl -julia> rt = Rotation(times=TimeRange(0.0, 0.5), pitch=15.0, roll=0.0, yaw=20.0) +julia> rt = Rotation(time=TimeRange(0.0, 0.5), pitch=15.0, roll=0.0, yaw=20.0) ``` """ -@with_kw struct Rotation{T<:Real, TS<:TimeScale{T}} <: SimpleMotion{T} - times :: TS +@with_kw struct Rotation{T<:Real, TS<:AbstractTimeSpan{T}} <: SimpleMotion{T} + time :: TS pitch :: T roll :: T yaw :: T end +RotationX(time::AbstractTimeSpan{T}, pitch::T) where {T<:Real} = Rotation(time, pitch, zero(T), zero(T)) +RotationY(time::AbstractTimeSpan{T}, roll::T) where {T<:Real} = Rotation(time, zero(T), roll, zero(T)) +RotationZ(time::AbstractTimeSpan{T}, yaw::T) where {T<:Real} = Rotation(time, zero(T), zero(T), yaw) + is_composable(motion::Rotation) = true function displacement_x!( @@ -68,7 +72,7 @@ function displacement_x!( z::AbstractArray{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) α = t_unit .* (motion.yaw) β = t_unit .* (motion.roll) γ = t_unit .* (motion.pitch) @@ -86,7 +90,7 @@ function displacement_y!( z::AbstractArray{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) α = t_unit .* (motion.yaw) β = t_unit .* (motion.roll) γ = t_unit .* (motion.pitch) @@ -104,7 +108,7 @@ function displacement_z!( z::AbstractArray{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) α = t_unit .* (motion.yaw) β = t_unit .* (motion.roll) γ = t_unit .* (motion.pitch) diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl index 88634f185..39949dda4 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl @@ -1,12 +1,12 @@ @doc raw""" - translation = Translation(times, dx, dy, dz) + translation = Translation(time, dx, dy, dz) Translation motion struct. It produces a linear translation of the phantom. Its fields are the final displacements in the three axes (dx, dy, dz) -and the start and end times of the translation. +and the start and end time of the translation. # Arguments -- `times`: (`::TimeScale{T<:Real}`, `[s]`) time scale +- `time`: (`::AbstractTimeSpan{T<:Real}`, `[s]`) time scale - `dx`: (`::Real`, `[m]`) translation in x - `dy`: (`::Real`, `[m]`) translation in y - `dz`: (`::Real`, `[m]`) translation in z @@ -16,16 +16,20 @@ and the start and end times of the translation. # Examples ```julia-repl -julia> tr = Translation(times=TimeRange(0.0, 0.5), dx=0.01, dy=0.02, dz=0.03) +julia> tr = Translation(time=TimeRange(0.0, 0.5), dx=0.01, dy=0.02, dz=0.03) ``` """ -@with_kw struct Translation{T<:Real, TS<:TimeScale{T}} <: SimpleMotion{T} - times :: TS +@with_kw struct Translation{T<:Real, TS<:AbstractTimeSpan{T}} <: SimpleMotion{T} + time :: TS dx :: T dy :: T dz :: T end +TranslationX(time::AbstractTimeSpan{T}, dx::T) where {T<:Real} = Translation(time, dx, zero(T), zero(T)) +TranslationY(time::AbstractTimeSpan{T}, dy::T) where {T<:Real} = Translation(time, zero(T), dy, zero(T)) +TranslationZ(time::AbstractTimeSpan{T}, dz::T) where {T<:Real} = Translation(time, zero(T), zero(T), dz) + function displacement_x!( ux::AbstractArray{T}, motion::Translation{T}, @@ -34,7 +38,7 @@ function displacement_x!( z::AbstractVector{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) ux .= t_unit .* motion.dx return nothing end @@ -47,7 +51,7 @@ function displacement_y!( z::AbstractVector{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) uy .= t_unit .* motion.dy return nothing end @@ -60,7 +64,7 @@ function displacement_z!( z::AbstractVector{T}, t::AbstractArray{T}, ) where {T<:Real} - t_unit = unit_time(t, motion.times) + t_unit = unit_time(t, motion.time) uz .= t_unit .* motion.dz return nothing end diff --git a/KomaMRIBase/src/timing/TimeScale.jl b/KomaMRIBase/src/timing/TimeScale.jl index 73adb1fc7..dad8167a0 100644 --- a/KomaMRIBase/src/timing/TimeScale.jl +++ b/KomaMRIBase/src/timing/TimeScale.jl @@ -1,12 +1,12 @@ -abstract type TimeScale{T<:Real} end +abstract type AbstractTimeSpan{T<:Real} end -@with_kw struct TimeRange{T<:Real} <: TimeScale{T} +@with_kw struct TimeRange{T<:Real} <: AbstractTimeSpan{T} t_start ::T t_end ::T = t_start @assert t_end >= t_start "t_end must be greater or equal than t_start" end -@with_kw struct Periodic{T<:Real} <: TimeScale{T} +@with_kw struct Periodic{T<:Real} <: AbstractTimeSpan{T} period::T asymmetry::T = eltype(period)(0.5) end @@ -46,8 +46,7 @@ function unit_time(t::AbstractArray{T}, ts::TimeRange{T}) where {T<:Real} if ts.t_start == ts.t_end return (t .>= ts.t_start) .* oneunit(T) else - tmp = max.((t .- ts.t_start) ./ (ts.t_end - ts.t_start), zero(T)) - return min.(tmp, oneunit(T)) + return min.(max.((t .- ts.t_start) ./ (ts.t_end - ts.t_start), zero(T)), oneunit(T)) end end diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl index 3b9b8aae8..9ec13403a 100644 --- a/KomaMRIBase/test/runtests.jl +++ b/KomaMRIBase/test/runtests.jl @@ -398,7 +398,7 @@ end t = collect(range(t_start, t_end, 11)) dx, dy, dz = [1.0, 0.0, 0.0] vx, vy, vz = [dx, dy, dz] ./ (t_end - t_start) - translation = MotionVector(Translation(TimeRange(t_start, t_end), dx, dy, dz)) + translation = MotionList(Translation(TimeRange(t_start, t_end), dx, dy, dz)) 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' @@ -406,7 +406,7 @@ end # ----- t_start = t_end -------- t_start = t_end = 0.0 t = [-0.5, -0.25, 0.0, 0.25, 0.5] - translation = MotionVector(Translation(TimeRange(t_start, t_end), dx, dy, dz)) + translation = MotionList(Translation(TimeRange(t_start, t_end), dx, dy, dz)) xt, yt, zt = get_spin_coords(translation, ph.x, ph.y, ph.z, t') @test xt == ph.x .+ dx*[0, 0, 1, 1, 1]' @test yt == ph.y .+ dy*[0, 0, 1, 1, 1]' @@ -420,7 +420,7 @@ end asymmetry = 0.5 dx, dy, dz = [1.0, 0.0, 0.0] vx, vy, vz = [dx, dy, dz] ./ (t_end - t_start) - periodictranslation = MotionVector(Translation(Periodic(period, asymmetry), dx, dy, dz)) + periodictranslation = MotionList(Translation(Periodic(period, asymmetry), dx, dy, dz)) 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' @@ -433,7 +433,7 @@ end pitch = 45.0 roll = 0.0 yaw = 45.0 - rotation = MotionVector(Rotation(TimeRange(t_start, t_end), pitch, roll, yaw)) + rotation = MotionList(Rotation(TimeRange(t_start, t_end), pitch, roll, yaw)) xt, yt, zt = get_spin_coords(rotation, ph.x, ph.y, ph.z, t') r = vcat(ph.x, ph.y, ph.z) R = rotz(π*yaw/180) * roty(π*roll/180) * rotx(π*pitch/180) @@ -444,7 +444,7 @@ end # ----- t_start = t_end -------- t_start = t_end = 0.0 t = [-0.5, -0.25, 0.0, 0.25, 0.5] - rotation = MotionVector(Rotation(TimeRange(t_start, t_end), pitch, roll, yaw)) + rotation = MotionList(Rotation(TimeRange(t_start, t_end), pitch, roll, yaw)) xt, yt, zt = get_spin_coords(rotation, ph.x, ph.y, ph.z, t') @test xt ≈ [ph.x ph.x rot_x rot_x rot_x] @test yt ≈ [ph.y ph.y rot_y rot_y rot_y] @@ -459,7 +459,7 @@ end pitch = 45.0 roll = 0.0 yaw = 45.0 - periodicrotation = MotionVector(Rotation(Periodic(period, asymmetry), pitch, roll, yaw)) + periodicrotation = MotionList(Rotation(Periodic(period, asymmetry), pitch, roll, yaw)) xt, yt, zt = get_spin_coords(periodicrotation, ph.x, ph.y, ph.z, t') r = vcat(ph.x, ph.y, ph.z) R = rotz(π*yaw/180) * roty(π*roll/180) * rotx(π*pitch/180) @@ -475,7 +475,7 @@ end circumferential_strain = -0.1 radial_strain = 0.0 longitudinal_strain = -0.1 - heartbeat = MotionVector(HeartBeat(TimeRange(t_start, t_end), circumferential_strain, radial_strain, longitudinal_strain)) + heartbeat = MotionList(HeartBeat(TimeRange(t_start, t_end), circumferential_strain, radial_strain, longitudinal_strain)) 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) @@ -485,7 +485,7 @@ end # ----- t_start = t_end -------- t_start = t_end = 0.0 t = [-0.5, -0.25, 0.0, 0.25, 0.5] - heartbeat = MotionVector(HeartBeat(TimeRange(t_start, t_end), circumferential_strain, radial_strain, longitudinal_strain)) + heartbeat = MotionList(HeartBeat(TimeRange(t_start, t_end), circumferential_strain, radial_strain, longitudinal_strain)) 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) @@ -505,7 +505,7 @@ end circumferential_strain = -0.1 radial_strain = 0.0 longitudinal_strain = -0.1 - periodicheartbeat = MotionVector(HeartBeat(Periodic(period, asymmetry), circumferential_strain, radial_strain, longitudinal_strain)) + periodicheartbeat = MotionList(HeartBeat(Periodic(period, asymmetry), circumferential_strain, radial_strain, longitudinal_strain)) xt, yt, zt = get_spin_coords(periodicheartbeat, ph.x, ph.y, ph.z, t') r = sqrt.(ph.x .^ 2 + ph.y .^ 2) θ = atan.(ph.y, ph.x) @@ -523,7 +523,7 @@ end dx = rand(Ns, Nt) dy = rand(Ns, Nt) dz = rand(Ns, Nt) - arbitrarymotion = MotionVector(Trajectory(TimeRange(t_start, t_end), dx, dy, dz)) + arbitrarymotion = MotionList(Trajectory(TimeRange(t_start, t_end), dx, dy, dz)) t = range(t_start, t_end, Nt) xt, yt, zt = get_spin_coords(arbitrarymotion, ph.x, ph.y, ph.z, t') @test xt == ph.x .+ dx @@ -538,7 +538,7 @@ end dx = rand(Ns, Nt) dy = rand(Ns, Nt) dz = rand(Ns, Nt) - arbitrarymotion = MotionVector(Trajectory(TimeRange(t_start, t_end), dx, dy, dz)) + arbitrarymotion = MotionList(Trajectory(TimeRange(t_start, t_end), dx, dy, dz)) t = range(t_start, t_end, Nt) xt, yt, zt = get_spin_coords(arbitrarymotion, ph.x, ph.y, ph.z, t') @test xt == ph.x .+ dx @@ -546,16 +546,16 @@ end @test zt == ph.z .+ dz end - simplemotion = MotionVector( - Translation(times=Periodic(period=0.5, asymmetry=0.5), dx=0.05, dy=0.05, dz=0.0), - Rotation(times=TimeRange(t_start=0.05, t_end=0.5), pitch=0.0, roll=0.0, yaw=π / 2) + simplemotion = MotionList( + Translation(time=Periodic(period=0.5, asymmetry=0.5), dx=0.05, dy=0.05, dz=0.0), + Rotation(time=TimeRange(t_start=0.05, t_end=0.5), pitch=0.0, roll=0.0, yaw=π / 2) ) Ns = length(obj1) Nt = 3 t_start = 0.0 t_end = 1.0 - arbitrarymotion = MotionVector(Trajectory(TimeRange(t_start, t_end), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt))) + arbitrarymotion = MotionList(Trajectory(TimeRange(t_start, t_end), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt), 0.01 .* rand(Ns, Nt))) # Test phantom subset obs1 = Phantom( diff --git a/KomaMRICore/src/simulation/Flow.jl b/KomaMRICore/src/simulation/Flow.jl index 6b74b8bfd..78af0d1df 100644 --- a/KomaMRICore/src/simulation/Flow.jl +++ b/KomaMRICore/src/simulation/Flow.jl @@ -5,20 +5,20 @@ function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{Complex{T}}, motion: return nothing end -function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{Complex{T}}, motion::MotionVector{T}, t::AbstractArray{T}) where {T<:Real} +function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{Complex{T}}, motion::MotionList{T}, t::AbstractArray{T}) where {T<:Real} for m in motion.motions reset_magnetization!(M, Mxy, m, t) end return nothing end -function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{Complex{T}}, motion::Motion{T}, t::AbstractArray{T}) where {T<:Real} +function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{Complex{T}}, motion::AbstractMotion{T}, t::AbstractArray{T}) where {T<:Real} return nothing end function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{Complex{T}}, motion::FlowTrajectory{T}, t::AbstractArray{T}) where {T<:Real} - itp = interpolate(motion.resetmag, Gridded(Constant{Previous}), Val(size(x,1))) - flags = resample(itp, unit_time(t, motion.times)) + itp = interpolate(motion.spin_reset, Gridded(Constant{Previous}), Val(size(x,1))) + flags = resample(itp, unit_time(t, motion.time)) reset = any(flags; dims=2) flags = .!(cumsum(flags; dims=2) .>= 1) Mxy .*= flags diff --git a/KomaMRICore/src/simulation/Functors.jl b/KomaMRICore/src/simulation/Functors.jl index 0354b763e..3217ab0a3 100644 --- a/KomaMRICore/src/simulation/Functors.jl +++ b/KomaMRICore/src/simulation/Functors.jl @@ -52,7 +52,7 @@ x = gpu(x, CUDABackend()) ``` """ gpu(x, backend::KA.GPU) = fmap(x -> adapt(backend, x), x; exclude=_isleaf) -adapt_storage(backend::KA.GPU, xs::MotionVector) = MotionVector(gpu.(xs.motions, Ref(backend))) +adapt_storage(backend::KA.GPU, xs::MotionList) = MotionList(gpu.(xs.motions, Ref(backend))) # To CPU """ @@ -76,7 +76,7 @@ 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::NoMotion) = NoMotion{T}() -adapt_storage(T::Type{<:Real}, xs::MotionVector) = MotionVector(paramtype.(T, xs.motions)) +adapt_storage(T::Type{<:Real}, xs::MotionList) = MotionList(paramtype.(T, xs.motions)) """ f32(m) diff --git a/KomaMRICore/test/runtests.jl b/KomaMRICore/test/runtests.jl index 69333e2c2..554cb6354 100644 --- a/KomaMRICore/test/runtests.jl +++ b/KomaMRICore/test/runtests.jl @@ -388,6 +388,49 @@ end @test NMRSE(sig, sig_jemris) < 1 #NMRSE < 1% end +@testitem "BlochSimple SimpleMotion" tags=[:important, :core, :motion] begin + using Suppressor + include("initialize_backend.jl") + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain_simple_motion() + + sim_params = Dict{String, Any}( + "gpu"=>USE_GPU, + "sim_method"=>KomaMRICore.BlochSimple(), + "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 "BlochSimple ArbitraryMotion" tags=[:important, :core, :motion] begin + using Suppressor + include("initialize_backend.jl") + include(joinpath(@__DIR__, "test_files", "utils.jl")) + + sig_jemris = signal_brain_motion_jemris() + seq = seq_epi_100x100_TE100_FOV230() + sys = Scanner() + obj = phantom_brain_arbitrary_motion() + + sim_params = Dict{String, Any}( + "gpu"=>USE_GPU, + "sim_method"=>KomaMRICore.BlochSimple(), + "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 "simulate_slice_profile" tags=[:core] begin using Suppressor include("initialize_backend.jl") diff --git a/KomaMRICore/test/test_files/utils.jl b/KomaMRICore/test/test_files/utils.jl index 7bc9d7caf..a194a1a16 100644 --- a/KomaMRICore/test/test_files/utils.jl +++ b/KomaMRICore/test/test_files/utils.jl @@ -18,7 +18,7 @@ end function phantom_brain_simple_motion() obj = phantom_brain() - obj.motion = MotionVector(Translation(times=TimeRange(0.0, 10.0), dx=0.0, dy=1.0, dz=0.0)) + obj.motion = MotionList(Translation(time=TimeRange(0.0, 10.0), dx=0.0, dy=1.0, dz=0.0)) return obj end @@ -30,7 +30,7 @@ function phantom_brain_arbitrary_motion() dx = zeros(Ns, 2) dz = zeros(Ns, 2) dy = [zeros(Ns,1) ones(Ns,1)] - obj.motion = MotionVector(Trajectory( + obj.motion = MotionList(Trajectory( TimeRange(t_start, t_end), dx, dy, diff --git a/KomaMRIFiles/src/Phantom/Phantom.jl b/KomaMRIFiles/src/Phantom/Phantom.jl index 51ea57af9..b695ed2ac 100644 --- a/KomaMRIFiles/src/Phantom/Phantom.jl +++ b/KomaMRIFiles/src/Phantom/Phantom.jl @@ -23,7 +23,7 @@ function read_phantom(filename::String) end end end - # Motion + # AbstractMotion motion_group = fid["motion"] import_motion!(phantom_fields, motion_group) @@ -56,10 +56,10 @@ end function import_motion!(phantom_fields::Array, motion_group::HDF5.Group) T = eltype(phantom_fields[2][2]) motion_type = read_attribute(motion_group, "type") - if motion_type == "MotionVector" - simple_motion_types = last.(split.(string.(reduce(vcat,(subtypes(subtypes(Motion)[2])))), ".")) - arbitrary_motion_types = last.(split.(string.(reduce(vcat,(subtypes(subtypes(Motion)[1])))), ".")) - motion_array = Motion{T}[] + if motion_type == "MotionList" + simple_motion_types = last.(split.(string.(reduce(vcat,(subtypes(subtypes(AbstractMotion)[2])))), ".")) + arbitrary_motion_types = last.(split.(string.(reduce(vcat,(subtypes(subtypes(AbstractMotion)[1])))), ".")) + motion_array = AbstractMotion{T}[] for key in keys(motion_group) type_group = motion_group[key] type_str = split(key, "_")[2] @@ -67,26 +67,26 @@ function import_motion!(phantom_fields::Array, motion_group::HDF5.Group) args = [] for smtype in subtypes(SimpleMotion) if type_str == last(split(string(smtype), ".")) - times = import_time_range(type_group["times"]) - type_fields = filter(x -> x != :times, fieldnames(smtype)) + time = import_time_range(type_group["time"]) + type_fields = filter(x -> x != :time, fieldnames(smtype)) for key in type_fields push!(args, read_attribute(type_group, string(key))) end - push!(motion_array, smtype(times, args...)) + push!(motion_array, smtype(time, args...)) end end for amtype in subtypes(ArbitraryMotion) if type_str == last(split(string(amtype), ".")) - times = import_time_range(type_group["times"]) - type_fields = filter(x -> x != :times, fieldnames(amtype)) + time = import_time_range(type_group["time"]) + type_fields = filter(x -> x != :time, fieldnames(amtype)) for key in type_fields push!(args, read(type_group[string(key)])) end - push!(motion_array, amtype(times, args...)) + push!(motion_array, amtype(time, args...)) end end end - return push!(phantom_fields, (:motion, MotionVector(motion_array))) + return push!(phantom_fields, (:motion, MotionList(motion_array))) elseif motion_type == "NoMotion" return push!(phantom_fields, (:motion, NoMotion{T}())) end @@ -94,7 +94,7 @@ end function import_time_range(times_group::HDF5.Group) time_scale_type = read_attribute(times_group, "type") - for tstype in subtypes(TimeScale) + for tstype in subtypes(AbstractTimeSpan) if time_scale_type == last(split(string(tstype), ".")) args = [] for key in filter(x -> x != :type, fieldnames(tstype)) @@ -147,13 +147,13 @@ function write_phantom( return close(fid) end -function export_motion!(motion_group::HDF5.Group, mv::MotionVector{T}) where {T<:Real} - HDF5.attributes(motion_group)["type"] = "MotionVector" +function export_motion!(motion_group::HDF5.Group, mv::MotionList{T}) where {T<:Real} + HDF5.attributes(motion_group)["type"] = "MotionList" for (counter, m) in enumerate(mv.motions) type_name = typeof(m).name.name type_group = create_group(motion_group, "$(counter)_$type_name") - export_time_range!(type_group, m.times) - type_fields = filter(x -> x != :times, fieldnames(typeof(m))) + export_time_range!(type_group, m.time) + type_fields = filter(x -> x != :time, fieldnames(typeof(m))) for field in type_fields field_value = getfield(m, field) if typeof(field_value) <: Number @@ -169,12 +169,12 @@ function export_motion!(motion_group::HDF5.Group, motion::NoMotion{T}) where {T< HDF5.attributes(motion_group)["type"] = "NoMotion" end -function export_time_range!(type_group::HDF5.Group, times::TimeScale) - times_name = typeof(times).name.name - times_group = create_group(type_group, "times") +function export_time_range!(type_group::HDF5.Group, time::AbstractTimeSpan) + times_name = typeof(time).name.name + times_group = create_group(type_group, "time") HDF5.attributes(times_group)["type"] = string(times_name) - for field in fieldnames(typeof(times)) - field_value = getfield(times, field) + for field in fieldnames(typeof(time)) + field_value = getfield(time, field) HDF5.attributes(times_group)[string(field)] = field_value end end \ No newline at end of file diff --git a/KomaMRIFiles/test/runtests.jl b/KomaMRIFiles/test/runtests.jl index dad477288..b9dbda971 100644 --- a/KomaMRIFiles/test/runtests.jl +++ b/KomaMRIFiles/test/runtests.jl @@ -64,14 +64,14 @@ using TestItems, TestItemRunner path = @__DIR__ filename = path * "/test_files/brain_simplemotion_w.phantom" obj1 = brain_phantom2D() - obj1.motion = MotionVector( + obj1.motion = MotionList( Rotation( - times=Periodic(period=1.0), + time=Periodic(period=1.0), yaw=45.0, pitch=0.0, roll=0.0), Translation( - times=TimeRange(t_start=0.0, t_end=0.5), + time=TimeRange(t_start=0.0, t_end=0.5), dx=0.0, dy=0.02, dz=0.0 @@ -90,7 +90,7 @@ using TestItems, TestItemRunner K = 10 t_start = 0.0 t_end = 1.0 - obj1.motion = MotionVector(Trajectory( + obj1.motion = MotionList(Trajectory( TimeRange(t_start, t_end), 0.01.*rand(Ns, K-1), 0.01.*rand(Ns, K-1), diff --git a/examples/3.tutorials/lit-05-SimpleMotion.jl b/examples/3.tutorials/lit-05-SimpleMotion.jl index 93b3da121..22fc9659d 100644 --- a/examples/3.tutorials/lit-05-SimpleMotion.jl +++ b/examples/3.tutorials/lit-05-SimpleMotion.jl @@ -19,8 +19,8 @@ obj.Δw .= 0 # hide # # In this example, we will add a [`Translation`](@ref) of 2 cm in x, with duration of 200 ms (v = 0.1 m/s): -obj.motion = SimpleMotion( - Translation(t_start=0.0, t_end=200e-3, dx=2e-2, dy=0.0, dz=0.0) +obj.motion = MotionList( + Translation(time=TimeRange(t_start=0.0, t_end=200e-3), dx=2e-2, dy=0.0, dz=0.0) ) p1 = plot_phantom_map(obj, :T2 ; height=450, intermediate_time_samples=4) # hide