Confused about distortion correction with sdcflows vs fsl's `applytopup`

[psadil]

I’m struggling to understand how sdcflows is applying a fieldmap in toy example, and how this compares to FSL. Help would be appreciated. In this dummy example, the image has only 5 voxels, 3 with signal. Before distortion, all voxels with signal have intensity 1, but distortion sends the third voxel’s signal to the second.

# /// script
# requires-python = ">=3.11"
# dependencies = [
#     "numpy",
#     "sdcflows",
# ]
#
# [tool.uv.sources]
# sdcflows = { git = "https://github.com/nipreps/sdcflows.git", rev = "656a51d4a734d2d8d66550e9e317386f677801a3" }
# ///
import numpy as np

import nibabel as nb
import nitransforms as nt
from sdcflows import transform


# the third voxel has been shifted 1 to the left
# so, there is signal "pileup" that would be nice to address
img = [0.0, 2.0, 0.0, 1.0, 0.0]

sdcflow’s result without the Jacobian seems straightforward

kwargs = {
    "data": img,
    "coordinates": (np.asarray([[0, 1, 2, 3, 4]], dtype=float)),
    "pe_info": (0, -1),
    "fmap_hz": np.asarray([0, 0, 1, 0, 0]),
    "hmc_xfm": None,
    "order": 1,
}


without_jacobian = transform._sdc_unwarp(jacobian=False, **kwargs)
print(without_jacobian)
# [0. 2. 2. 1. 0.]

For that output

• voxel 0,1: unchanged
• voxel 2: warped signal + the (warped) signal from voxel 1
• voxel 3,4: unchanged

With the Jacobian

with_jacobian = transform._sdc_unwarp(jacobian=True, **kwargs)
print(with_jacobian)
# [0.  1.  2.  1.5 0. ]

I understand that it’s the above multiplied by the jacobian of the coordinate transformation from distorted to undistorted space, but it’s just a bit surprising to me that voxel 4 has been modified, and also that the overall intensity of the image has increased. I guess the point is that the Jacobian only accounts for the remapping of the coordinates, and isn’t really there to handle signal (de)compression? And maybe modification of the 4th voxel should be understood as an artifact of the voxel gridding? Or is the lack of recover when the jacobian is used also just an artifact of the gridding?

But I’m very confused about why this is so different from FSL’s applytopup. Probably I’m inputting the data wrong?

# /// script
# requires-python = ">=3.11"
# dependencies = [
#     "nibabel",
#     "numpy",
# ]
# ///
import tempfile
from pathlib import Path
import subprocess
import numpy as np
import nibabel as nb

with tempfile.TemporaryDirectory() as tmpd:
    with tempfile.NamedTemporaryFile(suffix=".nii", dir=tmpd) as shifted:
        with tempfile.NamedTemporaryFile(suffix=".nii", dir=tmpd) as fmap:
            with tempfile.NamedTemporaryFile(dir=tmpd) as acqp:
                nb.nifti1.Nifti1Image(
                    np.asarray([0, 2, 0, 1, 0], float), affine=np.eye(4)
                ).to_filename(shifted.name)

                nb.nifti1.Nifti1Image(
                    np.asarray([0, 0, 1, 0, 0], dtype=float),
                    affine=np.eye(4),
                ).to_filename(fmap.name)
                Path(acqp.name).write_text("-1 0 0 1")

                subprocess.run(
                    [
                        "applytopup",
                        "--method=jac",
                        "--interp=trilinear",
                        f"--datain={acqp.name.removesuffix('.nii')}",
                        f"--imain={shifted.name.removesuffix('.nii')}",
                        f"--topup={fmap.name.removesuffix('.nii')}",
                        f"--out={tmpd}/clean",
                        "--inindex=1",
                    ]
                )
                out = nb.loadsave.load(Path(tmpd) / "clean.nii.gz").get_fdata()
                print(out)
# [0.         3.49996018 1.00007975 0.24970108 0.        ]

[psadil]

Oh, well, the FSL differences seem easy enough; I forgot that it takes coefficients and not the estimated fmap

[effigies]

A voxel shift map tells you where the voxel center is shifting to, since you can only have one point for data, but for the Jacobian, you care about where the voxel edges are, and those you can only get by looking at your neighbors.

Let’s look at your VSM and its Jacobian determinant:

>>> 1 + np.gradient([0, 0, -1, 0, 0])
array([1. , 0.5, 1. , 1.5, 1. ])

So what’s this telling us? Because the mean of the determinant array is 1, the size of the final field of view matches the original. Visually:

|----|----|----|----|----|  Original grid
|    |   /    /     |    |
|----|--|----|------|----|  "Actual" grid

The edge voxels are not shifted, as demanded by your VSM, and unless we really want to get fancy, our best estimate for the change in the size of a given voxel is by looking at its shift and those to the right and left of it. The central difference of [0, 0, -1] is (0 - 0)/2 + (-1 - 0)/2 = -0.5, so we’ve shrunk by half, and our resulting size is 0.5, and we conclude that voxel two’s intensity is about double what it should be.

Doing the same thing for voxel three, we have (-1 - 0)/2 + (0 - -1)/2 = 0, so both bounds of the voxel are shifted by a half, and we therefore think that the intensity that we get will be about right. And for voxel four, (0 - -1)/2 + (0 - 0)/2 = 0.5, so the voxel is growing by a half.

What I think you’re hoping for is to integrate the values in each voxel:

|____|¯¯¯¯|____|----|____|  Original grid with values (_ = 0, - = 1, ¯ = 2)

|    |   /    /     |    |  Boundary shifts

|____|¯¯|¯¯__|__----|____|  Contents resampled and joined

|    |   \    \     |    |  Regrid

|____|----|----|----|____|  Contents averaged to new values

But resampling is not actually an integration process. You’re taking a coordinate of your target voxel, and finding its location in a new image:

|____|¯¯¯¯|____|----|____|  Original grid with values (_ = 0, - = 1, ¯ = 2)

  🠛    🠛 ↳→→→↴   🠛    🠛      Voxel shifts

|____|¯¯¯¯|¯¯¯¯|----|____|  Values sampled into new locations

When you resample, you usually use an interpolating filter, such as a spline filter, which ensures that you have values at off-grid locations, and which incorporate neighgboring values and the averaging acts as a kind of integration.

However, because your VSM is entirely whole numbers, there’s nothing for an interpolating function to do. You’re not asking for values halfway between voxels 2 and 3; you’re asking for values at 1, 2, 2, 4, and 5. Thus, your resampling drops a 0 voxel and increases the average intensity of your image, and your Jacobian can only compensate a bit.

I suspect there’s some kind of Nyquist argument about how quickly a VSM map change per voxel without introducing discontinuities, but I’m not prepared to think through that at the moment.

[psadil]

This makes sense – thanks for the detailed response!

I’m now thinking about the shorthand that perhaps distortion correction is simply about addressing where voxel signals need to move (resampling), and not really about how much signal needs to move (what I think you’re referring to as integrating). In the case of distortion sending all of the third voxel’s signal to the second, I’d believe that there isn’t any information about how much to send back within the distortion map (no information about a 50/50 split vs any other ratio).

The point about dropping a whole (zero-intensity) voxel is interesting. I had thought of the third voxels’ value without the Jacobian (2. from [0. 2. 2. 1. 0.] ) as an addition (involving a zero intensity voxel), but I see that was a wishful longing for integration when this is just resampling.

I follow the argument that in some cases the sharing of signal by interpolation provides information about how much signal should end up in each voxel (acting as a kind of integration), but I admit to still feeling unsettled by being unable to recover the original signal in this toy example.