*The problem* of SDCFlows (and how to fix it)

[oesteban]

What would you like to see added in this software?

Instead of directly starting from coding, I will explain what I currently think has happened.
I have come to the realization that I've unintendedly been over-complicating the implementation of SDCFlows' interface.

Let me first thank you for your patience on this front, @mgxd, @effigies, @eilidhmacnicol and others. It took too much time for me to realize.

Why first share rather than code away?

I'm under the impression that SDCFlows needs better documentation and better communication from its architect (which would be me, in this case). The notebook describing the problem was a good move, but a second notebook describing correction has been missing and has contributed to the failure to set clear goals (and very likely would have made me realize of the following). Since I don't have time to write that notebook right now, I will explain what I just realized.

Daunted by the size of the forest, I missed that I should've cared for the first tree only

Since October 2022 I've been banging my head against the same wall, trying to interpolate a B-Spline field on a grid that is not aligned with the target grid in a manageable manner (memory, time). This problem became apparent with Mathias' data because there's a rotation of the image w.r.t. the fieldmap. This seemed essential as well to implement the vision of being able to correct each volume of the EPI series with the exact fieldmap realization (when the head moves, the deformation field changes because the phase-encoding axis doesn't while inhomogeneity of the field does with the head).

Well, after some thought, that is wrong -- or better said, pretty stupid. Instead:

Preconditions:

• We have a "target EPI" - a BOLD or DWI dataset, without having gone through HMC or SDC.
• We have also the list of HMC matrices
• We have estimated the fieldmap's coefficients with SDCFlows
• We have the "fieldmap-to-EPI" affine transform (or the reverse, and we invert) that aligns the target EPI and the fieldmap's "magnitude" (phasediff et al.) or "reference" (pepolar, syn).

The algorithm:

