Skip to content

Commit

Permalink
fix(interface): EstimateReferenceImage generates 4D references if…
Browse files Browse the repository at this point in the history
… 4D SBRefs are supplied as input.

This patch checks the dimensionality of SBRefs and generates an average if the SBRef is 4D.
Fixes nipreps/fmriprep#1579.
Intimately related to nipreps#253, although I'm not sure the proposed resampling there is absolutely necessary.
It could be detrimental for ``bbregister`` in co-registration.
It is very likely necessary for head-motion correction with either ``mcflirt`` or ``3dVolReg``.
  • Loading branch information
oesteban committed Apr 11, 2019
1 parent 69a228a commit 98d40bd
Showing 1 changed file with 28 additions and 18 deletions.
46 changes: 28 additions & 18 deletions niworkflows/interfaces/registration.py
Original file line number Diff line number Diff line change
Expand Up @@ -370,45 +370,55 @@ class EstimateReferenceImage(SimpleInterface):
output_spec = EstimateReferenceImageOutputSpec

def _run_interface(self, runtime):
in_nii = nb.load(self.inputs.in_file)
n_volumes_to_discard = _get_vols_to_discard(in_nii)
ref_name = self.inputs.in_file
ref_nii = nb.load(ref_name)
n_volumes_to_discard = _get_vols_to_discard(ref_nii)

self._results["n_volumes_to_discard"] = n_volumes_to_discard

out_ref_fname = os.path.join(runtime.cwd, "ref_bold.nii.gz")
if isdefined(self.inputs.sbref_file):
self._results['ref_image'] = self.inputs.sbref_file
return runtime
out_ref_fname = os.path.join(runtime.cwd, "ref_sbref.nii.gz")
ref_name = self.inputs.sbref_file
ref_nii = nb.squeeze_image(nb.load(ref_name))

# If reference is only 1 volume, return it directly
if ref_nii.get_data().ndim == 3:
ref_nii.header.extensions.clear()
ref_nii.to_filename(out_ref_fname)
self._results['ref_image'] = out_ref_fname
return runtime
else:
# Reset this variable as it no longer applies
# and value for the output is stored in self._results
n_volumes_to_discard = 0

# Slicing may induce inconsistencies with shape-dependent values in extensions.
# For now, remove all. If this turns out to be a mistake, we can select extensions
# that don't break pipeline stages.
in_nii.header.extensions.clear()

out_ref_fname = os.path.join(runtime.cwd, "ref_image.nii.gz")
ref_nii.header.extensions.clear()

if n_volumes_to_discard == 0:
if in_nii.shape[-1] > 40:
slice_fname = os.path.join(runtime.cwd, "slice.nii.gz")
nb.Nifti1Image(in_nii.dataobj[:, :, :, 20:40], in_nii.affine,
in_nii.header).to_filename(slice_fname)
else:
slice_fname = self.inputs.in_file
if ref_nii.shape[-1] > 40:
ref_name = os.path.join(runtime.cwd, "slice.nii.gz")
nb.Nifti1Image(ref_nii.dataobj[:, :, :, 20:40], ref_nii.affine,
ref_nii.header).to_filename(ref_name)

if self.inputs.mc_method == "AFNI":
res = afni.Volreg(in_file=slice_fname, args='-Fourier -twopass',
res = afni.Volreg(in_file=ref_name, args='-Fourier -twopass',
zpad=4, outputtype='NIFTI_GZ').run()
elif self.inputs.mc_method == "FSL":
res = fsl.MCFLIRT(in_file=slice_fname,
res = fsl.MCFLIRT(in_file=ref_name,
ref_vol=0, interpolation='sinc').run()
mc_slice_nii = nb.load(res.outputs.out_file)

median_image_data = np.median(mc_slice_nii.get_data(), axis=3)
else:
median_image_data = np.median(
in_nii.dataobj[:, :, :, :n_volumes_to_discard], axis=3)
ref_nii.dataobj[:, :, :, :n_volumes_to_discard], axis=3)

nb.Nifti1Image(median_image_data, in_nii.affine,
in_nii.header).to_filename(out_ref_fname)
nb.Nifti1Image(median_image_data, ref_nii.affine,
ref_nii.header).to_filename(out_ref_fname)

self._results["ref_image"] = out_ref_fname

Expand Down

0 comments on commit 98d40bd

Please # to comment.