Skip to content

Commit

Permalink
Solve motion-related bug
Browse files Browse the repository at this point in the history
  • Loading branch information
pvillacorta committed Oct 17, 2024
1 parent 6202b2a commit 86a504c
Show file tree
Hide file tree
Showing 12 changed files with 83 additions and 94 deletions.
3 changes: 2 additions & 1 deletion KomaMRIBase/src/KomaMRIBase.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ include("timing/KeyValuesCalculation.jl")
include("datatypes/Sequence.jl")
include("datatypes/sequence/Delay.jl")
# Motion
include("motion/AbstractMotion.jl")
include("motion/MotionList.jl")
include("motion/NoMotion.jl")
# Phantom
include("datatypes/Phantom.jl")
# Simulator
Expand Down
4 changes: 2 additions & 2 deletions KomaMRIBase/src/datatypes/Phantom.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ a property value representing a spin. This struct serves as an input for the sim
- `Dλ1`: (`::AbstractVector{T<:Real}`) spin Dλ1 (diffusion) parameter vector
- `Dλ2`: (`::AbstractVector{T<:Real}`) spin Dλ2 (diffusion) parameter vector
- `Dθ`: (`::AbstractVector{T<:Real}`) spin Dθ (diffusion) parameter vector
- `motion`: (`::AbstractMotion{T<:Real}`) motion
- `motion`: (`::Union{NoMotion, MotionList{T<:Real}}`) motion
# Returns
- `obj`: (`::Phantom`) Phantom struct
Expand Down Expand Up @@ -47,7 +47,7 @@ julia> obj.ρ
::AbstractVector{T} = zeros(eltype(x), size(x))
#Diff::Vector{DiffusionModel} #Diffusion map
#Motion
motion::AbstractMotion{T} = NoMotion{eltype(x)}()
motion::Union{NoMotion, MotionList{T}} = NoMotion()
end

const NON_STRING_PHANTOM_FIELDS = Iterators.filter(x -> fieldtype(Phantom, x) != String, fieldnames(Phantom))
Expand Down
11 changes: 0 additions & 11 deletions KomaMRIBase/src/motion/AbstractMotion.jl

This file was deleted.

Original file line number Diff line number Diff line change
@@ -1,3 +1,8 @@
include("motionlist/Action.jl")
include("motionlist/SpinSpan.jl")
include("motionlist/TimeSpan.jl")
include("motionlist/Motion.jl")

"""
motionlist = MotionList(motions...)
Expand Down Expand Up @@ -27,7 +32,7 @@ julia> motionlist = MotionList(
)
```
"""
struct MotionList{T<:Real} <: AbstractMotion{T}
struct MotionList{T<:Real}
motions::Vector{<:Motion{T}}
end

Expand All @@ -40,14 +45,14 @@ function Base.getindex(mv::MotionList{T}, p) where {T<:Real}
for m in mv.motions
m[p] !== nothing ? push!(motion_array_aux, m[p]) : nothing
end
return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion{T}()
return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion()
end
function Base.view(mv::MotionList{T}, p) where {T<:Real}
motion_array_aux = Motion{T}[]
for m in mv.motions
@view(m[p]) !== nothing ? push!(motion_array_aux, @view(m[p])) : nothing
end
return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion{T}()
return length(motion_array_aux) > 0 ? MotionList(motion_array_aux) : NoMotion()
end

