ENH: Add image.Rescale interface

Author: effigies

nipype/interfaces/image.py

@@ -9,6 +9,84 @@
 from .base import (SimpleInterface, TraitedSpec, BaseInterfaceInputSpec,
                    traits, File)
 
+
+class RescaleInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True,
+                   desc='Skull-stripped image to rescale')
+    ref_file = File(exists=True, mandatory=True,
+                    desc='Skull-stripped reference image')
+    invert = traits.Bool(desc='Invert contrast of rescaled image')
+    percentile = traits.Range(low=0., high=50., value=0., usedefault=True,
+                              desc='Percentile to use for reference to allow '
+                                   'for outliers - 1 indicates the 1st and '
+                                   '99th percentiles in the input file will '
+                                   'be mapped to the 99th and 1st percentiles '
+                                   'in the reference; 0 indicates minima and '
+                                   'maxima will be mapped')
+
+
+class RescaleOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='Rescaled image')
+
+
+class Rescale(SimpleInterface):
+    """Rescale an image
+
+    Rescales the non-zero portion of ``in_file`` to match the bounds of the '
+    non-zero portion of ``ref_file``.
+    Reference values in the input and reference images are defined by the
+    ``percentile`` parameter, and the reference values in each image are
+    identified and the remaining values are scaled accordingly.
+    In the case of ``percentile == 0``, the reference values are the maxima
+    and minima of each image.
+    If the ``invert`` parameter is set, the input file is inverted prior to
+    rescaling.
+
+    Examples
+    --------
+
+    To use a high-resolution T1w image as a registration target for a T2\*
+    image, it may be useful to invert the T1w image and rescale to the T2\*
+    range.
+    Using the 1st and 99th percentiles may reduce the impact of outlier
+    voxels.
+
+    >>> from nipype.interfaces.image import Rescale
+    >>> invert_t1w = Reorient(invert=True)
+    >>> invert_t1w.inputs.in_file = 'structural.nii'
+    >>> invert_t1w.inputs.ref_file = 'functional.nii'
+    >>> invert_t1w.inputs.percentile = 1.
+    >>> res = invert_t1w.run()  # doctest: +SKIP
+
+    """
+    input_spec = RescaleInputSpec
+    output_spec = RescaleOutputSpec
+
+    def _run_interface(self, runtime):
+        img = nb.load(self.inputs.in_file)
+        data = img.get_data()
+        ref_data = nb.load(self.inputs.ref_file).get_data()
+
+        in_mask = data > 0
+        ref_mask = ref_data > 0
+
+        q = [self.inputs.percentile, 100. - self.inputs.percentile]
+        in_low, in_high = np.percentile(data[in_mask], q)
+        ref_low, ref_high = np.percentile(ref_data[ref_mask], q)
+        scale_factor = (ref_high - ref_low) / (in_high - in_low)
+
+        signal = in_high - data if self.inputs.invert else data - in_low
+        out_data = in_mask * (signal * scale_factor + ref_low)
+
+        suffix = '_inv' if self.inputs.invert else '_rescaled'
+        out_file = fname_presuffix(self.inputs.in_file, suffix=suffix,
+                                   newpath=runtime.cwd)
+        img.__class__(out_data, img.affine, img.header).to_filename(out_file)
+
+        self._results['out_file'] = out_file
+        return runtime
+
+
 _axes = ('RL', 'AP', 'SI')
 _orientations = tuple(
     ''.join((x[i], y[j], z[k]))