1. Project the voxel positions of the target EPI (that you want to correct) into the fieldmap space using the "fieldmap-to-EPI" transform (that is, if you use antsApplyTransform, then the fieldmap is the moving and EPI is the reference). This is nitransforms.linear.Affine.from_file("fieldmap-to-EPI.h5").map((x, y, z)), where x, y, z are the coordinates in mm of each voxel of the EPI grid. This can be done easily at once just by rotating the affine of the EPI through the fieldmap-to-EPI transform (meaning, we just get another NIfTI file of the EPI that a decent viewer with good implementation of the orientation will render aligned with the fieldmap, despite them being on different grids). The latter is also doable with a nitransforms one-liner.
2. We interpolate (only this once) the field there (and in this case our B-Spline coefficients are aligned with the grid, so it's fast), and we have the fieldmap in Hz interpolated on the EPI grid. If we just replace the affine with the original affine of the EPI, our fieldmap is now on the target EPI without any more interpolations!
3. At this point, it would be convenient to convert the fieldmap from Hz into a voxel shift map (if I'm not wrong, this is just accounting for voxel zoom in the PE direction). (EDIT: see comment below)
4. Say you want to now resample your 60 volumes of BOLD with simultaneous correction of SD. For each volume you:
   4.1. Project all the voxels of the reference onto the volume
   4.2. Convert all those coordinates in mm into voxel coordinates
   4.3. Instead of directly interpolating the BOLD value on the projected location, now you pull up your voxel shift map, and at that voxel coordinate, you get the displacement of that voxel to account for distortion. You update the coordinates of each voxel with this amount (minus/plus) and now yes, interpolate because you'll be off-grid almost always.
   4.4. Now you have a data array on the same grid of the reference, and you can use its affine to build a new NIfTI where head motion AND susceptibility distortions are not present.

I hope the algorithm is clearly described. If not, it would be critical that you poked at me until we're all on the same page.

The good news is the above is trivial to implement with our current kit. And no one is doing this because of the fragmentation of code. I think this would be as trivial to implement for FSL as it is for us (in their case only with TOPUP, as opposed to us, where it would work for all estimation approaches). But I believe AFNI and others would not be as close.

Do you have any interest in helping implement the feature?

Yes

Additional information / screenshots

No response

[oesteban]

Re: step 3

The voxel zoom is relevant to convert to voxels to mm (reasonably). To convert from Hz (vox/s) to voxels we need the readout time (s). And we actually have it already there, exactly in the way I was mentioning above:

sdcflows/sdcflows/transform.py

Lines 197 to 200 in 891e81b

	# fmap * ro_time is the voxel-shift map (VSM)
	# The VSM is just the displacements field given in index coordinates
	# voxcoords is the deformation field, i.e., the target position of each voxel
	voxcoords[pe_axis, ...] += fmap.reshape(-1) * ro_time

[oesteban]

Okay, I will have to hold the reins a little because there's no way around reconstructing the field on a rotated grid. However, I think it would be cheap to start with an approximation and then consider if we want the very precise thing. See the above draft PR (and then new approx argument to the B0FieldTransform.fit().

[effigies]

Okay, just doing some math. Let's define the following spaces:

• Bold reference space: $\mathbf{B}$
• Bold volume space: $\mathbf{V}$
• Fieldmap space: $\mathbf{F}$

Each has an accompanying affine matrix $A_S$ that transforms an index $\vec i_S$ to a world coordinate $\vec c_S$. So:

$\vec c_B = A_B \vec i_B$
$\vec c_V = A_V \vec i_V$
$\vec c_F = A_F \vec i_F$

A transform $\mathcal{T}_{XY}$ is the affine transform from coordinates in space $\mathbf{X}$ to $\mathbf{Y}$. Assume we have invertible transforms:

$\mathcal T_{BF} : \mathbf{B} \rightarrow \mathbf{F}$
$\mathcal T_{VB} : \mathbf{V} \rightarrow \mathbf{B}$
$\vec c_F = \mathcal T_{BF} \vec c_B$
$\vec c_B = \mathcal T_{VB} \vec c_V$

and we will use the inverses when going in the opposite direction.

There is a fieldmap function $\Phi(\vec c_F)$ interpolated over a discrete field $\phi(\vec i_F)$ in Hz.

Finally, there is a constant readout field $\vec r(\vec c i_V)$ over $\mathbf V$ equal to the total readout time divided by the voxel size, constrained to the phase-encoding direction. This is most usefully recast as $\vec R_V(\vec c_v) = A_Vr(\vec i_V)$, and is constant, so we will write $\vec R_V$.

The problem is to correcting a bold volume in $\mathbf{V}$ simultaneously accounting for motion and distortion to resample in $\mathbf{B}$.

Given $\vec i_B$:

The field at that location is $\Phi(\mathcal{T}_{BF}A_B\vec i_B)$.

The uncorrected location in $\mathbf V$ is $\mathcal T_{VB}^{-1}A_B\vec i_B$. The corrected location is thus

$\mathcal T_{VB}^{-1}A_B\vec i_B + \Phi(\mathcal{T}_{BF}A_B\vec i_B)\cdot \vec R_V$

If we project $\vec R_V$ into $\mathbf{B}$ this becomes:

$\mathcal T_{VB}^{-1}(A_B\vec i_B + \Phi(\mathcal T_{BF}A_B\vec i_B)\cdot \mathcal T_{VB}\vec R_V)$

That is, at each voxel in the BOLD reference space, we interpolate the fieldmap in Hz and the per-volume readout vector in seconds to retrieve a "real coordinate" that we project into the volume space and interpolate.

[effigies]

Okay, rereading what you wrote, I think I disagree with steps 3 and 4. I think step three is to convert to a coordinate shift map:

readout_vec = np.zeros((3,), 'int8')
readout_vec["ijk".index(pedir[0])] = 1
if pedir[1:] == "-":
    readout_vec *= -1

volspace_readout = trt * (vol.affine @ readout_vec)
refspace_readout = mc_affine @ volspace_readout

coordinate_shift = refspace_readout * fieldmap

Now for step four, I simply interpolate from a shifted world location:

target_coord = boldref.affine @ target_index
volume_coord = inv(mc_affine) @ (target_coord + coordinate_shift * field[target_index])

(This is all pseudocode; we need to handle the 1s tacked onto each index/coordinate.)

[oesteban]

Yes, let's simplify further by dividing the problem into two halves.

Let's address first the most practical part of the problem: correcting the BOLD volume in $\mathbf{V}$, given a list of $\mathcal T_{BV}^t$ (where $t$ is the time index in $\mathbf{V}$) --estimated with, e.g., mcflirt--, the reference of head-motion estimation $\mathbf{B}$, and a voxel-shift map defined on $\mathbf{B}$, $\Theta_\mathbf{B}(\vec i_B)$.

$\Theta_\mathbf{B}(\vec i_B)$ is a mapping of the field in voxels, such that $\Theta_\mathbf{B}((20, 10, 6)) = 1.5$ is a displacement of 1.5 along the PE direction at that particular $(20, 10, 6)$ voxel. If the PE direction is j or j- (sign should be encoded in the mapping already), then that voxel's coordinates in the distorted $\mathbf{V}$ are $(20, 11.5, 6)$. The intent is to also clear out the affines to convert to physical space as often as the least possible.

Without distortion ($\Theta_\mathbf{B} = 0$), then the resampling problem is to find where every $\vec i_\mathbf{B}$ from the BOLD reference lands in the BOLD volume's grid ($\vec i_\mathbf{V}^t$). Then an interpolator will get those $\vec i_\mathbf{V}^t$ and calculate those off-grid values, and the interpolated values are then pulled back into a data grid on the BOLD reference space, so that the head-motion corrected $\vec i_\mathbf{B}^t$ are assigned these interpolated values:

$\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B} \vec i_\mathbf{B}$ (*)

Please note that the VOX/RAS conversions are necessary because $\mathcal T_{BV}^t$ are estimated in physical (RAS) coordinates.

Now, let's have a nonzero fieldmap $\Theta_\mathbf{B} = (0, \Delta \vec i, 0)$. My realization of today is that we just need to do:

$\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B} \vec i_\mathbf{B} + \Delta \vec i$,

