Tutorial for SimpleMotion creation and simulation

KomaMRIBase/src/datatypes/Phantom.jl

@@ -78,9 +78,7 @@ Base.:(≈)(obj1::Phantom, obj2::Phantom)      = reduce(&, [getfield(obj1, field
 Base.:(==)(m1::MotionModel, m2::MotionModel) = false
 Base.:(≈)(m1::MotionModel, m2::MotionModel)  = false
 
-"""
-Separate object spins in a sub-group
-"""
+"""Separate object spins in a sub-group"""
 Base.getindex(obj::Phantom, p::Union{AbstractRange,AbstractVector,Colon}) = begin
     fields = []
     for field in Iterators.filter(x -> !(x == :name), fieldnames(Phantom))
@@ -110,9 +108,7 @@ end
     return obj1
 end
 
-"""
-dims = get_dims(obj)
-"""
+"""dims = get_dims(obj)"""
 function get_dims(obj::Phantom)
     dims = Bool[]
     push!(dims, any(x -> x != 0, obj.x))
@@ -128,7 +124,7 @@ end
 """
     obj = heart_phantom(...)
 
-Heart-like LV phantom. The variable `circumferential_strain` and `radial_strain` are for streching (if positive) 
+Heart-like LV 2D phantom. The variable `circumferential_strain` and `radial_strain` are for streching (if positive) 
 or contraction (if negative). `rotation_angle` is for rotation.
 
 # Arguments

KomaMRIBase/src/datatypes/phantom/motion/ArbitraryMotion.jl

@@ -11,9 +11,37 @@ const LinearInterpolator = Interpolations.Extrapolation{
 } where {T<:Real,V<:AbstractVector{T}}
 
 """
-Arbitrary Motion
+    motion = ArbitraryMotion(period_durations, dx, dy, dz)
 
-x = x + ux
+ArbitraryMotion model. For this motion model, it is necessary to define 
+motion for each spin independently, in x (`dx`), y (`dy`) and z (`dz`).
+`dx`, `dy` and `dz` are three matrixes, of (``N_{spins}`` x ``N_{discrete\\,times}``) each.
+This means that each row corresponds to a spin trajectory over a set of discrete time instants.
+`period_durations` is a vector that contains the period for periodic (one element) or 
+pseudo-periodic (two or more elements) motion.
+The discrete time instants are calculated diving `period_durations` by ``N_{discrete\\,times}``.
+
+This motion model is useful for defining arbitrarly complex motion, specially
+for importing the spin trajectories from another source, like XCAT or a CFD.
+
+# Arguments
+- `period_durations`: (`Vector{T}`) 
+- `dx`: (`::Array{T,2}`) matrix for displacements in x
+- `dy`: (`::Array{T,2}`) matrix for displacements in y
+- `dz`: (`::Array{T,2}`) matrix for displacements in z
+
+# Returns
+- `motion`: (`::ArbitraryMotion`) ArbitraryMotion struct
+
+# Examples
+```julia-repl
+julia> motion = ArbitraryMotion(
+            [1.0], 
+            0.01.*rand(1000, 10), 
+            0.01.*rand(1000, 10), 
+            0.01.*rand(1000, 10)
+        )
+```
 """
 struct ArbitraryMotion{T<:Real,V<:AbstractVector{T}} <: MotionModel{T}
     period_durations::Vector{T}
@@ -27,27 +55,27 @@ end
 
 function ArbitraryMotion(
     period_durations::AbstractVector{T},
-    Δx::AbstractArray{T,2},
-    Δy::AbstractArray{T,2},
-    Δz::AbstractArray{T,2},
+    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(Δx)[1]
-    num_pieces = size(Δx)[2] + 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),Δx),1,length(period_durations)),zeros(Ns,1))
-    Δ[:,:,2] = hcat(repeat(hcat(zeros(Ns,1),Δy),1,length(period_durations)),zeros(Ns,1))
-    Δ[:,:,3] = hcat(repeat(hcat(zeros(Ns,1),Δz),1,length(period_durations)),zeros(Ns,1))
+    Δ[:,:,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, Δx, Δy, Δz, etpx, etpy, etpz)
+    return ArbitraryMotion(period_durations, dx, dy, dz, etpx, etpy, etpz)
 end
 
 function Base.getindex(

KomaMRIBase/src/datatypes/phantom/motion/SimpleMotion.jl

@@ -6,13 +6,25 @@ is_composable(motion_type::SimpleMotionType{T}) where {T<:Real} = false
 """
     motion = SimpleMotion(types)
 
-SimpleMotion model
+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
+as simple motions we want to combine.
 
 # Arguments
 - `types`: (`::Vector{<:SimpleMotionType{T}}`) vector of simple motion types
 
 # Returns
 - `motion`: (`::SimpleMotion`) SimpleMotion struct
+
+# Examples
+```julia-repl
+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}}
@@ -39,9 +51,10 @@ Base.:(==)(t1::SimpleMotionType, t2::SimpleMotionType) = false
 Base.:(≈)(t1::SimpleMotionType, t2::SimpleMotionType)  = false
 
 """
-    x, y, z = het_spin_coords(motion, x, y, z, t')
+    x, y, z = get_spin_coords(motion, x, y, z, t)
 
 Calculates the position of each spin at a set of arbitrary time instants, i.e. the time steps of the simulation. 
+For each dimension (x, y, z), the output matrix has ``N_{\text{spins}}`` rows and `length(t)` columns.
 
 # Arguments
 - `motion`: (`::MotionModel`) phantom motion
@@ -51,7 +64,7 @@ Calculates the position of each spin at a set of arbitrary time instants, i.e. t
 - `t`: (`::AbstractArray{T<:Real}`) horizontal array of time instants
 
 # Returns
-- `z, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time
+- `x, y, z`: (`::Tuple{AbstractArray, AbstractArray, AbstractArray}`) spin positions over time
 """
 function get_spin_coords(
     motion::SimpleMotion{T},

KomaMRIBase/src/datatypes/phantom/motion/simplemotion/HeartBeat.jl

@@ -5,14 +5,19 @@ HeartBeat struct. It produces a heartbeat-like motion, characterised by three ty
 Circumferential, Radial and Longitudinal
 
 # Arguments
-- `circumferential_strain`: (`::Real`, `=-0.3`) contraction parameter
-- `radial_strain`: (`::Real`, `=-0.3`) contraction parameter
-- `longitudinal_strain`: (`::Real`, `=1`) contraction parameter
+- `circumferential_strain`: (`::Real`) contraction parameter
+- `radial_strain`: (`::Real`) contraction parameter
+- `longitudinal_strain`: (`::Real`) contraction parameter
 - `t_start`: (`::Real`, `[s]`) initial time 
 - `t_end`: (`::Real`, `[s]`) final time 
 
 # Returns
 - `heartbeat`: (`::HeartBeat`) HeartBeat struct
+
+# Examples
+```julia-repl
+julia> hb = HeartBeat(circumferential_strain=-0.3, radial_strain=-0.2, longitudinal_strain=0.0, t_start=0.2, t_end=0.5)
+```
 """
 @with_kw struct HeartBeat{T<:Real} <: SimpleMotionType{T}
     circumferential_strain :: T

KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicHeartBeat.jl

@@ -5,9 +5,9 @@ HeartBeat struct. It produces a heartbeat-like motion, characterised by three ty
 Circumferential, Radial and Longitudinal
 
 # Arguments
-- `circumferential_strain`: (`::Real`, `=-0.3`) contraction parameter
-- `radial_strain`: (`::Real`, `=-0.3`) contraction parameter
-- `longitudinal_strain`: (`::Real`, `=1`) contraction parameter
+- `circumferential_strain`: (`::Real`) contraction parameter
+- `radial_strain`: (`::Real`) contraction parameter
+- `longitudinal_strain`: (`::Real`) contraction parameter
 - `period`: (`::Real`, `[s]`) period 
 - `asymmetry`: (`::Real`)  asymmetry factor, between 0 and 1
 

KomaMRIBase/src/datatypes/phantom/motion/simplemotion/PeriodicRotation.jl

@@ -68,7 +68,9 @@ function displacement_z(
     α = t_unit .* motion_type.pitch
     β = t_unit .* motion_type.roll
     γ = t_unit .* motion_type.yaw
-    return -sind.(β) .* x + cosd.(β) .* sind.(α) .* y .- z
+    return -sind.(β) .* x + 
+            cosd.(β) .* sind.(α) .* y +
+            cosd.(β) .* cosd.(α) .* z .- z
 end
 
 function times(motion_type::PeriodicRotation)

KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Rotation.jl

@@ -2,55 +2,39 @@
     rotation = Rotation(t_start, t_end, pitch, roll, yaw)
 
 Rotation motion struct. It produces a rotation of the phantom in the three axes: 
-x (pitch), y (roll), and z (yaw)
+x (pitch), y (roll), and z (yaw).
+We follow the RAS (Right-Anterior-Superior) orientation, 
+and the rotations are applied following the right-hand rule (counter-clockwise):
 
-```math
-\begin{align*}
-ux  &= cos(\gamma)cos(\beta)x \\ 
-    &+ (cos(\gamma)sin(\beta)sin(\alpha) - sin(\gamma)cos(\alpha))y\\ 
-    &+ (cos(\gamma)sin(\beta)cos(\alpha) + sin(\gamma)sin(\alpha))z\\
-    &- x
-\end{align*}
-```
-
-```math
-\begin{align*}
-uy  &= sin(\gamma)cos(\beta)x \\ 
-    &+ (sin(\gamma)sin(\beta)sin(\alpha) + cos(\gamma)cos(\alpha))y\\ 
-    &+ (sin(\gamma)sin(\beta)cos(\alpha) - cos(\gamma)sin(\alpha))z\\
-    &- y
-\end{align*}
-```
+![Head Rotation Axis](../assets/head-rotation-axis.svg)
 
+The applied rotation matrix is obtained as follows: 
 ```math
-\begin{align*}
-uz  &= -sin(\beta)x \\ 
-    &+ cos(\beta)sin(\alpha)y\\ 
-    &+ cos(\beta)cos(\alpha)z\\
-    &- z
-\end{align*}
-```
-
-where:
-
-```math
-\alpha  = \left\{\begin{matrix}
-0, & t <= t_start \\
-\frac{pitch}{t_end-t_start}(t-t_start), & t_start < t < t_end \\ 
-pitch, & t >= t_end
-\end{matrix}\right.
-,\qquad
-\beta  = \left\{\begin{matrix}
-0, & t <= t_start \\
-\frac{roll}{t_end-t_start}(t-t_start), & t_start < t < t_end \\ 
-roll, & t >= t_end
-\end{matrix}\right.
-,\qquad
-\gamma  = \left\{\begin{matrix}
-0, & t <= t_start \\
-\frac{yaw}{t_end-t_start}(t-t_start), & t_start < t < t_end \\ 
-yaw, & t >= t_end
-\end{matrix}\right.
+\begin{equation}
+\begin{aligned}
+R &= R_z(\alpha) R_y(\beta) R_x(\gamma) \\
+  &= \begin{bmatrix}
+\cos \alpha & -\sin \alpha & 0 \\
+\sin \alpha & \cos \alpha & 0 \\
+0 & 0 & 1
+\end{bmatrix}
+\begin{bmatrix}
+\cos \beta & 0 & \sin \beta \\
+0 & 1 & 0 \\
+-\sin \beta & 0 & \cos \beta
+\end{bmatrix}
+\begin{bmatrix}
+1 & 0 & 0 \\
+0 & \cos \gamma & -\sin \gamma \\
+0 & \sin \gamma & \cos \gamma
+\end{bmatrix} \\
+  &= \begin{bmatrix}
+\cos \alpha \cos \beta & \cos \alpha \sin \beta \sin \gamma - \sin \alpha \cos \gamma & \cos \alpha \sin \beta \cos \gamma + \sin \alpha \sin \gamma \\
+\sin \alpha \cos \beta & \sin \alpha \sin \beta \sin \gamma + \cos \alpha \cos \gamma & \sin \alpha \sin \beta \cos \gamma - \cos \alpha \sin \gamma \\
+-\sin \beta & \cos \beta \sin \gamma & \cos \beta \cos \gamma
+\end{bmatrix}
+\end{aligned}
+\end{equation}
 ```
 
 # Arguments
@@ -63,6 +47,10 @@ yaw, & t >= t_end
 # Returns
 - `rotation`: (`::Rotation`) Rotation struct
 
+# Examples
+```julia-repl
+julia> rt = Rotation(pitch=15.0, roll=0.0, yaw=20.0, t_start=0.1, t_end=0.5)
+```
 """
 @with_kw struct Rotation{T<:Real} <: SimpleMotionType{T}
     pitch      :: T
@@ -83,12 +71,12 @@ function displacement_x(
     t::AbstractArray{T},
 ) where {T<:Real}
     t_unit = unit_time(t, motion_type.t_start, motion_type.t_end)
-    α = t_unit .* (motion_type.pitch)
+    α = t_unit .* (motion_type.yaw)
     β = 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
+    γ = t_unit .* (motion_type.pitch)
+    return cosd.(α) .* cosd.(β) .* x +
+           (cosd.(α) .* sind.(β) .* sind.(γ) .- sind.(α) .* cosd.(γ)) .* y +
+           (cosd.(α) .* sind.(β) .* cosd.(γ) .+ sind.(α) .* sind.(γ)) .* z .- x
 end
 
 function displacement_y(
@@ -99,12 +87,12 @@ function displacement_y(
     t::AbstractArray{T},
 ) where {T<:Real}
     t_unit = unit_time(t, motion_type.t_start, motion_type.t_end)
-    α = 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
+    α = 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
 end
 
 function displacement_z(
@@ -115,10 +103,12 @@ function displacement_z(
     t::AbstractArray{T},
 ) where {T<:Real}
     t_unit = unit_time(t, motion_type.t_start, motion_type.t_end)
-    α = t_unit .* motion_type.pitch
-    β = t_unit .* motion_type.roll
-    γ = t_unit .* motion_type.yaw
-    return -sind.(β) .* x + cosd.(β) .* sind.(α) .* y .- z
+    α = t_unit .* (motion_type.yaw)
+    β = t_unit .* (motion_type.roll)
+    γ = t_unit .* (motion_type.pitch)
+    return -sind.(β) .* x + 
+            cosd.(β) .* sind.(γ) .* y +
+            cosd.(β) .* cosd.(γ) .* z .- z
 end
 
 times(motion_type::Rotation) = begin

KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl

@@ -1,31 +1,9 @@
 @doc raw"""
     translation = Translation(t_start, t_end, dx, dy, dz)
 
-Translation motion struct. It produces a translation of the phantom in the three directions x, y and z.
-
-```math
-ux=\left\{\begin{matrix}
-0, & t <= t_start\\
-\frac{dx}{t_end-t_start}(t-t_start), & t_start < t < t_end\\ 
-dx, & t >= t_end
-\end{matrix}\right.
-```
-
-```math
-uy=\left\{\begin{matrix}
-0, & t <= t_start\\
-\frac{dy}{t_end-t_start}(t-t_start), & t_start < t < t_end\\ 
-dy, & t >= t_end
-\end{matrix}\right.
-```
-
-```math
-uz=\left\{\begin{matrix}
-0, & t <= t_start\\
-\frac{dz}{t_end-t_start}(t-t_start), & t_start < t < t_end\\ 
-dz, & t >= t_end
-\end{matrix}\right.
-```
+Translation motion struct. It produces a linear translation of the phantom.
+Its fields are the final displacements in the three axes (dx, dy, dz) 
+and the start and end times of the translation.
 
 # Arguments
 - `dx`: (`::Real`, `[m]`) translation in x
@@ -37,6 +15,10 @@ dz, & t >= t_end
 # Returns
 - `translation`: (`::Translation`) Translation struct
 
+# Examples
+```julia-repl
+julia> tr = Translation(dx=0.01, dy=0.02, dz=0.03, t_start=0.0, t_end=0.5)
+```
 """
 @with_kw struct Translation{T<:Real} <: SimpleMotionType{T}
     dx         :: T