Skip to content

ENH: Add Rescale interface #2599

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jun 6, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
78 changes: 78 additions & 0 deletions nipype/interfaces/image.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is very unlikely the same percentile works well for both ends of the range. I would split this option or make it a tuple.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm. I expected that people would probably only ever want 0 or 1 (but left room for tweaking), but I can add more flexibility than this. This leads me to wonder if we want to permit different percentiles for the input and reference images. In which case we need 4 values, not 2.

WDYT? If possible, it'd be nice to allow a single knob, to keep it simple until somebody needs more refined control.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

WDYT?:

    percentile = traits.Either(
        traits.Range(low=0., high=50, value=0),  # One-knob behavior
        traits.Tuple(
            traits.Either(None, traits.Range(low=0., high=50, value=0)),
            traits.Either(None, traits.Range(low=50., high=100, value=100)),
        ), desc='...')

Then, if it is not a tuple, it has your one-knob mirror behavior. Otherwise you have a tuple, and you can always set one of the boundaries.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My fixation for this is that T1 images typically show background pixels under the 15% and, except when there are WM hyperintensities, you probably want 98% or higher for the upper limit.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, if I'm reading this correct, we don't need to be concerned with supporting different percentiles for the input and the reference?

My fixation for this is that T1 images typically show background pixels under the 15% and, except when there are WM hyperintensities, you probably want 98% or higher for the upper limit.

This is true for skull-stripped images?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nope, but I think we should aim for this interface to work well in any situation.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, if I'm reading this correct, we don't need to be concerned with supporting different percentiles for the input and the reference?

Okay, gotcha. I'm always thinking of the input. I guess we should specify that this interface expects the reference file to be "consistent" (as in already thresholded).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should aim for this interface to work well in any situation.

I don't think this interface will work at all well in non-skull-stripped cases. Calculating the magic numbers needed won't be significantly easier or more robust than even mediocre skull stripping, if I were to make a guess.

I guess we should specify that this interface expects the reference file to be "consistent" (as in already thresholded).

I do say:

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')

And:

Rescales the non-zero portion of ``in_file`` to match the bounds of the
non-zero portion of ``ref_file``.

But I can very explicitly say in the docstring that they should be skull-stripped. And would you suggest not applying the supplied percentiles to the reference image, as well?

I guess what's your vision for what this interface should do when percentiles != (0, 100)?

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 = Rescale(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]))
Expand Down
29 changes: 29 additions & 0 deletions nipype/interfaces/tests/test_auto_Rescale.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
from __future__ import unicode_literals
from ..image import Rescale


def test_Rescale_inputs():
input_map = dict(
ignore_exception=dict(
deprecated='1.0.0',
nohash=True,
usedefault=True,
),
in_file=dict(mandatory=True, ),
invert=dict(),
percentile=dict(usedefault=True, ),
ref_file=dict(mandatory=True, ),
)
inputs = Rescale.input_spec()

for key, metadata in list(input_map.items()):
for metakey, value in list(metadata.items()):
assert getattr(inputs.traits()[key], metakey) == value
def test_Rescale_outputs():
output_map = dict(out_file=dict(), )
outputs = Rescale.output_spec()

for key, metadata in list(output_map.items()):
for metakey, value in list(metadata.items()):
assert getattr(outputs.traits()[key], metakey) == value