which effectively is equivalent to applying $\mathcal T_{VB}^t$ to $\Theta_\mathbf{B}$ (i.e. moving the fieldmap with the head, please note how the latter BV has been inverted w.r.t. above), and then calculating the voxel shift there. Suppose we work with the field in RAS shifts/displacements. In that case, everything gets even more complicated because you also need to project the transformed displacements and discard the components where there's no deformation.

Please also note how the "traditional" correction works, where you first use topup/applytopup to create a corrected version of the dataset, and then apply head motion correction:

$\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B} (\vec i_\mathbf{B} + \Delta \vec i)$

Then the second part of the problem is as follows. When we compute the fieldmap, we do not get it with reference to the target images (or references of target images). We get them in their native (fieldmap) space.

However, I think I'm going to leave it for later because you just posted something :)

EDIT: corrected subindex in (*)

[oesteban]

Translating into nitransforms'esque:

$\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B} \vec i_\mathbf{V}$

has no direct API method, instead you can implement the physical to physical mapping:

$\vec c_\mathbf{V}^t = \mathcal T_{BV}^t \vec c_\mathbf{B}$

with:

vec_c_V_t = t_BV_xfm.map(vec_c_B)

However, the above ($\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B} \vec i_\mathbf{V}$) is internally calculated by apply (pseudocode):

volume_t_in_boldref = t_BV_xfm[t].apply(bold_volume[t])

In reality, apply will move all timepoints:

bold_in_boldref = t_BV_xfm.apply(bold_volume, reference=boldref)

Please note how pulling information through the transform in the image and mapping have I/O spaces "swapped".

[oesteban]

(if we can agree on this first half #345 (comment), let's dig into the second half of the problem)

[effigies]

Sorry, had to run some errands, but have now wrapped my head around your comment, @oesteban. The fundamental difference between our two proposals is:

Yours: $\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B} (\vec i_\mathbf{B} + \Delta \vec i_B^t)$

Mine: $\vec i_\mathbf{V}^t = A_\mathbf{V}^{-1} \mathcal T_{BV}^t (A_\mathbf{B} \vec i_\mathbf{B} + \Delta \vec c_B^t)$ where $\Delta \vec c_B^t = \Phi(\mathcal T_{BF}A_B\vec i_B)\cdot \mathcal T_{VB}^t\vec R_V$

