diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 28ca085f5c..0a8da2e84b 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -24,6 +24,43 @@ LOGGER = logging.getLogger('nipype.interface') +class _SubtractPhasesInputSpec(BaseInterfaceInputSpec): + in_phases = traits.List(File(exists=True), min=1, max=2, + desc='input phase maps') + in_meta = traits.List(traits.Dict(), min=1, max=2, + desc='metadata corresponding to the inputs') + + +class _SubtractPhasesOutputSpec(TraitedSpec): + phase_diff = File(exists=True, desc='phase difference map') + metadata = traits.Dict(desc='output metadata') + + +class SubtractPhases(SimpleInterface): + """Calculate a phase difference map.""" + + input_spec = _SubtractPhasesInputSpec + output_spec = _SubtractPhasesOutputSpec + + def _run_interface(self, runtime): + if len(self.inputs.in_phases) != len(self.inputs.in_meta): + raise ValueError( + 'Length of input phase-difference maps and metadata files ' + 'should match.') + + if len(self.inputs.in_phases) == 1: + self._results['phase_diff'] = self.inputs.in_phases[0] + self._results['metadata'] = self.inputs.in_meta[0] + return runtime + + self._results['phase_diff'], self._results['metadata'] = \ + _subtract_phases(self.inputs.in_phases, + self.inputs.in_meta, + newpath=runtime.cwd) + + return runtime + + class _FieldEnhanceInputSpec(BaseInterfaceInputSpec): in_file = File(exists=True, mandatory=True, desc='input fieldmap') in_mask = File(exists=True, desc='brain mask') @@ -649,3 +686,36 @@ def au2rads(in_file, newpath=None): out_file = fname_presuffix(in_file, suffix='_rads', newpath=newpath) nb.Nifti1Image(data, im.affine, hdr).to_filename(out_file) return out_file + + +def _subtract_phases(in_phases, in_meta, newpath=None): + # Discard traits with copy(), so that pop() works. + in_meta = (in_meta[0].copy(), in_meta[1].copy()) + echo_times = tuple([m.pop('EchoTime', None) for m in in_meta]) + if not all(echo_times): + raise ValueError( + 'One or more missing EchoTime metadata parameter ' + 'associated to one or more phase map(s).') + + if echo_times[0] > echo_times[1]: + in_phases = (in_phases[1], in_phases[0]) + in_meta = (in_meta[1], in_meta[0]) + echo_times = (echo_times[1], echo_times[0]) + + in_phases_nii = [nb.load(ph) for ph in in_phases] + sub_data = in_phases_nii[1].get_fdata(dtype='float32') - \ + in_phases_nii[0].get_fdata(dtype='float32') + + new_meta = in_meta[1].copy() + new_meta.update(in_meta[0]) + new_meta['EchoTime1'] = echo_times[0] + new_meta['EchoTime2'] = echo_times[1] + + hdr = in_phases_nii[0].header.copy() + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units('mm') + nii = nb.Nifti1Image(sub_data, in_phases_nii[0].affine, hdr) + out_phdiff = fname_presuffix(in_phases[0], suffix='_phdiff', + newpath=newpath) + nii.to_filename(out_phdiff) + return out_phdiff, new_meta diff --git a/sdcflows/workflows/phdiff.py b/sdcflows/workflows/phdiff.py index df97c342b6..c4c576f120 100644 --- a/sdcflows/workflows/phdiff.py +++ b/sdcflows/workflows/phdiff.py @@ -21,7 +21,7 @@ from nipype.pipeline import engine as pe from niworkflows.engine.workflows import LiterateWorkflow as Workflow -from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads +from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads, SubtractPhases from .gre import init_fmap_postproc_wf, init_magnitude_wf @@ -107,6 +107,9 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): # FSL PRELUDE will perform phase-unwrapping prelude = pe.MapNode(fsl.PRELUDE(), iterfield=['phase_file'], name='prelude') + calc_phdiff = pe.Node(SubtractPhases(), name='calc_phdiff', + run_without_submitting=True) + fmap_postproc_wf = init_fmap_postproc_wf(omp_nthreads=omp_nthreads, fmap_bspline=False) compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap') @@ -118,8 +121,11 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): ('outputnode.fmap_mask', 'mask_file')]), (split, phmap2rads, [('map_file', 'in_file')]), (phmap2rads, prelude, [('out_file', 'phase_file')]), - (split, fmap_postproc_wf, [('meta', 'inputnode.metadata')]), - (prelude, fmap_postproc_wf, [('unwrapped_phase_file', 'inputnode.fmap')]), + (prelude, calc_phdiff, [('unwrapped_phase_file', 'in_phases')]), + (split, calc_phdiff, [('meta', 'in_meta')]), + (calc_phdiff, fmap_postproc_wf, [ + ('phase_diff', 'inputnode.fmap'), + ('metadata', 'inputnode.metadata')]), (magnitude_wf, fmap_postproc_wf, [ ('outputnode.fmap_mask', 'inputnode.fmap_mask'), ('outputnode.fmap_ref', 'inputnode.fmap_ref')]), diff --git a/sdcflows/workflows/tests/test_phdiff.py b/sdcflows/workflows/tests/test_phdiff.py index 191259e3c0..03fa42bbe1 100644 --- a/sdcflows/workflows/tests/test_phdiff.py +++ b/sdcflows/workflows/tests/test_phdiff.py @@ -10,20 +10,77 @@ 'ds001600', 'testdata', ]) -def test_workflow(bids_layouts, tmpdir, output_path, dataset, workdir): +def test_phdiff(bids_layouts, tmpdir, output_path, dataset, workdir): """Test creation of the workflow.""" tmpdir.chdir() + extra_entities = {} + if dataset == 'ds001600': + extra_entities['acquisition'] = 'v4' + data = bids_layouts[dataset] wf = Workflow(name='phdiff_%s' % dataset) phdiff_wf = init_phdiff_wf(omp_nthreads=2) phdiff_wf.inputs.inputnode.magnitude = data.get( suffix=['magnitude1', 'magnitude2'], - acq='v4', + return_type='file', + extension=['.nii', '.nii.gz'], + **extra_entities) + + phdiff_files = data.get( + suffix='phasediff', + extension=['.nii', '.nii.gz'], + **extra_entities) + + phdiff_wf.inputs.inputnode.phasediff = [ + (ph.path, ph.get_metadata()) for ph in phdiff_files] + + if output_path: + from ...interfaces.reportlets import FieldmapReportlet + rep = pe.Node(FieldmapReportlet(reference_label='Magnitude'), 'simple_report') + rep.interface._always_run = True + + dsink = pe.Node(DerivativesDataSink( + base_directory=str(output_path), keep_dtype=True), name='dsink') + dsink.interface.out_path_base = 'sdcflows' + dsink.inputs.source_file = phdiff_files[0].path + + dsink_fmap = pe.Node(DerivativesDataSink( + base_directory=str(output_path), keep_dtype=True), name='dsink_fmap') + dsink_fmap.interface.out_path_base = 'sdcflows' + dsink_fmap.inputs.source_file = phdiff_files[0].path + + wf.connect([ + (phdiff_wf, rep, [ + ('outputnode.fmap', 'fieldmap'), + ('outputnode.fmap_ref', 'reference'), + ('outputnode.fmap_mask', 'mask')]), + (rep, dsink, [('out_report', 'in_file')]), + (phdiff_wf, dsink_fmap, [('outputnode.fmap', 'in_file')]), + ]) + else: + wf.add_nodes([phdiff_wf]) + + if workdir: + wf.base_dir = str(workdir) + + wf.run() + + +def test_phases(bids_layouts, tmpdir, output_path, workdir): + """Test creation of the workflow.""" + tmpdir.chdir() + + data = bids_layouts['ds001600'] + wf = Workflow(name='phases_ds001600') + phdiff_wf = init_phdiff_wf(omp_nthreads=2) + phdiff_wf.inputs.inputnode.magnitude = data.get( + suffix=['magnitude1', 'magnitude2'], + acquisition='v2', return_type='file', extension=['.nii', '.nii.gz']) - phdiff_files = data.get(suffix='phasediff', acq='v4', + phdiff_files = data.get(suffix=['phase1', 'phase2'], acquisition='v2', extension=['.nii', '.nii.gz']) phdiff_wf.inputs.inputnode.phasediff = [ @@ -39,12 +96,18 @@ def test_workflow(bids_layouts, tmpdir, output_path, dataset, workdir): dsink.interface.out_path_base = 'sdcflows' dsink.inputs.source_file = phdiff_files[0].path + dsink_fmap = pe.Node(DerivativesDataSink( + base_directory=str(output_path), keep_dtype=True), name='dsink_fmap') + dsink_fmap.interface.out_path_base = 'sdcflows' + dsink_fmap.inputs.source_file = phdiff_files[0].path + wf.connect([ (phdiff_wf, rep, [ ('outputnode.fmap', 'fieldmap'), ('outputnode.fmap_ref', 'reference'), ('outputnode.fmap_mask', 'mask')]), (rep, dsink, [('out_report', 'in_file')]), + (phdiff_wf, dsink_fmap, [('outputnode.fmap', 'in_file')]), ]) else: wf.add_nodes([phdiff_wf])