From e56852600a45cfa5ea49224d080cc5736ac74980 Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 9 Jan 2024 11:46:25 -0500 Subject: [PATCH 1/3] FIX: Derive field from SyN displacements using EPI affine --- sdcflows/interfaces/fmap.py | 2 ++ sdcflows/transform.py | 19 +++++++++++-------- sdcflows/workflows/fit/syn.py | 3 ++- 3 files changed, 15 insertions(+), 9 deletions(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 58dc4e73f4..6327460faf 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -161,6 +161,7 @@ def _run_interface(self, runtime): class _DisplacementsField2FieldmapInputSpec(BaseInterfaceInputSpec): transform = File(exists=True, mandatory=True, desc="input displacements field") + epi = File(exists=True, mandatory=True, desc="source EPI image") ro_time = traits.Float(mandatory=True, desc="total readout time") pe_dir = traits.Enum( "j-", "j", "i", "i-", "k", "k-", mandatory=True, desc="phase encoding direction" @@ -189,6 +190,7 @@ def _run_interface(self, runtime): ) fmapnii = disp_to_fmap( nb.load(self.inputs.transform), + nb.load(self.inputs.epi), ro_time=self.inputs.ro_time, pe_dir=self.inputs.pe_dir, itk_format=self.inputs.itk_transform, diff --git a/sdcflows/transform.py b/sdcflows/transform.py index a53b7224a0..fff9627e10 100644 --- a/sdcflows/transform.py +++ b/sdcflows/transform.py @@ -617,7 +617,7 @@ def fmap_to_disp(fmap_nii, ro_time, pe_dir, itk_format=True): return xyz_nii -def disp_to_fmap(xyz_nii, ro_time, pe_dir, itk_format=True): +def disp_to_fmap(xyz_nii, epi_nii, ro_time, pe_dir, itk_format=True): """ Convert a displacements field into a fieldmap in Hz. @@ -625,12 +625,14 @@ def disp_to_fmap(xyz_nii, ro_time, pe_dir, itk_format=True): Parameters ---------- - xyz_nii : :obj:`os.pathlike` - Path to a displacements field in NIfTI format. + xyz_nii : :class:`nibabel.Nifti1Image` + Displacements field in NIfTI format. + epi_nii : :class:`nibabel.Nifti1Image` + Source EPI image. ro_time : :obj:`float` - The total readout time in seconds. + The total readout time of the EPI image in seconds. pe_dir : :obj:`str` - The ``PhaseEncodingDirection`` metadata value. + The ``PhaseEncodingDirection`` metadata value of the EPI image. Returns ------- @@ -644,12 +646,13 @@ def disp_to_fmap(xyz_nii, ro_time, pe_dir, itk_format=True): # ITK displacement vectors are in LPS orientation xyz_deltas[:, (0, 1)] *= -1 - inv_aff = np.linalg.inv(xyz_nii.affine) - inv_aff[:3, 3] = 0 # Translations MUST NOT be applied. + # Use the EPI's affine to determine voxel spacing, axis ordering and flips + inv_aff = np.linalg.inv(epi_nii.affine) + inv_mat = inv_aff[:3, :3] # Convert displacements from mm to voxel units # Using the inverse affine accounts for reordering of axes, etc. - ijk_deltas = nb.affines.apply_affine(inv_aff, xyz_deltas).astype("float32") + ijk_deltas = (xyz_deltas @ inv_mat.T).astype("float32") pe_axis = "ijk".index(pe_dir[0]) vsm = ijk_deltas[:, pe_axis].reshape(xyz_nii.shape[:3]) scale_factor = -ro_time if pe_dir.endswith("-") else ro_time diff --git a/sdcflows/workflows/fit/syn.py b/sdcflows/workflows/fit/syn.py index 5827c5487d..5275990cac 100644 --- a/sdcflows/workflows/fit/syn.py +++ b/sdcflows/workflows/fit/syn.py @@ -247,7 +247,7 @@ def init_syn_sdc_wf( # Extract the corresponding fieldmap in Hz extract_field = pe.Node( - DisplacementsField2Fieldmap(demean=True), name="extract_field" + DisplacementsField2Fieldmap(), name="extract_field" ) unwarp = pe.Node(ApplyCoeffsField(), name="unwarp") @@ -314,6 +314,7 @@ def init_syn_sdc_wf( (epi_merge, syn, [("out", "moving_image")]), (moving_masks, syn, [("out", "moving_image_masks")]), (syn, extract_field, [(("forward_transforms", _pop), "transform")]), + (clip_epi, extract_field, [("out_file", "epi")]), (readout_time, extract_field, [("readout_time", "ro_time"), ("pe_direction", "pe_dir")]), (extract_field, zooms_field, [("out_file", "input_image")]), From 56454519b407e6b4b9cb95a66d41f6923e4bae2d Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 9 Jan 2024 11:50:08 -0500 Subject: [PATCH 2/3] FIX: Do not recenter exact displacements --- sdcflows/workflows/fit/syn.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sdcflows/workflows/fit/syn.py b/sdcflows/workflows/fit/syn.py index 5275990cac..a0783a6deb 100644 --- a/sdcflows/workflows/fit/syn.py +++ b/sdcflows/workflows/fit/syn.py @@ -267,14 +267,15 @@ def init_syn_sdc_wf( ) # Regularize with B-Splines - bs_filter = pe.Node(BSplineApprox(), name="bs_filter") + bs_filter = pe.Node( + BSplineApprox(recenter=False, debug=debug, extrapolate=not debug), + name="bs_filter", + ) bs_filter.interface._always_run = debug bs_filter.inputs.bs_spacing = ( [DEFAULT_LF_ZOOMS_MM, DEFAULT_HF_ZOOMS_MM] if not sloppy else [DEFAULT_ZOOMS_MM] ) - bs_filter.inputs.extrapolate = not debug - # fmt: off workflow.connect([ (inputnode, readout_time, [(("epi_ref", _pop), "in_file"), (("epi_ref", _pull), "metadata")]), @@ -330,8 +331,7 @@ def init_syn_sdc_wf( (bs_filter, outputnode, [("out_coeff", "fmap_coeff")]), (unwarp, outputnode, [("out_corrected", "fmap_ref"), ("out_field", "fmap")]), - ]) - # fmt: on + ]) # fmt:skip return workflow From 7880eea8d7c9763254656972922aec3bb69a36ca Mon Sep 17 00:00:00 2001 From: "Christopher J. Markiewicz" Date: Tue, 9 Jan 2024 16:53:56 -0500 Subject: [PATCH 3/3] TEST: Update test to pass explicit EPI image --- sdcflows/tests/test_transform.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sdcflows/tests/test_transform.py b/sdcflows/tests/test_transform.py index eb8bd09bb2..72e44e0759 100644 --- a/sdcflows/tests/test_transform.py +++ b/sdcflows/tests/test_transform.py @@ -327,11 +327,11 @@ def test_conversions(tmpdir, testdata_dir, pe_dir): ro_time=0.2, pe_dir=pe_dir, ), + fmap_nii, ro_time=0.2, pe_dir=pe_dir, ) - new_nii.to_filename("test.nii.gz") assert np.allclose( fmap_nii.get_fdata(dtype="float32"), new_nii.get_fdata(dtype="float32"),