I believe you're right that calculating $A_\mathbf{V}^{-1} \mathcal T_{BV}^t A_\mathbf{B}$ once and directly mapping indices will be more efficient, so for completeness:

$\Delta \vec i_B^t = A_B^{-1}\mathcal T_{VB}^t\vec R_V\cdot \Phi(\mathcal T_{BF}A_B\vec i_B)$

I think we agree on the math at this point; the only thing I'm not sure we agree on is what the intermediate representation of the fieldmap should be in $\mathbf{B}$. I think it should be in Hz, since by the equation above it's the multiplication of a single orientation vector by a scalar field, while if I'm reading you correctly, you were proposing having a reference voxel displacement field that needs to be rotated according to each volume.

That feels like an empirical question of which works better, so happy to move on to the second half.

[oesteban]

Thanks for initiating the conversation with the maths, in hindsight it's the only way that would've worked.

Unfortunately, it's a bit late for me. I'll try to write down the second half tomorrow.

(I'm on my phone now)

[effigies]

Okay, I decided to start writing this up into a single document with derivations that we could include in the docs, and worked myself around to your voxel displacement map:

---

We are looking to resample a volume $\mathbf V^t$ into the the BOLD reference volume
$\mathbf B$, given a transform $\mathcal T_{BV}^t$:

$$\vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B \vec i_B$$

We would like to find an additional displacement $\Delta \vec i_V^t$ that
would account for distortion induced shift:

$$ \vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B \vec i_B + \Delta \vec i_V^t $$

Stepping back from $\Delta \vec i_V^t$, let's consider the corresponding
location in world coordinates $\vec c_V^t$ and a corresponding shift
$\Delta \vec c_V^t$:

$$\Delta \vec c_V^t = \Phi_F(\mathcal{T}{BF}\mathcal T{VB}^t\vec c_V^t)\vec R_V$$

This is the fieldmap, resampled into the volume space $\mathbf V^t$, multiplied
by our constant vector field.
To get the voxel shift instead of the coordinate shift,

$$\begin{align}
\Delta \vec i_V^t &= \Phi_F(\mathcal{T}{BF}\mathcal T{VB}^t\vec c_V^t)M_V^{-1}\vec R_V \
&= \Phi_F(\mathcal{T}{BF}\mathcal T{VB}^tA_V\vec i_V^t)\tau_{ro}\vec o_V \
\end{align}$$

Using our earlier calculations, we can decompose this into the fieldmap resampled
into the reference space $\mathbf B$ and the orientation vector:

$$\begin{align}
\Delta \vec i_V^t &= \Phi_F(\mathcal{T}{BF}\mathcal T{VB}^tA_VA_V^{-1} \mathcal T_{BV}^t A_B \vec i_B)\tau_{ro}\vec o_V \
&= \Phi_F(\mathcal{T}{BF} A_B \vec i_B)\tau{ro}\vec o_V \
\end{align}$$

Hence, by resampling $\Phi_F$ once into $\mathbf B$,

$$\phi_B(\vec i_B) = \Phi_F(\mathcal{T}_{BF}A_B\vec i_B)$$

We are able to fully calculate the index-to-index mapping:

$$ \vec i_V^t = A_V^{-1} \mathcal T_{BV}^t A_B \vec i_B + \phi_B(\vec i_B)\tau_{ro}\vec o_V $$

[effigies]

I spent some time this morning to see if we could recast this final calculation as a traditional concatenation of transforms, but convinced myself that that would not work. We need to be able to create a voxel-shift-map in the target space, but operate on source indices.

However, I realized that this target space can be anything, not just the BOLD reference. If we want to do this correctly and consistently, we should resample the fieldmap to every target space that we do single-shot resampling to.

Here's the doc, btw: https://hackmd.io/@effigies/rethinking-resampling. It still assumes that our target space is the BOLD reference, but I'll aim to update it to be more general if we want to use this doc.