Skip to content

Commit

Permalink
Suggested changes
Browse files Browse the repository at this point in the history
  • Loading branch information
pvillacorta committed Aug 8, 2024
1 parent 596af23 commit fe9c8c2
Show file tree
Hide file tree
Showing 20 changed files with 180 additions and 127 deletions.
8 changes: 5 additions & 3 deletions KomaMRIBase/src/KomaMRIBase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions KomaMRIBase/src/datatypes/Phantom.jl
Original file line number Diff line number Diff line change
@@ -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")

Expand Down Expand Up @@ -54,7 +54,7 @@ julia> obj.ρ
::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"""
Expand Down Expand Up @@ -192,15 +192,15 @@ function heart_phantom(
Dλ1=Dλ1[ρ .!= 0],
Dλ2=Dλ2[ρ .!= 0],
=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,
Expand Down
47 changes: 24 additions & 23 deletions KomaMRIBase/src/datatypes/phantom/Motion.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand All @@ -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
Expand All @@ -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},
Expand All @@ -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)
Expand All @@ -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
6 changes: 3 additions & 3 deletions KomaMRIBase/src/datatypes/phantom/NoMotion.jl
Original file line number Diff line number Diff line change
@@ -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
Expand Down
12 changes: 6 additions & 6 deletions KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))])
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand All @@ -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

Expand Down
2 changes: 1 addition & 1 deletion KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
@@ -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
Original file line number Diff line number Diff line change
@@ -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}
Expand Down
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
20 changes: 12 additions & 8 deletions KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl
Original file line number Diff line number Diff line change
@@ -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).
Expand Down Expand Up @@ -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
Expand All @@ -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!(
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
Loading

0 comments on commit fe9c8c2

Please sign in to comment.