Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix Phantom-related bugs (#498 and more) #499

Merged
merged 8 commits into from
Nov 3, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
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/MotionList.jl")
include("motion/nomotion/NoMotion.jl")
# Phantom
include("datatypes/Phantom.jl")
# Simulator
Expand Down
8 changes: 4 additions & 4 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 All @@ -31,7 +31,7 @@ julia> obj.ρ
"""
@with_kw mutable struct Phantom{T<:Real}
name::String = "spins"
x :: AbstractVector{T}
x::AbstractVector{T} = @isdefined(T) ? T[] : Float64[]
y::AbstractVector{T} = zeros(eltype(x), size(x))
z::AbstractVector{T} = zeros(eltype(x), size(x))
ρ::AbstractVector{T} = ones(eltype(x), size(x))
Expand All @@ -47,7 +47,7 @@ julia> obj.ρ
Dθ::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 Expand Up @@ -115,7 +115,7 @@ function get_dims(obj::Phantom)
push!(dims, any(x -> x != 0, obj.x))
push!(dims, any(x -> x != 0, obj.y))
push!(dims, any(x -> x != 0, obj.z))
return dims
return sum(dims) > 0 ? dims : Bool[1, 0, 0]
end

"""
Expand Down
11 changes: 0 additions & 11 deletions KomaMRIBase/src/motion/AbstractMotion.jl

This file was deleted.

14 changes: 14 additions & 0 deletions KomaMRIBase/src/motion/motionlist/Motion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,20 @@ julia> motion = Motion(
spins ::AbstractSpinSpan = AllSpins()
end

# Main constructors
function Motion(action)
T = first(typeof(action).parameters)
return Motion(action, TimeRange(zero(T)), AllSpins())
end
function Motion(action, time::AbstractTimeSpan)
T = first(typeof(action).parameters)
return Motion(action, time, AllSpins())
end
function Motion(action, spins::AbstractSpinSpan)
T = first(typeof(action).parameters)
return Motion(action, TimeRange(zero(T)), spins)
end

# Custom constructors
"""
translate = Translate(dx, dy, dz, time, spins)
Expand Down
21 changes: 13 additions & 8 deletions KomaMRIBase/src/motion/motionlist/MotionList.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,8 @@
include("Action.jl")
include("SpinSpan.jl")
include("TimeSpan.jl")
include("Motion.jl")

"""
motionlist = MotionList(motions...)

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

""" Constructors """
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}()`"
MotionList(motions::Motion...) = 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}()`"

""" MotionList sub-group """
function Base.getindex(mv::MotionList{T}, p) where {T<:Real}
motion_array_aux = Motion{T}[]
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 @@ -135,20 +140,20 @@ end
"""
times = times(motion)
"""
function times(ml::MotionList{T}) where {T<:Real}
function times(ml::MotionList)
nodes = reduce(vcat, [times(m) 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`
Expand Down
30 changes: 9 additions & 21 deletions KomaMRIBase/src/motion/nomotion/NoMotion.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
nomotion = NoMotion{T<:Real}()
nomotion = NoMotion()

NoMotion struct. It is used to create static phantoms.

Expand All @@ -8,18 +8,18 @@ NoMotion struct. It is used to create static phantoms.

# Examples
```julia-repl
julia> nomotion = NoMotion{Float64}()
julia> nomotion = NoMotion()
```
"""
struct NoMotion{T<:Real} <: AbstractMotion{T} end
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, m2::NoMotion, Ns1, Ns2) = vcat(m2, m1, 0, Ns1)
function Base.vcat(m1::NoMotion{T}, m2, Ns1, Ns2) where {T}
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)
Expand All @@ -31,23 +31,11 @@ function Base.vcat(m1::NoMotion{T}, m2, Ns1, Ns2) where {T}
end

""" Compare two NoMotions """
Base.:(==)(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = true
Base.:(≈)(m1::NoMotion{T}, m2::NoMotion{T}) where {T<:Real} = true
Base.:(==)(m1::NoMotion, m2::NoMotion) = true
Base.:(≈)(m1::NoMotion, m2::NoMotion) = true

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

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

"""
sort_motions!(motionset)
"""
function sort_motions!(m::NoMotion)
return nothing
end
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
35 changes: 23 additions & 12 deletions KomaMRIPlots/src/ui/DisplayFunctions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1035,19 +1035,18 @@ function plot_phantom_map(
max_time_samples=100,
kwargs...,
)
function interpolate_times(motion)
function process_times(::NoMotion)
return [zero(eltype(obj.x))]
end

function process_times(motion::MotionList)
KomaMRIBase.sort_motions!(motion)
t = times(motion)
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()))
t = itp.(1:(length(t) + time_samples * (length(t) - 1)))
end
return t
end

function process_times(motion)
KomaMRIBase.sort_motions!(motion)
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]
Expand Down Expand Up @@ -1135,9 +1134,21 @@ function plot_phantom_map(
l = Layout(;title=obj.name*": "*string(key))

if view_2d # 2D
function get_displayed_dims(v)
if sum(v) == 1
idx = argmax(v)[1]
return [idx in [1, 3], idx in [1, 2], idx in [2,3]]
else
return v
end
end

# Layout config
dims = get_displayed_dims(KomaMRIBase.get_dims(obj))
axis = ["x", "y", "z"][dims]

l[:xaxis] = attr(
title="x",
title=axis[1],
range=[x0, xf],
ticksuffix=" cm",
backgroundcolor=plot_bgcolor,
Expand All @@ -1146,7 +1157,7 @@ function plot_phantom_map(
scaleanchor="y"
)
l[:yaxis] = attr(
title="y",
title=axis[2],
range=[x0, xf],
ticksuffix=" cm",
backgroundcolor=plot_bgcolor,
Expand All @@ -1158,8 +1169,8 @@ function plot_phantom_map(
# Add traces
for i in 1:length(t)
push!(traces, scattergl(
x=(x[:,i])*1e2,
y=(y[:,i])*1e2,
x=dims[1] ? (x[:,i])*1e2 : (y[:,i])*1e2,
y=dims[1] & dims[2] ? (y[:,i])*1e2 : (z[:,i])*1e2,
mode="markers",
marker=attr(color=getproperty(obj,key)*factor,
showscale=colorbar,
Expand Down Expand Up @@ -1244,7 +1255,7 @@ function plot_phantom_map(
)
for (i, t0) in enumerate(t)
],
currentvalue_prefix="x = ",
currentvalue_prefix="t = ",
currentvalue_suffix="ms",
)]
l[:margin] = attr(t=50, l=0, r=0)
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
Loading