Merge pull request #499 from JuliaHealth/fix-empty-phantom

Fix Phantom-related bugs (#498 and more)

KomaMRIBase/src/KomaMRIBase.jl

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

KomaMRIBase/src/datatypes/Phantom.jl

@@ -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
@@ -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))
@@ -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))
@@ -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
 
 """

KomaMRIBase/src/motion/motionlist/Motion.jl

@@ -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)

KomaMRIBase/src/motion/motionlist/MotionList.jl

@@ -1,3 +1,8 @@
+include("Action.jl")
+include("SpinSpan.jl")
+include("TimeSpan.jl")
+include("Motion.jl")
+
 """
     motionlist = MotionList(motions...)
 
@@ -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 """
@@ -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
@@ -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`

KomaMRIBase/src/motion/nomotion/NoMotion.jl

@@ -1,5 +1,5 @@
 """
-    nomotion = NoMotion{T<:Real}()
+    nomotion = NoMotion()
 
 NoMotion struct. It is used to create static phantoms.
 
@@ -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)
@@ -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

KomaMRICore/src/simulation/Functors.jl

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

KomaMRIPlots/src/ui/DisplayFunctions.jl

@@ -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]
@@ -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,
@@ -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,
@@ -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,
@@ -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)

docs/src/reference/2-koma-base.md

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

examples/3.tutorials/lit-05-SimpleMotion.jl

@@ -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:

examples/4.reproducible_notebooks/pluto-01-gradient-echo-spin-echo.jl

@@ -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)