""" Addition of MotionLists """
Expand Down Expand Up @@ -91,7 +96,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_{\t{spins}}`` rows and `length(t)` columns.
# Arguments
- `motionset`: (`::AbstractMotion{T<:Real}`) phantom motion
- `motion`: (`::Union{NoMotion, MotionList{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 Down Expand Up @@ -133,27 +138,27 @@ function get_spin_coords(
end

"""
times = times(motion)
times = times(motion, T)
"""
function times(ml::MotionList{T}) where {T<:Real}
nodes = reduce(vcat, [times(m) for m in ml.motions])
function times(ml::MotionList{T}, ::Type{T}) where {T<:Real}
nodes = reduce(vcat, [times(m, T) for m in ml.motions])
return unique(sort(nodes))
end

"""
sort_motions!(motionset)
sort_motions!(motion)
Sorts motions in a list according to their starting time. It modifies the original list.
If `motionset::NoMotion`, this function does nothing.
If `motionset::MotionList`, this function sorts its motions.
# Arguments
- `motionset`: (`::AbstractMotion{T<:Real}`) phantom motion
- `motion`: (`::Union{NoMotion, MotionList{T<:Real}}`) phantom motion
# Returns
- `nothing`
"""
function sort_motions!(m::MotionList)
sort!(m.motions; by=m -> times(m)[1])
function sort_motions!(m::MotionList{T}) where {T<:Real}
sort!(m.motions; by=m -> times(m, T)[1])
return nothing
end
53 changes: 53 additions & 0 deletions KomaMRIBase/src/motion/NoMotion.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
"""
nomotion = NoMotion()
NoMotion struct. It is used to create static phantoms.
# Returns
- `nomotion`: (`::NoMotion`) NoMotion struct
# Examples
```julia-repl
julia> nomotion = NoMotion()
```
"""
struct NoMotion end

Base.getindex(mv::NoMotion, p) = mv
Base.view(mv::NoMotion, p) = mv

""" Addition of NoMotions """
Base.vcat(m1::NoMotion, m2::NoMotion, Ns1, Ns2) = m1
Base.vcat(m1::MotionList, m2::NoMotion, Ns1, Ns2) = vcat(m2, m1, 0, Ns1)
function Base.vcat(m1::NoMotion, m2::MotionList{T}, Ns1, Ns2) where {T}
mv_aux = Motion{T}[]
for m in m2.motions
m_aux = copy(m)
m_aux.spins = expand(m_aux.spins, Ns2)
m_aux.spins = SpinRange(m_aux.spins.range .+ Ns1)
push!(mv_aux, m_aux)
end
return MotionList(mv_aux)
end

""" Compare two NoMotions """
Base.:(==)(m1::NoMotion, m2::NoMotion) = true
Base.:()(m1::NoMotion, m2::NoMotion) = true

function get_spin_coords(
mv::NoMotion, x::AbstractVector{T}, y::AbstractVector{T}, z::AbstractVector{T}, t
) where {T<:Real}
return x, y, z
end

"""
times = times(motion, T)
"""
times(::NoMotion, ::Type{T}) where {T<:Real} = [zero(T)]

"""
sort_motions!(nomotion)
"""
function sort_motions!(::NoMotion)
return nothing
end
2 changes: 1 addition & 1 deletion KomaMRIBase/src/motion/motionlist/Motion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -172,5 +172,5 @@ function Base.view(m::Motion, p)
end

# Auxiliary functions
times(m::Motion) = times(m.time)
times(m::Motion{T}, ::Type{T}) where {T<:Real} = times(m.time)
is_composable(m::Motion) = is_composable(m.action)
53 changes: 0 additions & 53 deletions KomaMRIBase/src/motion/nomotion/NoMotion.jl

This file was deleted.

1 change: 0 additions & 1 deletion KomaMRICore/src/simulation/Functors.jl
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,6 @@ f64(m) = paramtype(Float64, m)
# Koma motion-related adapts
adapt_storage(backend::KA.GPU, xs::MotionList) = MotionList(gpu.(xs.motions, Ref(backend)))
adapt_storage(T::Type{<:Real}, xs::MotionList) = MotionList(paramtype.(T, xs.motions))
adapt_storage(T::Type{<:Real}, xs::NoMotion) = NoMotion{T}()

#The functor macro makes it easier to call a function in all the parameters
# Phantom
Expand Down
11 changes: 6 additions & 5 deletions KomaMRIPlots/src/ui/DisplayFunctions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1035,8 +1035,8 @@ function plot_phantom_map(
max_time_samples=100,
kwargs...,
)
function interpolate_times(motion)
t = times(motion)
function interpolate_times(motion, T)
t = times(motion, T)
if length(t)>1
# Interpolate time points (as many as indicated by time_samples)
itp = interpolate((1:(time_samples + 1):(length(t) + time_samples * (length(t) - 1)), ), t, Gridded(Linear()))
Expand All @@ -1045,9 +1045,9 @@ function plot_phantom_map(
return t
end

function process_times(motion)
function process_times(motion, T)
KomaMRIBase.sort_motions!(motion)
t = interpolate_times(motion)
t = interpolate_times(motion, T)
# 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]
Expand Down Expand Up @@ -1122,7 +1122,8 @@ function plot_phantom_map(
cmin_key = get(kwargs, :cmin, factor * cmin_key)
cmax_key = get(kwargs, :cmax, factor * cmax_key)

t = process_times(obj.motion)
T = typeof(obj).parameters[1]
t = process_times(obj.motion, T)
x, y, z = get_spin_coords(obj.motion, obj.x, obj.y, obj.z, t')

x0 = -maximum(abs.([x y z])) * 1e2
Expand Down
8 changes: 1 addition & 7 deletions docs/src/reference/2-koma-base.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,19 +22,13 @@ heart_phantom

## `Motion`-related functions

### `AbstractMotion` types and related functions
```@docs
NoMotion
Motion
MotionList
get_spin_coords
```

### `Motion`

```@docs
Motion
```

### `AbstractAction` types

```@docs
Expand Down
2 changes: 1 addition & 1 deletion examples/3.tutorials/lit-05-SimpleMotion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ using PlotlyJS # hide
sys = Scanner() # hide

# It can also be interesting to see the effect of the patient's motion during an MRI scan.
# For this, Koma provides the ability to add `motion <: AbstractMotion` to the phantom.
# For this, Koma provides the ability to add `motion` to the phantom.
# In this tutorial, we will show how to add a [`Translate`](@ref) motion to a 2D brain phantom.

# First, let's load the 2D brain phantom used in the previous tutorials:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ In this excercise we will simplify this distribution, but we will obtain a simil
# (3.1) Create the new obj_t2star phantom
begin
# (3.1.1) Create an empty phantom
obj_t2star = Phantom{Float64}(x=[])
obj_t2star = Phantom()
# (3.1.2) Define the linear off-resonance distribution
Niso = 20
linear_offresonance_distribution = 2π .* range(-10, 10, Niso)
Expand Down

0 comments on commit 86a504c

Please sign in to comment.