diff --git a/KomaMRIBase/src/KomaMRIBase.jl b/KomaMRIBase/src/KomaMRIBase.jl index 3fd03d493..2e15501e3 100644 --- a/KomaMRIBase/src/KomaMRIBase.jl +++ b/KomaMRIBase/src/KomaMRIBase.jl @@ -31,6 +31,7 @@ include("datatypes/Phantom.jl") include("datatypes/simulation/DiscreteSequence.jl") include("timing/TimeStepCalculation.jl") include("timing/TrapezoidalIntegration.jl") +include("timing/UnitTime.jl") # Main export γ # gyro-magnetic ratio [Hz/T] @@ -51,8 +52,7 @@ export NoMotion, SimpleMotion, ArbitraryMotion export SimpleMotionType export Translation, Rotation, HeartBeat export PeriodicTranslation, PeriodicRotation, PeriodicHeartBeat -export get_spin_coords, sort_motions! -export LinearInterpolator +export get_spin_coords # Secondary export get_kspace, rotx, roty, rotz # Additionals diff --git a/KomaMRIBase/src/datatypes/Phantom.jl b/KomaMRIBase/src/datatypes/Phantom.jl index 08e33acfb..8223084f2 100644 --- a/KomaMRIBase/src/datatypes/Phantom.jl +++ b/KomaMRIBase/src/datatypes/Phantom.jl @@ -69,10 +69,7 @@ Base.getindex(x::Phantom, i::Integer) = x[i:i] """Compare two phantoms""" Base.:(==)(obj1::Phantom, obj2::Phantom) = reduce( &, - [ - getfield(obj1, field) == getfield(obj2, field) for - field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom)) - ], + [getfield(obj1, field) == getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))], ) Base.:(≈)(obj1::Phantom, obj2::Phantom) = reduce(&, [getfield(obj1, field) ≈ getfield(obj2, field) for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))]) Base.:(==)(m1::MotionModel, m2::MotionModel) = false @@ -88,7 +85,13 @@ Base.getindex(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begi end """Separate object spins in a sub-group (lightweigth).""" -Base.view(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = @views obj[p] +Base.view(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begin + fields = [] + for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom)) + push!(fields, (field, @view(getfield(obj, field)[p]))) + end + return Phantom(; name=obj.name, fields...) +end """Addition of phantoms""" +(obj1::Phantom, obj2::Phantom) = begin @@ -117,10 +120,6 @@ function get_dims(obj::Phantom) return dims end -function sort_motions!(motion::MotionModel) - return nothing -end - """ obj = heart_phantom(...) @@ -178,7 +177,7 @@ function heart_phantom( Dλ1=Dλ1[ρ .!= 0], Dλ2=Dλ2[ρ .!= 0], Dθ=Dθ[ρ .!= 0], - motion=SimpleMotion([ + motion=SimpleMotion( PeriodicHeartBeat(; period=period, asymmetry=asymmetry, @@ -189,7 +188,7 @@ function heart_phantom( PeriodicRotation(; period=period, asymmetry=asymmetry, yaw=rotation_angle, pitch=0.0, roll=0.0 ), - ]), + ), ) return phantom end diff --git a/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl index e3f3256dc..5191c8b43 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl @@ -2,13 +2,24 @@ # Interpolator{T,Degree,ETPType}, # Degree = Linear,Cubic.... # ETPType = Periodic, Flat... -const LinearInterpolator = Interpolations.Extrapolation{ - T, - 1, - Interpolations.GriddedInterpolation{T,1,V,Gridded{Linear{Throw{OnGrid}}},Tuple{V}}, - Gridded{Linear{Throw{OnGrid}}}, - Periodic{Nothing}, -} where {T<:Real,V<:AbstractVector{T}} + +const Interpolator1D = Interpolations.GriddedInterpolation{ + T,1,V,Itp,K +} where { + T<:Real, + V<:AbstractArray{T}, + Itp<:Interpolations.Gridded{Linear{Throw{OnGrid}}}, + K<:Tuple{AbstractVector{T}}, +} + +const Interpolator2D = Interpolations.GriddedInterpolation{ + T,2,V,Itp,K +} where { + T<:Real, + V<:AbstractArray{T}, + Itp<:Interpolations.Gridded{Linear{Throw{OnGrid}}}, + K<:Tuple{AbstractVector{T}, AbstractVector{T}}, +} """ motion = ArbitraryMotion(period_durations, dx, dy, dz) @@ -43,106 +54,82 @@ julia> motion = ArbitraryMotion( ) ``` """ -struct ArbitraryMotion{T<:Real,V<:AbstractVector{T}} <: MotionModel{T} - period_durations::Vector{T} - dx::Array{T,2} - dy::Array{T,2} - dz::Array{T,2} - ux::Vector{LinearInterpolator{T,V}} - uy::Vector{LinearInterpolator{T,V}} - uz::Vector{LinearInterpolator{T,V}} -end - -function ArbitraryMotion( - period_durations::AbstractVector{T}, - dx::AbstractArray{T,2}, - dy::AbstractArray{T,2}, - dz::AbstractArray{T,2}, -) where {T<:Real} - @warn "Note that ArbitraryMotion is under development so it is not optimized so far" maxlog = 1 - Ns = size(dx)[1] - num_pieces = size(dx)[2] + 1 - limits = times(period_durations, num_pieces) - - #! format: off - Δ = zeros(Ns,length(limits),4) - Δ[:,:,1] = hcat(repeat(hcat(zeros(Ns,1),dx),1,length(period_durations)),zeros(Ns,1)) - Δ[:,:,2] = hcat(repeat(hcat(zeros(Ns,1),dy),1,length(period_durations)),zeros(Ns,1)) - Δ[:,:,3] = hcat(repeat(hcat(zeros(Ns,1),dz),1,length(period_durations)),zeros(Ns,1)) - - etpx = [extrapolate(interpolate((limits,), Δ[i,:,1], Gridded(Linear())), Periodic()) for i in 1:Ns] - etpy = [extrapolate(interpolate((limits,), Δ[i,:,2], Gridded(Linear())), Periodic()) for i in 1:Ns] - etpz = [extrapolate(interpolate((limits,), Δ[i,:,3], Gridded(Linear())), Periodic()) for i in 1:Ns] - #! format: on - - return ArbitraryMotion(period_durations, dx, dy, dz, etpx, etpy, etpz) +struct ArbitraryMotion{T} <: MotionModel{T} + t_start::T + t_end::T + dx::AbstractArray{T} + dy::AbstractArray{T} + dz::AbstractArray{T} end function Base.getindex( motion::ArbitraryMotion, p::Union{AbstractRange,AbstractVector,Colon} ) - fields = [] - for field in fieldnames(ArbitraryMotion) - if field in (:dx, :dy, :dz) - push!(fields, getfield(motion, field)[p, :]) - elseif field in (:ux, :uy, :uz) - push!(fields, getfield(motion, field)[p]) - else - push!(fields, getfield(motion, field)) - end - end - return ArbitraryMotion(fields...) + return ArbitraryMotion(motion.t_start, motion.t_end, motion.dx[p,:], motion.dy[p,:], motion.dz[p,:]) +end +function Base.view( + motion::ArbitraryMotion, p::Union{AbstractRange,AbstractVector,Colon} +) + return ArbitraryMotion(motion.t_start, motion.t_end, @view(motion.dx[p,:]), @view(motion.dy[p,:]), @view(motion.dz[p,:])) end Base.:(==)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) == getfield(m2, field) for field in fieldnames(ArbitraryMotion)]) Base.:(≈)(m1::ArbitraryMotion, m2::ArbitraryMotion) = reduce(&, [getfield(m1, field) ≈ getfield(m2, field) for field in fieldnames(ArbitraryMotion)]) function Base.vcat(m1::ArbitraryMotion, m2::ArbitraryMotion) - fields = [] - @assert m1.period_durations == m2.period_durations "period_durations of both ArbitraryMotions must be the same" - for field in - Iterators.filter(x -> !(x == :period_durations), fieldnames(ArbitraryMotion)) - push!(fields, [getfield(m1, field); getfield(m2, field)]) - end - return ArbitraryMotion(m1.period_durations, fields...) + @assert (m1.t_start == m2.t_start) && (m1.t_end == m2.t_end) "Starting and ending times must be the same" + return ArbitraryMotion(m1.t_start, m1.t_end, [m1.dx; m2.dx], [m1.dy; m2.dy], [m1.dz; m2.dz]) end """ limits = times(obj.motion) """ function times(motion::ArbitraryMotion) - period_durations = motion.period_durations - num_pieces = size(motion.dx)[2] + 1 - return times(period_durations, num_pieces) + return range(motion.t_start, motion.t_end, length=size(motion.dx, 2)) +end + +function GriddedInterpolation(nodes, A, ITP) + return Interpolations.GriddedInterpolation{eltype(A), length(nodes), typeof(A), typeof(ITP), typeof(nodes)}(nodes, A, ITP) +end + +function interpolate(motion::ArbitraryMotion{T}, Ns::Val{1}) where {T<:Real} + _, Nt = size(motion.dx) + t = similar(motion.dx, Nt); copyto!(t, collect(range(zero(T), oneunit(T), Nt))) + itpx = GriddedInterpolation((t, ), motion.dx[:], Gridded(Linear())) + itpy = GriddedInterpolation((t, ), motion.dy[:], Gridded(Linear())) + itpz = GriddedInterpolation((t, ), motion.dz[:], Gridded(Linear())) + return itpx, itpy, itpz +end + +function interpolate(motion::ArbitraryMotion{T}, Ns::Val) where {T<:Real} + Ns, Nt = size(motion.dx) + id = similar(motion.dx, Ns); copyto!(id, collect(range(oneunit(T), T(Ns), Ns))) + t = similar(motion.dx, Nt); copyto!(t, collect(range(zero(T), oneunit(T), Nt))) + itpx = GriddedInterpolation((id, t), motion.dx, Gridded(Linear())) + itpy = GriddedInterpolation((id, t), motion.dy, Gridded(Linear())) + itpz = GriddedInterpolation((id, t), motion.dz, Gridded(Linear())) + return itpx, itpy, itpz +end + +function resample(itpx::Interpolator1D{T}, itpy::Interpolator1D{T}, itpz::Interpolator1D{T}, t::AbstractArray{T}) where {T<:Real} + return itpx.(t), itpy.(t), itpz.(t) end -function times(period_durations::AbstractVector, num_pieces::Int) - # Pre-allocating memory - limits = zeros(eltype(period_durations), num_pieces * length(period_durations) + 1) - - idx = 1 - for i in 1:length(period_durations) - segment_increment = period_durations[i] / num_pieces - cumulative_sum = limits[idx] # Start from the last computed value in limits - for j in 1:num_pieces - cumulative_sum += segment_increment - limits[idx + 1] = cumulative_sum - idx += 1 - end - end - return limits +function resample(itpx::Interpolator2D{T}, itpy::Interpolator2D{T}, itpz::Interpolator2D{T}, t::AbstractArray{T}) where {T<:Real} + Ns = size(itpx.coefs, 1) + id = similar(itpx.coefs, Ns) + copyto!(id, collect(range(oneunit(T), T(Ns), Ns))) + return itpx.(id, t), itpy.(id, t), itpz.(id, t) end -# TODO: Calculate interpolation functions "on the fly" function get_spin_coords( motion::ArbitraryMotion{T}, x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, - t::AbstractArray{T}, + t::AbstractArray{T} ) where {T<:Real} - xt = x .+ reduce(vcat, [etp.(t) for etp in motion.ux]) - yt = y .+ reduce(vcat, [etp.(t) for etp in motion.uy]) - zt = z .+ reduce(vcat, [etp.(t) for etp in motion.uz]) - return xt, yt, zt -end + motion_functions = interpolate(motion, Val(size(x,1))) + ux, uy, uz = resample(motion_functions..., unit_time(t, motion.t_start, motion.t_end)) + return x .+ ux, y .+ uy, z .+ uz +end \ No newline at end of file diff --git a/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl index 3a453d98c..4a6895f76 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/NoMotion.jl @@ -7,6 +7,7 @@ x = x struct NoMotion{T<:Real} <: MotionModel{T} end Base.getindex(motion::NoMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion +Base.view(motion::NoMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion Base.:(==)(m1::NoMotion, m2::NoMotion) = true Base.:(≈)(m1::NoMotion, m2::NoMotion) = true @@ -18,7 +19,7 @@ function get_spin_coords( x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, - t::AbstractArray{T}, + t::AbstractArray{T} ) where {T<:Real} return x, y, z end diff --git a/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl b/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl index 166832509..75d0b0576 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl @@ -8,38 +8,34 @@ is_composable(motion_type::SimpleMotionType{T}) where {T<:Real} = false SimpleMotion model. It allows for the definition of motion by means of simple parameters. The `SimpleMotion` struct is composed by only one field, called `types`, -which is a vector of simple motion types. This vector will contain as many elements +which is a tuple of simple motion types. This tuple will contain as many elements as simple motions we want to combine. # Arguments -- `types`: (`::Vector{<:SimpleMotionType{T}}`) vector of simple motion types +- `types`: (`::Tuple{Vararg{<:SimpleMotionType{T}}}`) vector of simple motion types # Returns - `motion`: (`::SimpleMotion`) SimpleMotion struct # Examples ```julia-repl -julia> motion = SimpleMotion([ +julia> motion = SimpleMotion( Translation(dx=0.01, dy=0.02, dz=0.0, t_start=0.0, t_end=0.5), Rotation(pitch=15.0, roll=0.0, yaw=20.0, t_start=0.1, t_end=0.5), HeartBeat(circumferential_strain=-0.3, radial_strain=-0.2, longitudinal_strain=0.0, t_start=0.2, t_end=0.5) - ]) + ) ``` """ struct SimpleMotion{T<:Real} <: MotionModel{T} - types::Vector{<:SimpleMotionType{T}} + types::Tuple{Vararg{<:SimpleMotionType{T}}} end +SimpleMotion(types...) = SimpleMotion(types) + Base.getindex(motion::SimpleMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion -function Base.getindex( - motion::SimpleMotion, - p::Union{AbstractRange,AbstractVector,Colon}, - q::Union{AbstractRange,AbstractVector,Colon}, -) - return motion -end +Base.view(motion::SimpleMotion, p::Union{AbstractRange,AbstractVector,Colon}) = motion -Base.vcat(m1::SimpleMotion, m2::SimpleMotion) = SimpleMotion([m1.types; m2.types]) +Base.vcat(m1::SimpleMotion, m2::SimpleMotion) = SimpleMotion(m1.types..., m2.types...) Base.:(==)(m1::SimpleMotion, m2::SimpleMotion) = reduce(&, [m1.types[i] == m2.types[i] for i in 1:length(m1.types)]) Base.:(≈)(m1::SimpleMotion, m2::SimpleMotion) = reduce(&, [m1.types[i] ≈ m2.types[i] for i in 1:length(m1.types)]) @@ -71,21 +67,31 @@ function get_spin_coords( x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, - t::AbstractArray{T}, + t::AbstractArray{T} ) where {T<:Real} + motion = sort_motion(motion) + # Buffers for positions: xt, yt, zt = x .+ 0*t, y .+ 0*t, z .+ 0*t + # Buffers for displacements: + ux, uy, uz = similar(xt), similar(yt), similar(zt) + # Composable motions: they need to be run sequentially for motion in Iterators.filter(is_composable, motion.types) - aux = xt .+ 0, yt .+ 0, zt .+ 0 - xt .+= displacement_x(motion, aux..., t) - yt .+= displacement_y(motion, aux..., t) - zt .+= displacement_z(motion, aux..., t) + displacement_x!(ux, motion, xt, yt, zt, t) + displacement_y!(uy, motion, xt, yt, zt, t) + displacement_z!(uz, motion, xt, yt, zt, t) + xt .+= ux + yt .+= uy + zt .+= uz end # Additive motions: these motions can be run in parallel for motion in Iterators.filter(!is_composable, motion.types) - xt .+= displacement_x(motion, x, y, z, t) - yt .+= displacement_y(motion, x, y, z, t) - zt .+= displacement_z(motion, x, y, z, t) + displacement_x!(ux, motion, x, y, z, t) + displacement_y!(uy, motion, x, y, z, t) + displacement_z!(uz, motion, x, y, z, t) + xt .+= ux + yt .+= uy + zt .+= uz end return xt, yt, zt end @@ -96,38 +102,16 @@ function times(motion::SimpleMotion) return nodes end -function sort_motions!(motion::SimpleMotion) - return sort!(motion.types; by=m -> times(m)[1]) -end +# Sort Motion +include("simplemotion/sorting.jl") -# --------- Simple Motion Types: ------------- +# Simple Motion Types: # Non-periodic types: defined by an initial time (t_start), an end time (t_end) and a displacement include("simplemotion/Translation.jl") include("simplemotion/Rotation.jl") include("simplemotion/HeartBeat.jl") - -function unit_time(t::AbstractArray{T}, t_start::T, t_end::T) where {T<:Real} - if t_start == t_end - return (t .>= t_start) .* one(T) # The problem with this is that it returns a BitVector (type stability issues) - else - return min.(max.((t .- t_start) ./ (t_end - t_start), zero(T)), one(T)) - end -end - # Periodic types: defined by the period, the temporal symmetry and a displacement (amplitude) include("simplemotion/PeriodicTranslation.jl") include("simplemotion/PeriodicRotation.jl") include("simplemotion/PeriodicHeartBeat.jl") -function unit_time_triangular(t, period, asymmetry) - t_rise = period * asymmetry - t_fall = period * (1 - asymmetry) - t_relative = mod.(t, period) - t_unit = - ifelse.( - t_relative .< t_rise, - t_relative ./ t_rise, - 1 .- (t_relative .- t_rise) ./ t_fall, - ) - return t_unit -end diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl index 4db9f8a7c..c2b8d5181 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl @@ -1,5 +1,5 @@ @doc raw""" - heartbeat = HeartBeat(t_start, t_end, dx, dy, dz) + heartbeat = HeartBeat(circumferential_strain, radial_strain, longitudinal_strain, t_start, t_end) HeartBeat struct. It produces a heartbeat-like motion, characterised by three types of strain: Circumferential, Radial and Longitudinal @@ -25,12 +25,13 @@ julia> hb = HeartBeat(circumferential_strain=-0.3, radial_strain=-0.2, longitudi longitudinal_strain::T = typeof(circumferential_strain)(0.0) t_start::T = typeof(circumferential_strain)(0.0) t_end::T = typeof(circumferential_strain)(0.0) - @assert t_end >= t_start "t_end must be major or equal than t_start" + @assert t_end >= t_start "t_end must be greater or equal than t_start" end is_composable(motion_type::HeartBeat) = true -function displacement_x( +function displacement_x!( + ux::AbstractArray{T}, motion_type::HeartBeat{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -47,10 +48,12 @@ function displacement_x( neg = (r .+ Δr) .< 0 Δr = (.!neg) .* Δr Δr .-= neg .* r - return Δr .* cos.(θ) + ux .= Δr .* cos.(θ) + return nothing end -function displacement_y( +function displacement_y!( + uy::AbstractArray{T}, motion_type::HeartBeat{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -67,10 +70,12 @@ function displacement_y( neg = (r .+ Δr) .< 0 Δr = (.!neg) .* Δr Δr .-= neg .* r - return Δr .* sin.(θ) + uy .= Δr .* sin.(θ) + return nothing end -function displacement_z( +function displacement_z!( + uz::AbstractArray{T}, motion_type::HeartBeat{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -78,7 +83,8 @@ function displacement_z( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) - return t_unit .* (z .* motion_type.longitudinal_strain) + uz .= t_unit .* (z .* motion_type.longitudinal_strain) + return nothing end times(motion_type::HeartBeat) = begin diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl index d26a3f394..2a9c02537 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl @@ -1,5 +1,5 @@ @doc raw""" - periodic_heartbeat = PeriodicHeartBeat(t_start, t_end, dx, dy, dz) + periodic_heartbeat = PeriodicHeartBeat(circumferential_strain, radial_strain, longitudinal_strain, period, asymmetry) HeartBeat struct. It produces a heartbeat-like motion, characterised by three types of strain: Circumferential, Radial and Longitudinal @@ -24,7 +24,8 @@ end is_composable(motion_type::PeriodicHeartBeat) = true -function displacement_x( +function displacement_x!( + ux::AbstractArray{T}, motion_type::PeriodicHeartBeat{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -41,10 +42,12 @@ function displacement_x( neg = (r .+ Δr) .< 0 Δr = (.!neg) .* Δr Δr .-= neg .* r - return Δr .* cos.(θ) + ux .= Δr .* cos.(θ) + return nothing end -function displacement_y( +function displacement_y!( + uy::AbstractArray{T}, motion_type::PeriodicHeartBeat{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -61,10 +64,12 @@ function displacement_y( neg = (r .+ Δr) .< 0 Δr = (.!neg) .* Δr Δr .-= neg .* r - return Δr .* sin.(θ) + uy .= Δr .* sin.(θ) + return nothing end -function displacement_z( +function displacement_z!( + uz::AbstractArray{T}, motion_type::PeriodicHeartBeat{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -72,7 +77,8 @@ function displacement_z( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) - return t_unit .* (z .* motion_type.longitudinal_strain) + uz .= t_unit .* (z .* motion_type.longitudinal_strain) + return nothing end function times(motion_type::PeriodicHeartBeat) diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl index e2ff92e41..8bb8cc5f1 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl @@ -1,5 +1,5 @@ @doc raw""" - periodic_rotation = PeriodicRotation(period, asymmetry, pitch, roll, yaw) + periodic_rotation = PeriodicRotation(pitch, roll, yaw, period, asymmetry) PeriodicRotation motion struct. It produces a rotation of the phantom in the three axes: x (pitch), y (roll), and z (yaw) @@ -25,7 +25,8 @@ end is_composable(motion_type::PeriodicRotation) = true -function displacement_x( +function displacement_x!( + ux::AbstractArray{T}, motion_type::PeriodicRotation{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -36,12 +37,14 @@ function displacement_x( α = t_unit .* motion_type.pitch β = t_unit .* motion_type.roll γ = t_unit .* motion_type.yaw - return cosd.(γ) .* cosd.(β) .* x + - (cosd.(γ) .* sind.(β) .* sind.(α) .- sind.(γ) .* cosd.(α)) .* y + - (cosd.(γ) .* sind.(β) .* cosd.(α) .+ sind.(γ) .* sind.(α)) .* z .- x + ux .= cosd.(γ) .* cosd.(β) .* x + + (cosd.(γ) .* sind.(β) .* sind.(α) .- sind.(γ) .* cosd.(α)) .* y + + (cosd.(γ) .* sind.(β) .* cosd.(α) .+ sind.(γ) .* sind.(α)) .* z .- x + return nothing end -function displacement_y( +function displacement_y!( + uy::AbstractArray{T}, motion_type::PeriodicRotation{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -52,12 +55,14 @@ function displacement_y( α = t_unit .* motion_type.pitch β = t_unit .* motion_type.roll γ = t_unit .* motion_type.yaw - return sind.(γ) .* cosd.(β) .* x + - (sind.(γ) .* sind.(β) .* sind.(α) .+ cosd.(γ) .* cosd.(α)) .* y + - (sind.(γ) .* sind.(β) .* cosd.(α) .- cosd.(γ) .* sind.(α)) .* z .- y + uy .= sind.(γ) .* cosd.(β) .* x + + (sind.(γ) .* sind.(β) .* sind.(α) .+ cosd.(γ) .* cosd.(α)) .* y + + (sind.(γ) .* sind.(β) .* cosd.(α) .- cosd.(γ) .* sind.(α)) .* z .- y + return nothing end -function displacement_z( +function displacement_z!( + uz::AbstractArray{T}, motion_type::PeriodicRotation{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -68,9 +73,10 @@ function displacement_z( α = t_unit .* motion_type.pitch β = t_unit .* motion_type.roll γ = t_unit .* motion_type.yaw - return -sind.(β) .* x + - cosd.(β) .* sind.(α) .* y + - cosd.(β) .* cosd.(α) .* z .- z + uz .= -sind.(β) .* x + + cosd.(β) .* sind.(α) .* y + + cosd.(β) .* cosd.(α) .* z .- z + return nothing end function times(motion_type::PeriodicRotation) diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl index 4347ad33d..7dcebaac8 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicTranslation.jl @@ -1,5 +1,5 @@ @doc raw""" - periodic_translation = PeriodicTranslation(period, asymmetry, dx, dy, dz) + periodic_translation = PeriodicTranslation(dx, dy, dz, period, asymmetry) PeriodicTranslation motion struct. It produces a periodic translation of the phantom in the three directions x, y and z. The amplitude of the oscillation will be defined by dx, dy and dz @@ -23,7 +23,8 @@ The amplitude of the oscillation will be defined by dx, dy and dz asymmetry::T = typeof(dx)(0.5) end -function displacement_x( +function displacement_x!( + ux::AbstractArray{T}, motion_type::PeriodicTranslation{T}, x::AbstractVector{T}, y::AbstractVector{T}, @@ -31,10 +32,12 @@ function displacement_x( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) - return t_unit .* motion_type.dx + ux .= t_unit .* motion_type.dx + return nothing end -function displacement_y( +function displacement_y!( + uy::AbstractArray{T}, motion_type::PeriodicTranslation{T}, x::AbstractVector{T}, y::AbstractVector{T}, @@ -42,10 +45,12 @@ function displacement_y( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) - return t_unit .* motion_type.dy + uy .= t_unit .* motion_type.dy + return nothing end -function displacement_z( +function displacement_z!( + uz::AbstractArray{T}, motion_type::PeriodicTranslation{T}, x::AbstractVector{T}, y::AbstractVector{T}, @@ -53,7 +58,8 @@ function displacement_z( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time_triangular(t, motion_type.period, motion_type.asymmetry) - return t_unit .* motion_type.dz + uz .= t_unit .* motion_type.dz + return nothing end function times(motion_type::PeriodicTranslation) diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl index 6b48a6917..2216e7d1d 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(t_start, t_end, pitch, roll, yaw) + rotation = Rotation(pitch, roll, yaw, t_start, t_end) Rotation motion struct. It produces a rotation of the phantom in the three axes: x (pitch), y (roll), and z (yaw). @@ -58,12 +58,13 @@ julia> rt = Rotation(pitch=15.0, roll=0.0, yaw=20.0, t_start=0.1, t_end=0.5) yaw :: T t_start::T = typeof(pitch)(0.0) t_end = typeof(pitch)(0.0) - @assert t_end >= t_start "t_end must be major or equal than t_start" + @assert t_end >= t_start "t_end must be greater or equal than t_start" end is_composable(motion_type::Rotation) = true -function displacement_x( +function displacement_x!( + ux::AbstractArray{T}, motion_type::Rotation{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -74,12 +75,14 @@ function displacement_x( α = t_unit .* (motion_type.yaw) β = t_unit .* (motion_type.roll) γ = t_unit .* (motion_type.pitch) - return cosd.(α) .* cosd.(β) .* x + - (cosd.(α) .* sind.(β) .* sind.(γ) .- sind.(α) .* cosd.(γ)) .* y + - (cosd.(α) .* sind.(β) .* cosd.(γ) .+ sind.(α) .* sind.(γ)) .* z .- x + ux .= cosd.(α) .* cosd.(β) .* x + + (cosd.(α) .* sind.(β) .* sind.(γ) .- sind.(α) .* cosd.(γ)) .* y + + (cosd.(α) .* sind.(β) .* cosd.(γ) .+ sind.(α) .* sind.(γ)) .* z .- x + return nothing end -function displacement_y( +function displacement_y!( + uy::AbstractArray{T}, motion_type::Rotation{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -90,12 +93,14 @@ function displacement_y( α = t_unit .* (motion_type.yaw) β = t_unit .* (motion_type.roll) γ = t_unit .* (motion_type.pitch) - return sind.(α) .* cosd.(β) .* x + - (sind.(α) .* sind.(β) .* sind.(γ) .+ cosd.(α) .* cosd.(γ)) .* y + - (sind.(α) .* sind.(β) .* cosd.(γ) .- cosd.(α) .* sind.(γ)) .* z .- y + uy .= sind.(α) .* cosd.(β) .* x + + (sind.(α) .* sind.(β) .* sind.(γ) .+ cosd.(α) .* cosd.(γ)) .* y + + (sind.(α) .* sind.(β) .* cosd.(γ) .- cosd.(α) .* sind.(γ)) .* z .- y + return nothing end -function displacement_z( +function displacement_z!( + uz::AbstractArray{T}, motion_type::Rotation{T}, x::AbstractArray{T}, y::AbstractArray{T}, @@ -106,9 +111,10 @@ function displacement_z( α = t_unit .* (motion_type.yaw) β = t_unit .* (motion_type.roll) γ = t_unit .* (motion_type.pitch) - return -sind.(β) .* x + + uz .= -sind.(β) .* x + cosd.(β) .* sind.(γ) .* y + cosd.(β) .* cosd.(γ) .* z .- z + return nothing end times(motion_type::Rotation) = begin diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl index f0832e3c1..bbdca7317 100644 --- a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl @@ -1,5 +1,5 @@ @doc raw""" - translation = Translation(t_start, t_end, dx, dy, dz) + translation = Translation(dx, dy, dz, t_start, t_end) Translation motion struct. It produces a linear translation of the phantom. Its fields are the final displacements in the three axes (dx, dy, dz) @@ -26,10 +26,11 @@ julia> tr = Translation(dx=0.01, dy=0.02, dz=0.03, t_start=0.0, t_end=0.5) dz :: T t_start::T = typeof(dx)(0.0) t_end::T = typeof(dx)(0.0) - @assert t_end >= t_start "t_end must be major or equal than t_start" + @assert t_end >= t_start "t_end must be greater or equal than t_start" end -function displacement_x( +function displacement_x!( + ux::AbstractArray{T}, motion_type::Translation{T}, x::AbstractVector{T}, y::AbstractVector{T}, @@ -37,10 +38,12 @@ function displacement_x( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) - return t_unit .* motion_type.dx + ux .= t_unit .* motion_type.dx + return nothing end -function displacement_y( +function displacement_y!( + uy::AbstractArray{T}, motion_type::Translation{T}, x::AbstractVector{T}, y::AbstractVector{T}, @@ -48,10 +51,12 @@ function displacement_y( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) - return t_unit .* motion_type.dy + uy .= t_unit .* motion_type.dy + return nothing end -function displacement_z( +function displacement_z!( + uz::AbstractArray{T}, motion_type::Translation{T}, x::AbstractVector{T}, y::AbstractVector{T}, @@ -59,7 +64,8 @@ function displacement_z( t::AbstractArray{T}, ) where {T<:Real} t_unit = unit_time(t, motion_type.t_start, motion_type.t_end) - return t_unit .* motion_type.dz + uz .= t_unit .* motion_type.dz + return nothing end times(motion_type::Translation) = begin diff --git a/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/sorting.jl b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/sorting.jl new file mode 100644 index 000000000..17a4931a5 --- /dev/null +++ b/KomaMRIBase/src/datatypes/phantom/motion/simplemotion/sorting.jl @@ -0,0 +1,41 @@ +""" + sorted_motion = sort_motion(motion) + +Sorts, with respect to time, the motion types of a `SimpleMotion` instance. +No allocations, since it uses the TupleTools.sort method +""" +function sort_motion(motion::SimpleMotion) + return SimpleMotion(_sort(motion.types, isless, m -> times(m)[1])) +end + + +""" + _sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) -> ::Tuple + +Sorts the tuple `t`. Extracted from TupleTools.jl +""" +@inline function _sort(t::Tuple, lt=isless, by=identity, rev::Bool=false) + t1, t2 = _split(t) + t1s = _sort(t1, lt, by, rev) + t2s = _sort(t2, lt, by, rev) + return _merge(t1s, t2s, lt, by, rev) +end +_sort(t::Tuple{Any}, lt=isless, by=identity, rev::Bool=false) = t +_sort(t::Tuple{}, lt=isless, by=identity, rev::Bool=false) = t + +function _split(t::Tuple) + N = length(t) + M = N >> 1 + return ntuple(i -> t[i], M), ntuple(i -> t[i + M], N - M) +end + +function _merge(t1::Tuple, t2::Tuple, lt, by, rev) + if lt(by(first(t1)), by(first(t2))) != rev + return (first(t1), _merge(tail(t1), t2, lt, by, rev)...) + else + return (first(t2), _merge(t1, tail(t2), lt, by, rev)...) + end +end +_merge(::Tuple{}, t2::Tuple, lt, by, rev) = t2 +_merge(t1::Tuple, ::Tuple{}, lt, by, rev) = t1 +_merge(::Tuple{}, ::Tuple{}, lt, by, rev) = () \ No newline at end of file diff --git a/KomaMRIBase/src/timing/UnitTime.jl b/KomaMRIBase/src/timing/UnitTime.jl new file mode 100644 index 000000000..0dff6bb7e --- /dev/null +++ b/KomaMRIBase/src/timing/UnitTime.jl @@ -0,0 +1,79 @@ +""" + t_unit = unit_time(t, t_start, t_end) + +The `unit_time` function normalizes a given array of time values t +to a unit interval [0, 1] based on a specified start time t_start and end time t_end. +This function is used for non-periodic motions, where each element of t is transformed +to fit within the range [0, 1] based on the provided start and end times. + +![Unit Time](../assets/unit-time.svg) + +# Arguments +- `t`: (`::AbstractArray{T<:Real}`, `[s]`) array of time values to be normalized +- `t_start`: (`::T`, `[s]`) start time for normalization +- `t_end`: (`::T`, `[s]`) end time for normalization + +# Returns +- `t_unit`: (`::AbstractArray{T<:Real}`, `[s]`) array of normalized time values + +# Examples +```julia-repl +julia> t_unit = KomaMRIBase.unit_time([0.0, 1.0, 2.0, 3.0, 4.0, 5.0], 1.0, 4.0) +6-element Vector{Float64}: + 0.0 + 0.0 + 0.3333333333333333 + 0.6666666666666666 + 1.0 + 1.0 +""" +function unit_time(t::AbstractArray{T}, t_start::T, t_end::T) where {T<:Real} + if t_start == t_end + return (t .>= t_start) .* oneunit(T) + else + return min.(max.((t .- t_start) ./ (t_end - t_start), zero(T)), oneunit(T)) + end +end + +""" + t_unit = unit_time_triangular(t, period, asymmetry) + +The `unit_time_triangular` function normalizes a given array +of time values t to a unit interval [0, 1] for periodic motions, +based on a specified period and an asymmetry factor. +This function is useful for creating triangular waveforms +or normalizing time values in periodic processes. + +![Unit Time Triangular](../assets/unit-time-triangular.svg) + +# Arguments +- `t`: (`::AbstractArray{T<:Real}`, `[s]`) array of time values to be normalized +- `period`: (`::T`, `[s]`) the period of the triangular waveform +- `asymmetry`: (`::T`) asymmetry factor, a value in the range (0, 1) indicating the fraction of the period in the rising part of the triangular wave + +# Returns +- `t_unit`: (`::AbstractArray{T<:Real}`, `[s]`) array of normalized time values + +# Examples +```julia-repl +julia> t_unit = KomaMRIBase.unit_time_triangular([0.0, 1.0, 2.0, 3.0, 4.0, 5.0], 4.0, 0.5) +6-element Vector{Float64}: + 0.0 + 0.5 + 1.0 + 0.5 + 0.0 + 0.5 +""" +function unit_time_triangular(t::AbstractArray{T}, period::T, asymmetry::T) where {T<:Real} + t_rise = period * asymmetry + t_fall = period * (oneunit(T) - asymmetry) + t_relative = mod.(t, period) + t_unit = + ifelse.( + t_relative .< t_rise, + t_relative ./ t_rise, + oneunit(T) .- (t_relative .- t_rise) ./ t_fall, + ) + return t_unit +end diff --git a/KomaMRIBase/test/runtests.jl b/KomaMRIBase/test/runtests.jl index f15bfaca2..a9e251611 100644 --- a/KomaMRIBase/test/runtests.jl +++ b/KomaMRIBase/test/runtests.jl @@ -382,82 +382,180 @@ end @test size(obj1) == size(ρ) @test length(obj1) == length(ρ) - # Test obtaining spin psositions - @testset "SimpleMotion" begin + # Test obtaining spin positions + @testset "NoMotion" begin + ph = Phantom(x=[1.0, 2.0], y=[1.0, 2.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) + xt, yt, zt = get_spin_coords(ph.motion, ph.x, ph.y, ph.z, t') + @test xt == ph.x + @test yt == ph.y + @test zt == ph.z + end + @testset "Translation" begin ph = Phantom(x=[1.0], y=[1.0]) t_start=0.0; t_end=1.0 t = collect(range(t_start, t_end, 11)) - period = 2.0 - asymmetry = 0.5 - # Translation dx, dy, dz = [1.0, 0.0, 0.0] vx, vy, vz = [dx, dy, dz] ./ (t_end - t_start) - translation = SimpleMotion([Translation(dx, dy, dz, t_start, t_end)]) + translation = SimpleMotion(Translation(dx, dy, dz, t_start, t_end)) xt, yt, zt = get_spin_coords(translation, ph.x, ph.y, ph.z, t') @test xt == ph.x .+ vx.*t' @test yt == ph.y .+ vy.*t' @test zt == ph.z .+ vz.*t' - # PeriodicTranslation - periodictranslation = SimpleMotion([PeriodicTranslation(dx, dy, dz, period, asymmetry)]) + # ----- t_start = t_end -------- + t_start = t_end = 0.0 + t = [-0.5, -0.25, 0.0, 0.25, 0.5] + translation = SimpleMotion(Translation(dx, dy, dz, t_start, t_end)) + xt, yt, zt = get_spin_coords(translation, ph.x, ph.y, ph.z, t') + @test xt == ph.x .+ dx*[0, 0, 1, 1, 1]' + @test yt == ph.y .+ dy*[0, 0, 1, 1, 1]' + @test zt == ph.z .+ dz*[0, 0, 1, 1, 1]' + end + @testset "PeriodicTranslation" begin + ph = Phantom(x=[1.0], y=[1.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) + period = 2.0 + asymmetry = 0.5 + dx, dy, dz = [1.0, 0.0, 0.0] + vx, vy, vz = [dx, dy, dz] ./ (t_end - t_start) + periodictranslation = SimpleMotion(PeriodicTranslation(dx, dy, dz, period, asymmetry)) xt, yt, zt = get_spin_coords(periodictranslation, ph.x, ph.y, ph.z, t') @test xt == ph.x .+ vx.*t' @test yt == ph.y .+ vy.*t' @test zt == ph.z .+ vz.*t' - # Rotation (2D) - pitch = 0.0 + end + @testset "Rotation" begin + ph = Phantom(x=[1.0], y=[1.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) + pitch = 45.0 roll = 0.0 yaw = 45.0 - rotation = SimpleMotion([Rotation(pitch, roll, yaw, t_start, t_end)]) + rotation = SimpleMotion(Rotation(pitch, roll, yaw, t_start, t_end)) xt, yt, zt = get_spin_coords(rotation, ph.x, ph.y, ph.z, t') - @test xt[:,end] == ph.x .* cosd(yaw) - ph.y .* sind(yaw) - @test yt[:,end] == ph.x .* sind(yaw) + ph.y .* cosd(yaw) - @test zt[:,end] == ph.z - # PeriodicRotation (2D) - periodicrotation = SimpleMotion([PeriodicRotation(pitch, roll, yaw, period, asymmetry)]) + r = vcat(ph.x, ph.y, ph.z) + R = rotz(π*yaw/180) * roty(π*roll/180) * rotx(π*pitch/180) + rot_x, rot_y, rot_z = R*r + @test xt[end ,end] ≈ rot_x + @test yt[end ,end] ≈ rot_y + @test zt[end ,end] ≈ rot_z + # ----- t_start = t_end -------- + t_start = t_end = 0.0 + t = [-0.5, -0.25, 0.0, 0.25, 0.5] + rotation = SimpleMotion(Rotation(pitch, roll, yaw, t_start, t_end)) + xt, yt, zt = get_spin_coords(rotation, ph.x, ph.y, ph.z, t') + @test xt ≈ [ph.x ph.x rot_x rot_x rot_x] + @test yt ≈ [ph.y ph.y rot_y rot_y rot_y] + @test zt ≈ [ph.z ph.z rot_z rot_z rot_z] + end + @testset "PeriodicRotation" begin + ph = Phantom(x=[1.0], y=[1.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) + period = 2.0 + asymmetry = 0.5 + pitch = 45.0 + roll = 0.0 + yaw = 45.0 + periodicrotation = SimpleMotion(PeriodicRotation(pitch, roll, yaw, period, asymmetry)) xt, yt, zt = get_spin_coords(periodicrotation, ph.x, ph.y, ph.z, t') - @test xt[:,end] == ph.x .* cosd(yaw) - ph.y .* sind(yaw) - @test yt[:,end] == ph.x .* sind(yaw) + ph.y .* cosd(yaw) - @test zt[:,end] == ph.z - # HeartBeat + r = vcat(ph.x, ph.y, ph.z) + R = rotz(π*yaw/180) * roty(π*roll/180) * rotx(π*pitch/180) + rot_x, rot_y, rot_z = R*r + @test xt[end ,end] ≈ rot_x + @test yt[end ,end] ≈ rot_y + @test zt[end ,end] ≈ rot_z + end + @testset "HeartBeat" begin + ph = Phantom(x=[1.0], y=[1.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) circumferential_strain = -0.1 radial_strain = 0.0 longitudinal_strain = -0.1 - heartbeat = SimpleMotion([HeartBeat(circumferential_strain, radial_strain, longitudinal_strain, t_start, t_end)]) + heartbeat = SimpleMotion(HeartBeat(circumferential_strain, radial_strain, longitudinal_strain, t_start, t_end)) xt, yt, zt = get_spin_coords(heartbeat, ph.x, ph.y, ph.z, t') r = sqrt.(ph.x .^ 2 + ph.y .^ 2) θ = atan.(ph.y, ph.x) @test xt[:,end] == ph.x .* (1 .+ circumferential_strain * maximum(r) .* cos.(θ)) @test yt[:,end] == ph.y .* (1 .+ circumferential_strain * maximum(r) .* sin.(θ)) @test zt[:,end] == ph.z .* (1 .+ longitudinal_strain) - # PeriodicHeartBeat - periodicheartbeat = SimpleMotion([PeriodicHeartBeat(circumferential_strain, radial_strain, longitudinal_strain, period, asymmetry)]) + # ----- t_start = t_end -------- + t_start = t_end = 0.0 + t = [-0.5, -0.25, 0.0, 0.25, 0.5] + heartbeat = SimpleMotion(HeartBeat(circumferential_strain, radial_strain, longitudinal_strain, t_start, t_end)) xt, yt, zt = get_spin_coords(heartbeat, ph.x, ph.y, ph.z, t') + r = sqrt.(ph.x .^ 2 + ph.y .^ 2) + θ = atan.(ph.y, ph.x) + dx = (1 .+ circumferential_strain * maximum(r) .* cos.(θ)) + dy = (1 .+ circumferential_strain * maximum(r) .* sin.(θ)) + dz = (1 .+ longitudinal_strain) + @test xt == [ph.x ph.x ph.x .* dx ph.x .* dx ph.x .* dx] + @test yt == [ph.y ph.y ph.y .* dy ph.y .* dy ph.y .* dy] + @test zt == [ph.z ph.z ph.z .* dz ph.z .* dz ph.z .* dz] + end + @testset "PeriodicHeartBeat" begin + ph = Phantom(x=[1.0], y=[1.0]) + t_start=0.0; t_end=1.0 + t = collect(range(t_start, t_end, 11)) + period = 2.0 + asymmetry = 0.5 + circumferential_strain = -0.1 + radial_strain = 0.0 + longitudinal_strain = -0.1 + periodicheartbeat = SimpleMotion(PeriodicHeartBeat(circumferential_strain, radial_strain, longitudinal_strain, period, asymmetry)) + 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) @test xt[:,end] == ph.x .* (1 .+ circumferential_strain * maximum(r) .* cos.(θ)) @test yt[:,end] == ph.y .* (1 .+ circumferential_strain * maximum(r) .* sin.(θ)) @test zt[:,end] == ph.z .* (1 .+ longitudinal_strain) end @testset "ArbitraryMotion" begin + # 1 spin ph = Phantom(x=[1.0], y=[1.0]) Ns = length(ph) - period_durations = [1.0] - num_pieces = 10 - dx = dy = dz = rand(Ns, num_pieces - 1) - arbitrarymotion = @suppress ArbitraryMotion(period_durations, dx, dy, dz) + t_start = 0.0 + t_end = 1.0 + Nt = 10 + dx = rand(Ns, Nt) + dy = rand(Ns, Nt) + dz = rand(Ns, Nt) + arbitrarymotion = @suppress ArbitraryMotion(t_start, t_end, dx, dy, dz) + t = times(arbitrarymotion) + xt, yt, zt = get_spin_coords(arbitrarymotion, ph.x, ph.y, ph.z, t') + @test xt == ph.x .+ dx + @test yt == ph.y .+ dy + @test zt == ph.z .+ dz + # More than 1 spin + ph = Phantom(x=[1.0, 2.0], y=[1.0, 2.0]) + Ns = length(ph) + t_start = 0.0 + t_end = 1.0 + Nt = 10 + dx = rand(Ns, Nt) + dy = rand(Ns, Nt) + dz = rand(Ns, Nt) + arbitrarymotion = @suppress ArbitraryMotion(t_start, t_end, dx, dy, dz) t = times(arbitrarymotion) xt, yt, zt = get_spin_coords(arbitrarymotion, ph.x, ph.y, ph.z, t') - @test xt[:,2:end-1] == ph.x .+ dx - @test yt[:,2:end-1] == ph.y .+ dy - @test zt[:,2:end-1] == ph.z .+ dz + @test xt == ph.x .+ dx + @test yt == ph.y .+ dy + @test zt == ph.z .+ dz end - simplemotion = SimpleMotion([ + simplemotion = SimpleMotion( PeriodicTranslation(dx=0.05, dy=0.05, dz=0.0, period=0.5, asymmetry=0.5), Rotation(pitch=0.0, roll=0.0, yaw=π / 2, t_start=0.05, t_end=0.5), - ]) + ) Ns = length(obj1) - K = 10 - arbitrarymotion = @suppress ArbitraryMotion([1.0], 0.01 .* rand(Ns, K - 1), 0.01 .* rand(Ns, K - 1), 0.01 .* rand(Ns, K - 1)) + Nt = 3 + t_start = 0.0 + t_end = 1.0 + arbitrarymotion = @suppress ArbitraryMotion(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/ext/KomaAMDGPUExt.jl b/KomaMRICore/ext/KomaAMDGPUExt.jl index 112dc2285..d0b5f1350 100644 --- a/KomaMRICore/ext/KomaAMDGPUExt.jl +++ b/KomaMRICore/ext/KomaAMDGPUExt.jl @@ -10,12 +10,6 @@ KomaMRICore.set_device!(::ROCBackend, dev_idx::Integer) = AMDGPU.device_id!(dev_ KomaMRICore.set_device!(::ROCBackend, dev::AMDGPU.HIPDevice) = AMDGPU.device!(dev) KomaMRICore.device_name(::ROCBackend) = AMDGPU.HIP.name(AMDGPU.device()) -function Adapt.adapt_storage( - ::ROCBackend, x::Vector{KomaMRICore.LinearInterpolator{T,V}} -) where {T<:Real,V<:AbstractVector{T}} - return AMDGPU.rocconvert.(x) -end - function KomaMRICore._print_devices(::ROCBackend) devices = [ Symbol("($(i-1)$(i == 1 ? "*" : " "))") => AMDGPU.HIP.name(d) for diff --git a/KomaMRICore/ext/KomaCUDAExt.jl b/KomaMRICore/ext/KomaCUDAExt.jl index 9f8ba96e7..cc97a0a96 100644 --- a/KomaMRICore/ext/KomaCUDAExt.jl +++ b/KomaMRICore/ext/KomaCUDAExt.jl @@ -9,12 +9,6 @@ KomaMRICore.isfunctional(::CUDABackend) = CUDA.functional() KomaMRICore.set_device!(::CUDABackend, val) = CUDA.device!(val) KomaMRICore.device_name(::CUDABackend) = CUDA.name(CUDA.device()) -function Adapt.adapt_storage( - ::CUDABackend, x::Vector{KomaMRICore.LinearInterpolator{T,V}} -) where {T<:Real,V<:AbstractVector{T}} - return Adapt.adapt.(CuArray, x) -end - function KomaMRICore._print_devices(::CUDABackend) devices = [ Symbol("($(i-1)$(i == 1 ? "*" : " "))") => CUDA.name(d) for diff --git a/KomaMRICore/ext/KomaMetalExt.jl b/KomaMRICore/ext/KomaMetalExt.jl index bc452f31e..d50af7901 100644 --- a/KomaMRICore/ext/KomaMetalExt.jl +++ b/KomaMRICore/ext/KomaMetalExt.jl @@ -12,12 +12,6 @@ KomaMRICore.set_device!(::MetalBackend, device_index::Integer) = device_index == KomaMRICore.set_device!(::MetalBackend, dev::Metal.MTLDevice) = Metal.device!(dev) KomaMRICore.device_name(::MetalBackend) = String(Metal.current_device().name) -function Adapt.adapt_storage( - ::MetalBackend, x::Vector{KomaMRICore.LinearInterpolator{T,V}} -) where {T<:Real,V<:AbstractVector{T}} - return Metal.mtl.(x) -end - function KomaMRICore._print_devices(::MetalBackend) @info "Metal device type: $(KomaMRICore.device_name(MetalBackend()))" end diff --git a/KomaMRICore/ext/KomaoneAPIExt.jl b/KomaMRICore/ext/KomaoneAPIExt.jl index 03cc270de..0bec755fc 100644 --- a/KomaMRICore/ext/KomaoneAPIExt.jl +++ b/KomaMRICore/ext/KomaoneAPIExt.jl @@ -9,12 +9,6 @@ KomaMRICore.isfunctional(::oneAPIBackend) = oneAPI.functional() KomaMRICore.set_device!(::oneAPIBackend, val) = oneAPI.device!(val) KomaMRICore.device_name(::oneAPIBackend) = oneAPI.properties(oneAPI.device()).name -function Adapt.adapt_storage( - ::oneAPIBackend, x::Vector{KomaMRICore.LinearInterpolator{T,V}} -) where {T<:Real,V<:AbstractVector{T}} - return Adapt.adapt.(oneArray, a) -end - function KomaMRICore._print_devices(::oneAPIBackend) devices = [ Symbol("($(i-1)$(i == 1 ? "*" : " "))") => oneAPI.properties(d).name for diff --git a/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl b/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl index f60a8b278..a6dabbd0a 100644 --- a/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl +++ b/KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl @@ -19,7 +19,6 @@ function initialize_spins_state( Mxy = zeros(T, Nspins) Mz = obj.ρ Xt = Mag{T}(Mxy, Mz) - sort_motions!(obj.motion) return Xt, obj end diff --git a/KomaMRICore/src/simulation/Functors.jl b/KomaMRICore/src/simulation/Functors.jl index 4aac6628f..0d698b78f 100644 --- a/KomaMRICore/src/simulation/Functors.jl +++ b/KomaMRICore/src/simulation/Functors.jl @@ -71,36 +71,13 @@ x = x |> cpu """ cpu(x) = fmap(x -> adapt(KA.CPU(), x), x, exclude=_isleaf) -#MotionModel structs -adapt_storage(::KA.GPU, x::NoMotion) = x -adapt_storage(::KA.GPU, x::SimpleMotion) = x -function adapt_storage(backend::KA.GPU, xs::ArbitraryMotion) - fields = [] - for field in fieldnames(ArbitraryMotion) - if field in (:ux, :uy, :uz) - push!(fields, adapt(backend, getfield(xs, field))) - else - push!(fields, getfield(xs, field)) - end - end - return KomaMRICore.ArbitraryMotion(fields...) -end - #Precision paramtype(T::Type{<:Real}, m) = fmap(x -> adapt(T, x), m) adapt_storage(T::Type{<:Real}, xs::Real) = convert(T, xs) adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Real}) = convert.(T, xs) adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Complex}) = convert.(Complex{T}, xs) adapt_storage(T::Type{<:Real}, xs::AbstractArray{<:Bool}) = xs -adapt_storage(T::Type{<:Real}, xs::SimpleMotion) = SimpleMotion(paramtype(T, xs.types)) adapt_storage(T::Type{<:Real}, xs::NoMotion) = NoMotion{T}() -function adapt_storage(T::Type{<:Real}, xs::ArbitraryMotion) - fields = [] - for field in fieldnames(ArbitraryMotion) - push!(fields, paramtype(T, getfield(xs, field))) - end - return ArbitraryMotion(fields...) -end """ f32(m) @@ -123,16 +100,21 @@ See also [`f32`](@ref). f64(m) = paramtype(Float64, m) #The functor macro makes it easier to call a function in all the parameters +# Phantom @functor Phantom - +# SimpleMotion +@functor SimpleMotion @functor Translation @functor Rotation @functor HeartBeat @functor PeriodicTranslation @functor PeriodicRotation @functor PeriodicHeartBeat - +# ArbitraryMotion +@functor ArbitraryMotion +# Spinor @functor Spinor +# DiscreteSequence @functor DiscreteSequence export gpu, cpu, f32, f64 \ No newline at end of file diff --git a/KomaMRICore/test/runtests.jl b/KomaMRICore/test/runtests.jl index 315d37f1a..77246c48e 100644 --- a/KomaMRICore/test/runtests.jl +++ b/KomaMRICore/test/runtests.jl @@ -295,8 +295,7 @@ end sig_jemris = signal_brain_motion_jemris() seq = seq_epi_100x100_TE100_FOV230() sys = Scanner() - obj = phantom_brain() - obj.motion = SimpleMotion([Translation(t_end=10.0, dx=0.0, dy=1.0, dz=0.0)]) + obj = phantom_brain_simple_motion() sim_params = Dict{String, Any}( "gpu"=>USE_GPU, "sim_method"=>KomaMRICore.Bloch(), @@ -316,16 +315,7 @@ end sig_jemris = signal_brain_motion_jemris() seq = seq_epi_100x100_TE100_FOV230() sys = Scanner() - obj = phantom_brain() - Ns = length(obj) - period_durations=[20.0] - dx = dz = zeros(Ns, 1) - dy = 1.0 .* ones(Ns, 1) - obj.motion = @suppress ArbitraryMotion( - period_durations, - dx, - dy, - dz) + obj = phantom_brain_arbitrary_motion() sim_params = Dict{String, Any}( "gpu"=>USE_GPU, "sim_method"=>KomaMRICore.Bloch(), diff --git a/KomaMRICore/test/test_files/utils.jl b/KomaMRICore/test/test_files/utils.jl index d10066810..0a5f878a6 100644 --- a/KomaMRICore/test/test_files/utils.jl +++ b/KomaMRICore/test/test_files/utils.jl @@ -16,6 +16,29 @@ function signal_brain_motion_jemris() return sig end +function phantom_brain_simple_motion() + obj = phantom_brain() + obj.motion = SimpleMotion(Translation(t_end=10.0, dx=0.0, dy=1.0, dz=0.0)) + return obj +end + +function phantom_brain_arbitrary_motion() + obj = phantom_brain() + Ns = length(obj) + t_start = 0.0 + t_end = 10.0 + dx = zeros(Ns, 2) + dz = zeros(Ns, 2) + dy = [zeros(Ns,1) ones(Ns,1)] + obj.motion = ArbitraryMotion( + t_start, + t_end, + dx, + dy, + dz) + return obj +end + function phantom_sphere() path = @__DIR__ fid = h5open(joinpath(path, "phantom_sphere.h5"), "r") diff --git a/KomaMRIFiles/src/Phantom/Phantom.jl b/KomaMRIFiles/src/Phantom/Phantom.jl index cd2bf15d8..8409d792f 100644 --- a/KomaMRIFiles/src/Phantom/Phantom.jl +++ b/KomaMRIFiles/src/Phantom/Phantom.jl @@ -5,13 +5,13 @@ Reads a (.phantom) file and creates a Phantom structure from it """ function read_phantom(filename::String) fid = HDF5.h5open(filename, "r") - fields = [] + phantom_fields = [] version = read_attribute(fid, "Version") dims = read_attribute(fid, "Dims") Ns = read_attribute(fid, "Ns") # Name name = read_attribute(fid, "Name") - push!(fields, (:name, name)) + push!(phantom_fields, (:name, name)) # Position and contrast for key in ["position", "contrast"] group = fid[key] @@ -19,16 +19,16 @@ function read_phantom(filename::String) param = group[label] values = read_param(param) if values != "Default" - push!(fields, (Symbol(label), values)) + push!(phantom_fields, (Symbol(label), values)) end end end # Motion motion_group = fid["motion"] model = read_attribute(motion_group, "model") - import_motion!(fields, Ns, Symbol(model), motion_group) + import_motion!(phantom_fields, Ns, Symbol(model), motion_group) - obj = Phantom(; fields...) + obj = Phantom(; phantom_fields...) close(fid) return obj end @@ -54,16 +54,16 @@ function read_param(param::HDF5.Group) return values end -function import_motion!(fields::Array, Ns::Int, model::Symbol, motion_group::HDF5.Group) - return import_motion!(fields, Ns, Val(model), motion_group) +function import_motion!(phantom_fields::Array, Ns::Int, model::Symbol, motion_group::HDF5.Group) + return import_motion!(phantom_fields, Ns, Val(model), motion_group) end function import_motion!( - fields::Array, Ns::Int, model::Val{:NoMotion}, motion_group::HDF5.Group + phantom_fields::Array, Ns::Int, model::Val{:NoMotion}, motion_group::HDF5.Group ) return nothing end function import_motion!( - fields::Array, Ns::Int, model::Val{:SimpleMotion}, motion_group::HDF5.Group + phantom_fields::Array, Ns::Int, model::Val{:SimpleMotion}, motion_group::HDF5.Group ) types_group = motion_group["types"] types = SimpleMotionType[] @@ -81,29 +81,17 @@ function import_motion!( end end end - return push!(fields, (:motion, SimpleMotion(vcat(types...)))) + return push!(phantom_fields, (:motion, SimpleMotion((types...)))) end function import_motion!( - fields::Array, Ns::Int, model::Val{:ArbitraryMotion}, motion_group::HDF5.Group + phantom_fields::Array, Ns::Int, model::Val{:ArbitraryMotion}, motion_group::HDF5.Group ) - dur = read(motion_group["period_durations"]) - K = read_attribute(motion_group, "K") - dx = zeros(Ns, K - 1) - dy = zeros(Ns, K - 1) - dz = zeros(Ns, K - 1) - for key in HDF5.keys(motion_group) - if key != "period_durations" - values = read_param(motion_group[key]) - if key == "dx" - dx = values - elseif key == "dy" - dy = values - elseif key == "dz" - dz = values - end - end - end - return push!(fields, (:motion, ArbitraryMotion(dur, dx, dy, dz))) + t_start = read(motion_group["t_start"]) + t_end = read(motion_group["t_end"]) + dx = read(motion_group["dx"]) + dy = read(motion_group["dy"]) + dz = read(motion_group["dz"]) + return push!(phantom_fields, (:motion, ArbitraryMotion(t_start, t_end, dx, dy, dz))) end """ @@ -158,20 +146,17 @@ function export_motion!(motion_group::HDF5.Group, motion::SimpleMotion) for (counter, sm_type) in enumerate(motion.types) simple_motion_type = typeof(sm_type).name.name type_group = create_group(types_group, "$(counter)_$simple_motion_type") - fields = fieldnames(typeof(sm_type)) - for field in fields + phantom_fields = fieldnames(typeof(sm_type)) + for field in phantom_fields HDF5.attributes(type_group)[string(field)] = getfield(sm_type, field) end end end function export_motion!(motion_group::HDF5.Group, motion::ArbitraryMotion) HDF5.attributes(motion_group)["model"] = "ArbitraryMotion" - HDF5.attributes(motion_group)["K"] = size(motion.dx)[2] + 1 - motion_group["period_durations"] = motion.period_durations - - for key in ["dx", "dy", "dz"] - d_group = create_group(motion_group, key) - HDF5.attributes(d_group)["type"] = "Explicit" #TODO: consider "Indexed" type - d_group["values"] = getfield(motion, Symbol(key)) - end + motion_group["t_start"] = motion.t_start + motion_group["t_end"] = motion.t_end + motion_group["dx"] = motion.dx + motion_group["dy"] = motion.dy + motion_group["dz"] = motion.dz end \ No newline at end of file diff --git a/KomaMRIFiles/test/runtests.jl b/KomaMRIFiles/test/runtests.jl index deccbf82d..02d9a432f 100644 --- a/KomaMRIFiles/test/runtests.jl +++ b/KomaMRIFiles/test/runtests.jl @@ -58,10 +58,13 @@ using TestItems, TestItemRunner write_phantom(obj1, filename) obj2 = read_phantom(filename) @test obj1 == obj2 + end + @testset "SimpleMotion" begin # SimpleMotion + path = @__DIR__ filename = path * "/test_files/brain_simplemotion.phantom" obj1 = brain_phantom2D() - obj1.motion = SimpleMotion([ + obj1.motion = SimpleMotion( PeriodicRotation( period=1.0, yaw=45.0, @@ -73,17 +76,24 @@ using TestItems, TestItemRunner dx=0.0, dy=0.02, dz=0.0 - )]) + ) + ) write_phantom(obj1, filename) obj2 = read_phantom(filename) @test obj1 == obj2 + end + @testset "ArbitraryMotion" begin # ArbitraryMotion + path = @__DIR__ filename = path * "/test_files/brain_arbitrarymotion.phantom" obj1 = brain_phantom2D() Ns = length(obj1) K = 10 + t_start = 0.0 + t_end = 1.0 obj1.motion = ArbitraryMotion( - [1.0], + t_start, + t_end, 0.01.*rand(Ns, K-1), 0.01.*rand(Ns, K-1), 0.01.*rand(Ns, K-1)) diff --git a/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom b/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom index e33788427..2010fcf20 100644 Binary files a/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom and b/KomaMRIFiles/test/test_files/brain_arbitrarymotion.phantom differ diff --git a/KomaMRIFiles/test/test_files/brain_nomotion.phantom b/KomaMRIFiles/test/test_files/brain_nomotion.phantom index a487a64b8..cf0d1bae2 100644 Binary files a/KomaMRIFiles/test/test_files/brain_nomotion.phantom and b/KomaMRIFiles/test/test_files/brain_nomotion.phantom differ diff --git a/KomaMRIFiles/test/test_files/brain_simplemotion.phantom b/KomaMRIFiles/test/test_files/brain_simplemotion.phantom index c72c99584..cfc0a7d34 100644 Binary files a/KomaMRIFiles/test/test_files/brain_simplemotion.phantom and b/KomaMRIFiles/test/test_files/brain_simplemotion.phantom differ diff --git a/KomaMRIPlots/src/ui/DisplayFunctions.jl b/KomaMRIPlots/src/ui/DisplayFunctions.jl index 9bc6c550a..4d181c409 100644 --- a/KomaMRIPlots/src/ui/DisplayFunctions.jl +++ b/KomaMRIPlots/src/ui/DisplayFunctions.jl @@ -792,8 +792,8 @@ function plot_image( image; height=600, width=nothing, - zmin=minimum(abs.(image[:])), - zmax=maximum(abs.(image[:])), + zmin=minimum(image[:]), + zmax=maximum(image[:]), darkmode=false, title="", ) @@ -1032,36 +1032,29 @@ function plot_phantom_map( frame_duration_ms=250, kwargs..., ) - function process_times(motion::SimpleMotion) + + function interpolate_times(motion::MotionModel) t = times(motion) # Interpolate time points (as many as indicated by intermediate_time_samples) - itp = interpolate( - ( - 1:(intermediate_time_samples + 1):(length(t) + intermediate_time_samples * (length(t) - 1)), - ), - t, - Gridded(Linear()), - ) - return itp.(1:(length(t) + intermediate_time_samples * (length(t) - 1))) + itp = interpolate((1:(intermediate_time_samples + 1):(length(t) + intermediate_time_samples * (length(t) - 1)), ), t, Gridded(Linear())) + t = itp.(1:(length(t) + intermediate_time_samples * (length(t) - 1))) + return t + end + + function process_times(motion::SimpleMotion) + motion = KomaMRIBase.sort_motion(motion) + return interpolate_times(motion) end function process_times(motion::ArbitraryMotion) - t = times(motion) - # Interpolate time points (as many as indicated by intermediate_time_samples) - itp = interpolate( - ( - 1:(intermediate_time_samples + 1):(length(t) + intermediate_time_samples * (length(t) - 1)), - ), - t, - Gridded(Linear()), - ) - t = itp.(1:(length(t) + intermediate_time_samples * (length(t) - 1))) + t = interpolate_times(motion) # Decimate time points so their number is smaller than max_time_samples ss = length(t) > max_time_samples ? length(t) ÷ max_time_samples : 1 return t[1:ss:end] end - process_times(motion::MotionModel) = times(motion) + function process_times(motion::MotionModel) + return times(motion) + end - # IDEA: subsample phantoms which are too large function decimate_uniform_phantom(ph::Phantom, num_points::Int) dimx, dimy, dimz = KomaMRIBase.get_dims(ph) ss = Int(ceil((length(ph) / num_points)^(1 / sum(KomaMRIBase.get_dims(ph))))) @@ -1077,7 +1070,7 @@ function plot_phantom_map( if length(ph) > max_spins ph = decimate_uniform_phantom(ph, max_spins) - @warn "For performance reasons, the number of displayed spins was capped to `max_spins`=$(max_spins)." maxlog=1 + @warn "For performance reasons, the number of displayed spins was capped to `max_spins`=$(max_spins)." end path = @__DIR__ @@ -1129,7 +1122,6 @@ function plot_phantom_map( cmin_key = get(kwargs, :cmin, factor * cmin_key) cmax_key = get(kwargs, :cmax, factor * cmax_key) - sort_motions!(ph.motion) t = process_times(ph.motion) x, y, z = get_spin_coords(ph.motion, ph.x, ph.y, ph.z, t') diff --git a/docs/src/assets/unit-time-triangular.svg b/docs/src/assets/unit-time-triangular.svg new file mode 100644 index 000000000..5a91bd907 --- /dev/null +++ b/docs/src/assets/unit-time-triangular.svg @@ -0,0 +1,355 @@ + + + + + + + + + + + + + + + + Periodic + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + tunit + t + + + + + period + + + + + + + + period(1-asymmetry) + + + + period asymmetry + · + + 1 + + + + + + diff --git a/docs/src/assets/unit-time.svg b/docs/src/assets/unit-time.svg new file mode 100644 index 000000000..c7cc6d3f9 --- /dev/null +++ b/docs/src/assets/unit-time.svg @@ -0,0 +1,405 @@ + + + + + + + + + + + + + + + + Non-periodic + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + tunit + t + + + + tstart + tend + + 1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/3.tutorials/lit-05-SimpleMotion.jl b/examples/3.tutorials/lit-05-SimpleMotion.jl index c26a8aa3f..93b3da121 100644 --- a/examples/3.tutorials/lit-05-SimpleMotion.jl +++ b/examples/3.tutorials/lit-05-SimpleMotion.jl @@ -19,9 +19,9 @@ 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([ +obj.motion = SimpleMotion( Translation(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 #md savefig(p1, "../assets/5-phantom1.html") # hide diff --git a/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl b/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl index e2f63f0d9..7e3c9933c 100644 --- a/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl +++ b/examples/5.koma_paper/comparison_accuracy/d.brain_motion/run_koma.jl @@ -1,7 +1,7 @@ using KomaMRI, Suppressor, MAT obj = @suppress read_phantom_jemris("../../../2.phantoms/brain.h5") -obj.uy = (x,y,z,t)-> 0.1f0 * t # Hacer que el fichero brain_motion.phantom tenga este movimiento (con SimpleMotion y ArbitraryMotion) +obj.uy = (x,y,z,t)-> 0.1f0 * t seq = @suppress read_seq("../sequences/EPI/epi_100x100_TE100_FOV230.seq") sys = Scanner() #Time