JuliaHealth / KomaMRI.jl

Koma is a Pulseq-compatible framework to efficiently simulate Magnetic Resonance Imaging (MRI) acquisitions. The main focus of this package is to simulate general scenarios that could arise in pulse sequence development.
https://JuliaHealth.github.io/KomaMRI.jl/dev/
MIT License
109 stars 18 forks source link

New motion approach to combine SimpleMotion and ArbitraryMotion into the same phantom #410

Open pvillacorta opened 2 months ago

pvillacorta commented 2 months ago

Idea for new motion related types and structs redistribution, in order to be able to combine them and make them more alike:

abstract type Motion end

abstract type SimpleMotion <: Motion end
abstract type ArbitraryMotion <: Motion end

Specific motion structures:

Translation <: SimpleMotion
Rotation <: SimpleMotion
...
PeriodicHeartBeat <: SimpleMotion
Trajectory <: ArbitraryMotion
PeriodicTrajectory <: ArbitraryMotion
FlowTrajectory <: ArbitraryMotion
PeriodicFlowTrajectory <: ArbitraryMotion

That way, phantom struct will look like this:

@with_kw mutable struct Phantom{T<:Real}
    ...
    motion::Vector{<:Motion{T}} = []
    ...
end

and motion will be initialized, for example, like this:

motion = [Translation(...), Rotation(...), Trajectory(...), HeartBeat(...)]
cncastillo commented 2 months ago

This is a very good idea! We just need to think about how this will impact:

pvillacorta commented 1 month ago

Another possible approach regarding periodicity:

struct Periodic{T<:Real} <: TimeScale{T} period::T asymmetry::T end

struct TimeRange{T<:Real} <: TimeScale{T} t_start::T t_end::T end

- KomaMRIBase/src/datatypes/phantom/motion/simplemotion/Translation.jl:
```julia
struct Translation{T<:Real, TS<:TimeScale{T}} <: SimpleMotionType{T}
   dx::T
   dy::T
   dz::T
   times::TS
end

This way we could avoid re-defining structs and methods for periodic types. A SimpleMotiontype, for example, could be instantiated as follows: tr = Translation(dx, dy, dz, Aperiodic(t_start, t_end))

pvillacorta commented 1 month ago

Regarding flow:

function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{T}, motion<:Vector{Motion{T}})
   for m in motion
      reset_magnetization!(M, Mxy, m)
   end
   return nothing
end

function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{T}, motion<:Motion{T})
   return nothing
end

function reset_magnetization!(M::Mag{T}, Mxy::AbstractArray{T}, motion::FlowTrajectory{T})
   # do things with flow (create interpolation functions, resolve them, etc):
   Mxy .= ...
   M.z .= ...
   return nothing
end

This will be called from simulation functions:

function run_spin_precession!(...)
   ...
   #Mxy precession and relaxation, and Mz relaxation
   tp = cumsum(seq.Δt) # t' = t - t0
   dur = sum(seq.Δt)   # Total length, used for signal relaxation
   Mxy = [M.xy M.xy .* exp.(1im .* ϕ .- tp' ./ p.T2)] #This assumes Δw and T2 are constant in time
   M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1))
   reset_magnetization!(M, Mxy, phantom.motion)
   M.xy .= Mxy[:, end]
   #Acquired signal
   sig .= transpose(sum(Mxy[:, findall(seq.ADC)]; dims=1)) #<--- TODO: add coil sensitivities
return nothing