From 6f6095dc29413422c883657c2a4453265793880d Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 10 Jan 2024 09:57:00 -0500 Subject: [PATCH] FIX: Derive field from SyN displacements using EPI affine (#421) Theoretically, a fieldmap in Hz is the number of voxels/second that signal is shifted from the true location. SyN-SDC calculates the shift in mm from the true location. In order to find the equivalent fieldmap, we need the size of the voxels, readout time and phase-encoding direction of the EPI image that is to be corrected. Previously, we were deriving the fieldmap from the readout time and phase-encoding direction of the EPI image, but the voxel size of the displacement field. This causes an overestimate of the number of voxels to be shifted, and thus exaggerated correction. --- sdcflows/interfaces/fmap.py | 2 ++ sdcflows/tests/test_transform.py | 2 +- sdcflows/transform.py | 19 +++++++++++-------- sdcflows/workflows/fit/syn.py | 13 +++++++------ 4 files changed, 21 insertions(+), 15 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/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"), 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..a0783a6deb 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") @@ -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")]), @@ -314,6 +315,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")]), @@ -329,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