Magnetization issues for long sequences in flow phantoms

[cncastillo]

Feature Request

After talking to @mcgrathcm, we realized it is not clear what happens for FlowPath motions when the periodic argument in TimeCurve is true.

We talked with @pvillacorta, and we realized that some problems can be solved by adding more "key points" to the simulation, to consider the the periodic case.

Questions for @pvillacorta:

• What additional points are we adding to the simulation for Path, and FlowPath? Is this documented?

Questions for @mcgrathcm:

• Can you describe which points were not clear in the docs?

Related #516

[mcgrathcm]

For me it was generally unclear what occurs at the end of FlowPath, specifically in the case were the defined motion is not truly periodic, like for flow.

Since you are not enforcing any periodicity on the FlowPath itself, what occurs at the end seems ambiguous. Do the particles reappear at their starting point? Is their magnetization preserved or reset? Is there some smart mapping to nearest starting point? If so, what happens if there are spatial particle density changes that might interfere with a nice mapping?

Ideally, for flow simulation, there should be a step that remaps magnetization, such that if I took a snapshot of all particle before and after the end of the period, the overall magnetization distribution should looks the same, but the actual particle positions should reflect either the last or the first timepoint.

Finally, I think it is especially unclear for users who have not familiarized themselves with the phantom format, as they might assume that everything is handled under the hood, and might misconstrue results from the periodic resetting with real MR artifacts.

[pvillacorta]

Since you are not enforcing any periodicity on the FlowPath itself, what occurs at the end seems ambiguous. Do the particles reappear at their starting point? Is their magnetization preserved or reset?

Yes, particles re-appear at their starting point, and their magnetization is preserved (this is one of the things that we might have to study more deeply).

Is there some smart mapping to nearest starting point? If so, what happens if there are spatial particle density changes that might interfere with a nice mapping?

There is no smart mapping. This could be investigated, but it is not implemented yet. For now, this re-mapping is delegated to the phantom design stage (in case that you are referring to a spatial re-mapping). The only difference between Path and FlowPath is that the latter contains the spin_reset matrix

Ideally, for flow simulation, there should be a step that remaps magnetization, such that if I took a snapshot of all particle before and after the end of the period, the overall magnetization distribution should looks the same, but the actual particle positions should reflect either the last or the first timepoint.

• This is related with what is colored in yellow in @cncastillo's drawing. We had a meeting two days ago, where we set this issue as future work (short- to medium-term), since it needs a bit more of research.

What additional points are we adding to the simulation for Path, and FlowPath? Is this documented?

No additional simulation points are added for Path.
Instead, we have the following implementation of the add_jump_times! function for FlowPath actions:

function add_jump_times!(t, a::FlowPath, tc::TimeCurve) 
    jump_times = (tc.t_end - tc.t_start)/(size(a.spin_reset)[2]-1) * (getindex.(findall(a.spin_reset .== 1), 2) .- 1) .- 1e-6
    append!(t, jump_times)
end

which adds to the simulation the time points at which a spin reset is produced. The substraction of 1us (.- 1e-6) to these time points ensures that we are evaluating within the time interval in which interpolation will return 1 (related #511). This way, we do not miss any spin reset.

However, this add_jump_times! function certainly presents some problems. I will put them as checkboxes so I can address them in a PR (#638):

• Until now, I hadn't realized that, if periodic == true, no time points are added for next periods. I must extend the function with push!(jump_times, jump_times .+ n*period) or something similar.
• Since a big "jump" or spatial discontinuity may be produced when spins re-appear at their starting point at the beginning of each period, we might well ensure that no phase is accumulated during a simulation time step (Δt in the drawing) due to this position change. To achieve this, we must include in jump_times two more time points (as depicted in blue in the drawing): tf and tf+eps

[mcgrathcm]

After speaking with @cncastillo, this makes sense to me.
I do wonder if there is a way to achieve the same result without adding additional simulation points, but for now this is a nice solution.

Also, we did run into another issue (#625), that is related.

[cncastillo]

This -1us seems a little hacky I would try to use some constant like MIN_RISE_TIME instead.