Skip to content

Commit

Permalink
fix+enh: wrong transform creating NIfTI object + move to transform po…
Browse files Browse the repository at this point in the history
…ints convention

The wrong affine was assigned to the NIfTI object carrying the resampled
fieldmap, which did not fail if no transform between fieldmap and EPI
was given.

This PR also changes the naming of variables from BIDS' ``mode-image``
convention to ``mode-points`` convention, as the latter is a definition
closer to what we are actually doing (transforming points).
  • Loading branch information
oesteban committed Aug 8, 2023
1 parent 01ef011 commit 53663d4
Show file tree
Hide file tree
Showing 3 changed files with 30 additions and 30 deletions.
24 changes: 12 additions & 12 deletions sdcflows/interfaces/bspline.py
Original file line number Diff line number Diff line change
Expand Up @@ -300,12 +300,12 @@ class _ApplyCoeffsFieldInputSpec(BaseInterfaceInputSpec):
)
fmap2data_xfm = InputMultiObject(
File(exists=True),
desc="the transform by which the fieldmap can be resampled on the target EPI's grid.",
desc="the transform by which the target EPI can be resampled on the fieldmap's grid.",
xor="data2fmap_xfm",
)
data2fmap_xfm = InputMultiObject(
File(exists=True),
desc="the transform by which the target EPI can be resampled on the fieldmap's grid.",
desc="the transform by which the fieldmap can be resampled on the target EPI's grid.",
xor="fmap2data_xfm",
)
in_xfms = traits.List(
Expand Down Expand Up @@ -375,19 +375,19 @@ class ApplyCoeffsField(SimpleInterface):
def _run_interface(self, runtime):
from sdcflows.transform import B0FieldTransform

fmap2data_xfm = None
data2fmap_xfm = None

if isdefined(self.inputs.fmap2data_xfm):
fmap2data_xfm = Affine.from_filename(
self.inputs.fmap2data_xfm if not isinstance(self.inputs.fmap2data_xfm, list)
else self.inputs.fmap2data_xfm[0],
if isdefined(self.inputs.data2fmap_xfm):
data2fmap_xfm = Affine.from_filename(
self.inputs.data2fmap_xfm if not isinstance(self.inputs.data2fmap_xfm, list)
else self.inputs.data2fmap_xfm[0],
fmt="itk"
).matrix
elif isdefined(self.inputs.data2fmap_xfm):
elif isdefined(self.inputs.fmap2data_xfm):
# Same, but inverting direction
fmap2data_xfm = (~Affine.from_filename(
self.inputs.data2fmap_xfm if not isinstance(self.inputs.data2fmap_xfm, list)
else self.inputs.data2fmap_xfm[0],
data2fmap_xfm = (~Affine.from_filename(
self.inputs.fmap2data_xfm if not isinstance(self.inputs.fmap2data_xfm, list)
else self.inputs.fmap2data_xfm[0],
fmt="itk"
)).matrix

Expand Down Expand Up @@ -416,7 +416,7 @@ def _run_interface(self, runtime):
self.inputs.pe_dir,
self.inputs.ro_time,
xfms=self.inputs.in_xfms if isdefined(self.inputs.in_xfms) else None,
fmap2data_xfm=fmap2data_xfm,
xfm_data2fmap=data2fmap_xfm,
approx=self.inputs.approx,
num_threads=(
self.inputs.num_threads if isdefined(self.inputs.num_threads) else None
Expand Down
34 changes: 17 additions & 17 deletions sdcflows/transform.py
Original file line number Diff line number Diff line change
Expand Up @@ -236,7 +236,7 @@ class B0FieldTransform:
def fit(
self,
target_reference: nb.spatialimages.SpatialImage,
fmap2data_xfm: np.ndarray = None,
xfm_data2fmap: np.ndarray = None,
approx: bool = True,
) -> bool:
r"""
Expand All @@ -251,15 +251,15 @@ def fit(
The image object containing a reference grid (same as that of the data
to be resampled). If a 4D dataset is provided, then the fourth dimension
will be dropped.
fmap2data_xfm : :obj:`numpy.ndarray`
xfm_data2fmap : :obj:`numpy.ndarray`
Transform that maps coordinates on the `target_reference` onto the
fieldmap reference (that is, the linear transform through which the fieldmap can
be resampled in register with the `target_reference`).
In other words, `fmap2data_xfm` is the result of calling a registration tool
In other words, `xfm_data2fmap` is the result of calling a registration tool
such as ANTs configured for a linear transform with at most 12 degrees of freedom
and called with the image carrying `target_affine` as reference and the fieldmap
reference as moving.
The result of such a registration framework is an affine (our `fmap2data_xfm` here)
The result of such a registration framework is an affine (our `xfm_data2fmap` here)
that maps coordinates in reference (target) RAS onto the fieldmap RAS.
approx : :obj:`bool`
If ``True``, do not reconstruct the B-Spline field directly on the target
Expand All @@ -277,14 +277,14 @@ def fit(
if isinstance(target_reference, (str, bytes, Path)):
target_reference = nb.load(target_reference)

approx &= fmap2data_xfm is not None # Approximate iff fmap2data_xfm is defined
fmap2data_xfm = fmap2data_xfm if fmap2data_xfm is not None else np.eye(4)
approx &= xfm_data2fmap is not None # Approximate iff xfm_data2fmap is defined
xfm_data2fmap = xfm_data2fmap if xfm_data2fmap is not None else np.eye(4)
target_affine = target_reference.affine.copy()

# Project the reference's grid onto the fieldmap's
target_reference = target_reference.__class__(
target_reference.dataobj,
fmap2data_xfm @ target_affine,
xfm_data2fmap @ target_affine,
target_reference.header,
)

Expand Down Expand Up @@ -347,7 +347,7 @@ def fit(
hdr["cal_min"] = -hdr["cal_max"]

# Cache
self.mapped = nb.Nifti1Image(fmap, target_affine, hdr)
self.mapped = nb.Nifti1Image(fmap, target_reference.affine, hdr)

if approx:
from nitransforms.linear import Affine
Expand All @@ -364,7 +364,7 @@ def apply(
pe_dir: Union[str, Sequence[str]],
ro_time: Union[float, Sequence[float]],
xfms: Sequence[np.ndarray] = None,
fmap2data_xfm: np.ndarray = None,
xfm_data2fmap: np.ndarray = None,
approx: bool = True,
order: int = 3,
mode: str = "constant",
Expand Down Expand Up @@ -394,15 +394,15 @@ def apply(
Therefore, each of these matrices express the transform of every
voxel's RAS (physical) coordinates in the image used as reference
for realignment into the coordinates of each of the EPI series volume.
fmap2data_xfm : :obj:`numpy.ndarray`
Transform that maps coordinates on the `target_reference` onto the
xfm_data2fmap : :obj:`numpy.ndarray`
Transform that maps coordinates on the ``target_reference`` onto the
fieldmap reference (that is, the linear transform through which the fieldmap can
be resampled in register with the `target_reference`).
In other words, `fmap2data_xfm` is the result of calling a registration tool
be resampled in register with the ``target_reference``).
In other words, ``xfm_data2fmap`` is the result of calling a registration tool
such as ANTs configured for a linear transform with at most 12 degrees of freedom
and called with the image carrying `target_affine` as reference and the fieldmap
and called with the image carrying ``target_affine`` as reference and the fieldmap
reference as moving.
The result of such a registration framework is an affine (our `fmap2data_xfm` here)
The result of such a registration framework is an affine (our ``xfm_data2fmap`` here)
that maps coordinates in reference (target) RAS onto the fieldmap RAS.
approx : :obj:`bool`
If ``True``, do not reconstruct the B-Spline field directly on the target
Expand Down Expand Up @@ -456,9 +456,9 @@ def apply(
)
else:
# Generate warp field (before ensuring positive cosines)
self.fit(moving, fmap2data_xfm=fmap2data_xfm, approx=approx)
self.fit(moving, xfm_data2fmap=xfm_data2fmap, approx=approx)

# Squeeze non-spatial dimensions
# Squeeze n33on-spatial dimensions
newshape = moving.shape[:3] + tuple(dim for dim in moving.shape[3:] if dim > 1)
data = np.reshape(moving.dataobj, newshape)
ndim = min(data.ndim, 3)
Expand Down
2 changes: 1 addition & 1 deletion sdcflows/workflows/apply/tests/test_correction.py
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ def test_unwarp_wf(tmpdir, datadir, workdir, outdir, with_affine):
]

if with_affine:
workflow.inputs.inputnode.data2fmap_xfm = str(
workflow.inputs.inputnode.fmap2data_xfm = str(
str(derivs_path / "sub-101006_from-sbrefLR_to-fieldmapref_mode-image_xfm.mat")
)

Expand Down

0 comments on commit 53663d4

Please sign in to comment.