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. ]
```