|
| 1 | +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- |
| 2 | +# vi: set ft=python sts=4 ts=4 sw=4 et: |
| 3 | +""" |
| 4 | +Interfaces to deal with the various types of fieldmap sources. |
| 5 | +
|
| 6 | + .. testsetup:: |
| 7 | +
|
| 8 | + >>> tmpdir = getfixture('tmpdir') |
| 9 | + >>> tmp = tmpdir.chdir() # changing to a temporary directory |
| 10 | + >>> nb.Nifti1Image(np.zeros((90, 90, 60)), None, None).to_filename( |
| 11 | + ... tmpdir.join('epi.nii.gz').strpath) |
| 12 | +
|
| 13 | +
|
| 14 | +""" |
| 15 | + |
| 16 | +from nipype.interfaces.base import ( |
| 17 | + BaseInterfaceInputSpec, |
| 18 | + TraitedSpec, |
| 19 | + File, |
| 20 | + isdefined, |
| 21 | + traits, |
| 22 | + SimpleInterface, |
| 23 | +) |
| 24 | + |
| 25 | + |
| 26 | +class _GetReadoutTimeInputSpec(BaseInterfaceInputSpec): |
| 27 | + in_file = File(exists=True, desc="EPI image corresponding to the metadata") |
| 28 | + metadata = traits.Dict(mandatory=True, desc="metadata corresponding to the inputs") |
| 29 | + |
| 30 | + |
| 31 | +class _GetReadoutTimeOutputSpec(TraitedSpec): |
| 32 | + readout_time = traits.Float |
| 33 | + |
| 34 | + |
| 35 | +class GetReadoutTime(SimpleInterface): |
| 36 | + """Calculate the readout time from available metadata.""" |
| 37 | + |
| 38 | + input_spec = _GetReadoutTimeInputSpec |
| 39 | + output_spec = _GetReadoutTimeOutputSpec |
| 40 | + |
| 41 | + def _run_interface(self, runtime): |
| 42 | + self._results["readout_time"] = get_trt( |
| 43 | + self.inputs.metadata, |
| 44 | + self.inputs.in_file if isdefined(self.inputs.in_file) else None, |
| 45 | + ) |
| 46 | + return runtime |
| 47 | + |
| 48 | + |
| 49 | +def get_trt(in_meta, in_file=None): |
| 50 | + r""" |
| 51 | + Extract the *total readout time* :math:`t_\text{RO}` from BIDS. |
| 52 | +
|
| 53 | + Calculate the *total readout time* for an input |
| 54 | + :abbr:`EPI (echo-planar imaging)` scan. |
| 55 | +
|
| 56 | + There are several procedures to calculate the total |
| 57 | + readout time. The basic one is that a ``TotalReadoutTime`` |
| 58 | + field is set in the JSON sidecar. The following examples |
| 59 | + use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the |
| 60 | + j-axis encoding direction. |
| 61 | +
|
| 62 | + >>> meta = {'TotalReadoutTime': 0.02596} |
| 63 | + >>> get_trt(meta) |
| 64 | + 0.02596 |
| 65 | +
|
| 66 | + If the *effective echo spacing* :math:`t_\text{ees}` |
| 67 | + (``EffectiveEchoSpacing`` BIDS field) is provided, then the |
| 68 | + total readout time can be calculated reading the number |
| 69 | + of voxels along the readout direction :math:`T_\text{ro}` |
| 70 | + and the parallel acceleration factor of the EPI :math:`f_\text{acc}`. |
| 71 | +
|
| 72 | + .. math :: |
| 73 | +
|
| 74 | + T_\text{ro} = t_\text{ees} \, (N_\text{PE} / f_\text{acc} - 1) |
| 75 | +
|
| 76 | + >>> meta = {'EffectiveEchoSpacing': 0.00059, |
| 77 | + ... 'PhaseEncodingDirection': 'j-', |
| 78 | + ... 'ParallelReductionFactorInPlane': 2} |
| 79 | + >>> get_trt(meta, in_file='epi.nii.gz') |
| 80 | + 0.02596 |
| 81 | +
|
| 82 | + Some vendors, like Philips, store different parameter names: |
| 83 | +
|
| 84 | + >>> meta = {'WaterFatShift': 8.129, |
| 85 | + ... 'MagneticFieldStrength': 3, |
| 86 | + ... 'PhaseEncodingDirection': 'j-', |
| 87 | + ... 'ParallelReductionFactorInPlane': 2} |
| 88 | + >>> get_trt(meta, in_file='epi.nii.gz') |
| 89 | + 0.018721183563864822 |
| 90 | +
|
| 91 | + """ |
| 92 | + import nibabel as nb |
| 93 | + |
| 94 | + # Use case 1: TRT is defined |
| 95 | + trt = in_meta.get("TotalReadoutTime", None) |
| 96 | + if trt is not None: |
| 97 | + return trt |
| 98 | + |
| 99 | + # All other cases require the parallel acc and npe (N vox in PE dir) |
| 100 | + acc = float(in_meta.get("ParallelReductionFactorInPlane", 1.0)) |
| 101 | + npe = nb.load(in_file).shape[_get_pe_index(in_meta)] |
| 102 | + etl = npe // acc |
| 103 | + |
| 104 | + # Use case 2: TRT is defined |
| 105 | + ees = in_meta.get("EffectiveEchoSpacing", None) |
| 106 | + if ees is not None: |
| 107 | + return ees * (etl - 1) |
| 108 | + |
| 109 | + # Use case 3 (philips scans) |
| 110 | + wfs = in_meta.get("WaterFatShift", None) |
| 111 | + if wfs is not None: |
| 112 | + fstrength = in_meta["MagneticFieldStrength"] |
| 113 | + wfd_ppm = 3.4 # water-fat diff in ppm |
| 114 | + g_ratio_mhz_t = 42.57 # gyromagnetic ratio for proton (1H) in MHz/T |
| 115 | + wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t |
| 116 | + return wfs / wfs_hz |
| 117 | + |
| 118 | + raise ValueError("Unknown total-readout time specification") |
| 119 | + |
| 120 | + |
| 121 | +def get_ees(in_meta, in_file=None): |
| 122 | + r""" |
| 123 | + Extract the *effective echo spacing* :math:`t_\text{ees}` from BIDS. |
| 124 | +
|
| 125 | + Calculate the *effective echo spacing* :math:`t_\text{ees}` |
| 126 | + for an input :abbr:`EPI (echo-planar imaging)` scan. |
| 127 | +
|
| 128 | +
|
| 129 | + There are several procedures to calculate the effective |
| 130 | + echo spacing. The basic one is that an ``EffectiveEchoSpacing`` |
| 131 | + field is set in the JSON sidecar. The following examples |
| 132 | + use an ``'epi.nii.gz'`` file-stub which has 90 pixels in the |
| 133 | + j-axis encoding direction. |
| 134 | +
|
| 135 | + >>> meta = {'EffectiveEchoSpacing': 0.00059, |
| 136 | + ... 'PhaseEncodingDirection': 'j-'} |
| 137 | + >>> get_ees(meta) |
| 138 | + 0.00059 |
| 139 | +
|
| 140 | + If the *total readout time* :math:`T_\text{ro}` (``TotalReadoutTime`` |
| 141 | + BIDS field) is provided, then the effective echo spacing can be |
| 142 | + calculated reading the number of voxels :math:`N_\text{PE}` along the |
| 143 | + readout direction and the parallel acceleration |
| 144 | + factor of the EPI |
| 145 | +
|
| 146 | + .. math :: |
| 147 | +
|
| 148 | + = T_\text{ro} \, (N_\text{PE} / f_\text{acc} - 1)^{-1} |
| 149 | +
|
| 150 | + where :math:`N_y` is the number of pixels along the phase-encoding direction |
| 151 | + :math:`y`, and :math:`f_\text{acc}` is the parallel imaging acceleration factor |
| 152 | + (:abbr:`GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition)`, |
| 153 | + :abbr:`ARC (Autocalibrating Reconstruction for Cartesian imaging)`, etc.). |
| 154 | +
|
| 155 | + >>> meta = {'TotalReadoutTime': 0.02596, |
| 156 | + ... 'PhaseEncodingDirection': 'j-', |
| 157 | + ... 'ParallelReductionFactorInPlane': 2} |
| 158 | + >>> get_ees(meta, in_file='epi.nii.gz') |
| 159 | + 0.00059 |
| 160 | +
|
| 161 | + Some vendors, like Philips, store different parameter names (see |
| 162 | + http://dbic.dartmouth.edu/pipermail/mrusers/attachments/20141112/eb1d20e6/attachment.pdf |
| 163 | + ): |
| 164 | +
|
| 165 | + >>> meta = {'WaterFatShift': 8.129, |
| 166 | + ... 'MagneticFieldStrength': 3, |
| 167 | + ... 'PhaseEncodingDirection': 'j-', |
| 168 | + ... 'ParallelReductionFactorInPlane': 2} |
| 169 | + >>> get_ees(meta, in_file='epi.nii.gz') |
| 170 | + 0.00041602630141921826 |
| 171 | +
|
| 172 | + """ |
| 173 | + import nibabel as nb |
| 174 | + from sdcflows.interfaces.epi import _get_pe_index |
| 175 | + |
| 176 | + # Use case 1: EES is defined |
| 177 | + ees = in_meta.get("EffectiveEchoSpacing", None) |
| 178 | + if ees is not None: |
| 179 | + return ees |
| 180 | + |
| 181 | + # All other cases require the parallel acc and npe (N vox in PE dir) |
| 182 | + acc = float(in_meta.get("ParallelReductionFactorInPlane", 1.0)) |
| 183 | + npe = nb.load(in_file).shape[_get_pe_index(in_meta)] |
| 184 | + etl = npe // acc |
| 185 | + |
| 186 | + # Use case 2: TRT is defined |
| 187 | + trt = in_meta.get("TotalReadoutTime", None) |
| 188 | + if trt is not None: |
| 189 | + return trt / (etl - 1) |
| 190 | + |
| 191 | + # Use case 3 (philips scans) |
| 192 | + wfs = in_meta.get("WaterFatShift", None) |
| 193 | + if wfs is not None: |
| 194 | + fstrength = in_meta["MagneticFieldStrength"] |
| 195 | + wfd_ppm = 3.4 # water-fat diff in ppm |
| 196 | + g_ratio_mhz_t = 42.57 # gyromagnetic ratio for proton (1H) in MHz/T |
| 197 | + wfs_hz = fstrength * wfd_ppm * g_ratio_mhz_t |
| 198 | + return wfs / (wfs_hz * etl) |
| 199 | + |
| 200 | + raise ValueError("Unknown effective echo-spacing specification") |
| 201 | + |
| 202 | + |
| 203 | +def _get_pe_index(meta): |
| 204 | + pe = meta["PhaseEncodingDirection"] |
| 205 | + try: |
| 206 | + return {"i": 0, "j": 1, "k": 2}[pe[0]] |
| 207 | + except KeyError: |
| 208 | + raise RuntimeError('"%s" is an invalid PE string' % pe) |
0 commit comments