[ENH] Add phase1/phase2 SDC support #1359
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -21,7 +21,7 @@ | |
| from nipype.utils.filemanip import fname_presuffix | ||
| from nipype.interfaces.base import ( | ||
| BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits, | ||
| SimpleInterface, InputMultiPath) | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
| LOGGER = logging.getLogger('nipype.interface') | ||
|
|
||
|
|
@@ -32,10 +32,13 @@ class FieldEnhanceInputSpec(BaseInterfaceInputSpec): | |
| in_magnitude = File(exists=True, desc='input magnitude') | ||
| unwrap = traits.Bool(False, usedefault=True, desc='run phase unwrap') | ||
| despike = traits.Bool(True, usedefault=True, desc='run despike filter') | ||
| bspline_smooth = traits.Bool(True, usedefault=True, | ||
| desc='run 3D bspline smoother') | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| mask_erode = traits.Int(1, usedefault=True, desc='mask erosion iterations') | ||
| despike_threshold = traits.Float(0.2, usedefault=True, | ||
| desc='mask erosion iterations') | ||
| num_threads = traits.Int(1, usedefault=True, nohash=True, | ||
| desc='number of jobs') | ||
|
|
||
|
|
||
| class FieldEnhanceOutputSpec(TraitedSpec): | ||
|
|
@@ -68,7 +71,8 @@ def _run_interface(self, runtime): | |
|
|
||
| # Dilate mask | ||
| if self.inputs.mask_erode > 0: | ||
| struc = sim.iterate_structure( | ||
| sim.generate_binary_structure(3, 2), 1) | ||
| mask = sim.binary_erosion( | ||
| mask, struc, | ||
| iterations=self.inputs.mask_erode | ||
|
|
@@ -108,7 +112,8 @@ def _run_interface(self, runtime): | |
|
|
||
| nslices = 0 | ||
| try: | ||
| errorslice = np.squeeze( | ||
| np.argwhere(errormask.sum(0).sum(0) > 0)) | ||
| nslices = errorslice[-1] - errorslice[0] | ||
| except IndexError: # mask is empty, do not refine | ||
| pass | ||
|
|
@@ -119,7 +124,8 @@ def _run_interface(self, runtime): | |
| diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk, | ||
| datanii.affine, datanii.header) | ||
|
|
||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| bspobj2 = fbsp.BSplineFieldmap(diffmapnii, | ||
| knots_zooms=[24., 24., 4.], | ||
| njobs=self.inputs.num_threads) | ||
| bspobj2.fit() | ||
| smoothed2 = bspobj2.get_smoothed().get_data() | ||
|
|
@@ -207,6 +213,38 @@ def _run_interface(self, runtime): | |
| return runtime | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
|
|
||
| class Phases2FieldmapInputSpec(BaseInterfaceInputSpec): | ||
| phase_files = InputMultiPath(File(exists=True), mandatory=True, | ||
| desc='list of phase1, phase2 files') | ||
| magnitude_files = InputMultiPath(File(exists=True), mandatory=True, | ||
| desc='list of phase1, phase2 files') | ||
| metadatas = traits.List(traits.Dict, mandatory=True, | ||
| desc='list of phase1, phase2 metadata dicts') | ||
|
|
||
|
|
||
| class Phases2FieldmapOutputSpec(TraitedSpec): | ||
| out_file = File(desc='the output fieldmap') | ||
| phasediff_metadata = traits.Dict(desc='the phasediff metadata') | ||
|
|
||
|
|
||
| class Phases2Fieldmap(SimpleInterface): | ||
| """ | ||
| Convert a phase1, phase2 into a difference map | ||
| """ | ||
| input_spec = Phases2FieldmapInputSpec | ||
| output_spec = Phases2FieldmapOutputSpec | ||
|
|
||
| def _run_interface(self, runtime): | ||
| # Get the echo times | ||
| fmap_file, merged_metadata = phases2fmap( | ||
| self.inputs.phase_files, self.inputs.magnitude_files, | ||
| self.inputs.metadatas | ||
| ) | ||
| self._results['phasediff_metadata'] = merged_metadata | ||
| self._results['out_file'] = fmap_file | ||
| return runtime | ||
|
|
||
|
|
||
| def _despike2d(data, thres, neigh=None): | ||
| """ | ||
| despiking as done in FSL fugue | ||
|
|
@@ -252,14 +290,16 @@ def _unwrap(fmap_data, mag_file, mask=None): | |
|
|
||
| nb.Nifti1Image(fmap_data, magnii.affine).to_filename('fmap_rad.nii.gz') | ||
| nb.Nifti1Image(mask, magnii.affine).to_filename('fmap_mask.nii.gz') | ||
| nb.Nifti1Image(magnii.get_data(), magnii.affine | ||
| ).to_filename('fmap_mag.nii.gz') | ||
|
|
||
| # Run prelude | ||
| res = PRELUDE(phase_file='fmap_rad.nii.gz', | ||
| magnitude_file='fmap_mag.nii.gz', | ||
| mask_file='fmap_mask.nii.gz').run() | ||
|
|
||
| unwrapped = nb.load(res.outputs.unwrapped_phase_file | ||
| ).get_data() * (fmapmax / pi) | ||
| return unwrapped | ||
|
|
||
|
|
||
|
|
@@ -290,8 +330,9 @@ def get_ees(in_meta, in_file=None): | |
|
|
||
| = T_\\text{ro} \\, (N_\\text{PE} / f_\\text{acc} - 1)^{-1} | ||
|
|
||
| where :math:`N_y` is the number of pixels along the phase-encoding | ||
| direction :math:`y`, and :math:`f_\\text{acc}` is the parallel imaging | ||
| acceleration factor | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| (:abbr:`GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition)`, | ||
| :abbr:`ARC (Autocalibrating Reconstruction for Cartesian imaging)`, etc.). | ||
|
|
||
|
|
@@ -471,7 +512,8 @@ def phdiff2fmap(in_file, delta_te, newpath=None): | |
|
|
||
| .. math:: | ||
|
|
||
| \Delta B_0 (\text{T}^{-1}) = \frac{\Delta\Theta}{\ | ||
| 2\pi\gamma\Delta\text{TE}} | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
|
|
||
| In this case, we do not take into account the gyromagnetic ratio of the | ||
|
|
@@ -497,6 +539,56 @@ def phdiff2fmap(in_file, delta_te, newpath=None): | |
| return out_file | ||
|
|
||
|
|
||
| def phases2fmap(phase_files, magnitude_files, metadatas, newpath=None): | ||
| """calculates a phasediff from two phase images. Assumes monopolar | ||
| readout and that the second image in the list corresponds to the | ||
| longer echo time""" | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| import numpy as np | ||
| import nibabel as nb | ||
| from nipype.utils.filemanip import fname_presuffix | ||
| from copy import deepcopy | ||
|
|
||
| phasediff_file = fname_presuffix(phase_files[0], suffix='_phasediff', | ||
| newpath=newpath) | ||
| echo_times = [meta.get("EchoTime") for meta in metadatas] | ||
| if None in echo_times or echo_times[0] == echo_times[1]: | ||
| raise RuntimeError() | ||
| # Determine the order of subtraction | ||
| short_echo_index = echo_times.index(min(echo_times)) | ||
| long_echo_index = echo_times.index(max(echo_times)) | ||
|
|
||
| magnitude_image = magnitude_files[short_echo_index] | ||
| short_phase_image = phase_files[short_echo_index] | ||
| long_phase_image = phase_files[long_echo_index] | ||
|
|
||
| mag_img = nb.load(magnitude_image) | ||
| mag = mag_img.get_data().astype(np.float) | ||
| image0 = nb.load(short_phase_image) | ||
| phase0 = image0.get_data().astype(np.float) | ||
| image1 = nb.load(long_phase_image) | ||
| phase1 = image1.get_data().astype(np.float) | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
|
|
||
| # Calculate fieldmaps | ||
| rad0 = np.pi * phase0 / 2048 - np.pi | ||
| rad1 = np.pi * phase1 / 2048 - np.pi | ||
|
There was a problem hiding this comment. Will phases always range from 0 to 4096? There was a problem hiding this comment. We're trying to recreate these scanning sequences to see if this is still true for the current Siemens product. This has been true for all the older versions of the sequence I've run into. There was a problem hiding this comment. Thanks for the diligence on this. Do you know if any other scanners produce phase1/phase2 fieldmaps? It would be good to see if this convention extends beyond Siemens. Perhaps @neurolabusc might have the experience or know someone who does. There was a problem hiding this comment. I have a web page that is relevant. For Siemens images I think the scaled range is -4096..4092 (scl_slope = 2.000000; scl_inter = -4096.000000). For Philips I think the scaled range is ~-500..496 for some systems and -pi..+pi for others. For GE I think the convention is to save real/imaginary images so you would use fslcomplex to get the phase map. As a good estimate, you can assume min..max should be scaled very close to -pi..+pi. There was a problem hiding this comment. Thanks for the rundown. If the slope and intercept are 2, -4096, and the original range is 0..4096, then this should be: rad0 = np.pi * phase0 / 4096
rad1 = np.pi * phase1 / 4096This is because I think we can ignore the GE case for now; that seems like it should probably not be considered "phase1"/"phase2" fieldmaps, and it's a job for BIDS to establish what the correct naming scheme is to indicate those files. I think in the short term, we could have some a heuristic where we check for values above 500 and pi, to establish the correct denominator. In the long term, I'm leery of putting too many heuristics for interpreting data into fMRIPrep that are not in the BIDS spec. We may want to consider a metadata field that unambiguously indicates the correct interpretation. (Though the heuristic would probably need to stay, to keep processing datasets that conform with BIDS 1.0.) There was a problem hiding this comment. I see. So I guess we need to watch for cases where the scale and intercept are missing. There was a problem hiding this comment. We scanned all the available GRE sequences on a prisma and uploaded them on openneuro in case you want to take a look. Here's the link: https://openneuro.org/datasets/ds001600. These were converted with a recent I'm open to suggestions on how to handle these robustly. Also, sorry about the formatting mess. I ran a misconfigured yapf and am still putting things back. There was a problem hiding this comment. @mattcieslak can you send me the raw DICOM data and permission to distribute the DICOMs on the dcm2niix wiki? This would provide a nice reference set for all DICOM-to-BIDS conversion tools to test, and it would provide nice validation for any future changes or features. There was a problem hiding this comment. Despite being acquired on a E11C Prisma that supports high dynamic range (16-bit DICOM), all the samples you acquired appear to use 12-bit precision (0..4095). My guess is 12-bit is more than sufficient for the SNR provided, but I would make sure that any code generates an error if image intensities are outside this range. This will detect any sequences that provide more than 12-bit precision.
There was a problem hiding this comment. @neurolabusc sent a link to your mailbox.sc.edu from my pennmedicine email |
||
| a = mag * np.cos(rad0) | ||
| b = mag * np.sin(rad0) | ||
| c = mag * np.cos(rad1) | ||
| d = mag * np.sin(rad1) | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| fmap = -np.arctan2(b*c-a*d, a*c+b*d) | ||
|
|
||
| phasediff_nii = nb.Nifti1Image(fmap, image0.affine) | ||
| phasediff_nii.set_data_dtype(np.float32) | ||
| phasediff_nii.to_filename(phasediff_file) | ||
|
|
||
| merged_metadata = deepcopy(metadatas[0]) | ||
| del merged_metadata['EchoTime'] | ||
| merged_metadata['EchoTime1'] = float(echo_times[short_echo_index]) | ||
| merged_metadata['EchoTime2'] = float(echo_times[long_echo_index]) | ||
|
|
||
| return phasediff_file, merged_metadata | ||
|
|
||
|
|
||
| def _delta_te(in_values, te1=None, te2=None): | ||
| """Read :math:`\Delta_\text{TE}` from BIDS metadata dict""" | ||
| if isinstance(in_values, float): | ||
|
|
@@ -520,13 +612,15 @@ def _delta_te(in_values, te1=None, te2=None): | |
|
|
||
| # For convienience if both are missing we should give one error about them | ||
| if te1 is None and te2 is None: | ||
| raise RuntimeError('EchoTime1 and EchoTime2 metadata fields not found.' | ||
| ' Please consult the BIDS specification.') | ||
| if te1 is None: | ||
| raise RuntimeError( | ||
| 'EchoTime1 metadata field not found. Please consult the BIDS ' | ||
| 'specification.') | ||
| if te2 is None: | ||
| raise RuntimeError( | ||
| 'EchoTime2 metadata field not found. Please consult the BIDS ' | ||
| 'specification.') | ||
|
|
||
| return abs(float(te2) - float(te1)) | ||
mattcieslak marked this conversation as resolved.
Show resolved
Hide resolved
|
||
Uh oh!
There was an error while loading. Please reload this page.