Merge pull request #853 from kaczmarj/add/fs-conform

NF: Conformation function and CLI tool to apply shape, orientation and zooms

Author: effigies

nibabel/affines.py

@@ -323,3 +323,43 @@ def obliquity(affine):
     vs = voxel_sizes(affine)
     best_cosines = np.abs(affine[:-1, :-1] / vs).max(axis=1)
     return np.arccos(best_cosines)
+
+
+def rescale_affine(affine, shape, zooms, new_shape=None):
+    """ Return a new affine matrix with updated voxel sizes (zooms)
+
+    This function preserves the rotations and shears of the original
+    affine, as well as the RAS location of the central voxel of the
+    image.
+
+    Parameters
+    ----------
+    affine : (N, N) array-like
+        NxN transform matrix in homogeneous coordinates representing an affine
+        transformation from an (N-1)-dimensional space to an (N-1)-dimensional
+        space. An example is a 4x4 transform representing rotations and
+        translations in 3 dimensions.
+    shape : (N-1,) array-like
+        The extent of the (N-1) dimensions of the original space
+    zooms : (N-1,) array-like
+        The size of voxels of the output affine
+    new_shape : (N-1,) array-like, optional
+        The extent of the (N-1) dimensions of the space described by the
+        new affine. If ``None``, use ``shape``.
+
+    Returns
+    -------
+    affine : (N, N) array
+        A new affine transform with the specified voxel sizes
+
+    """
+    shape = np.array(shape, copy=False)
+    new_shape = np.array(new_shape if new_shape is not None else shape)
+
+    s = voxel_sizes(affine)
+    rzs_out = affine[:3, :3] * zooms / s
+
+    # Using xyz = A @ ijk, determine translation
+    centroid = apply_affine(affine, (shape - 1) // 2)
+    t_out = centroid - rzs_out @ ((new_shape - 1) // 2)
+    return from_matvec(rzs_out, t_out)

nibabel/cmdline/conform.py

@@ -0,0 +1,59 @@
+#!python
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+"""
+Conform neuroimaging volume to arbitrary shape and voxel size.
+"""
+
+import argparse
+from pathlib import Path
+
+from nibabel import __version__
+from nibabel.loadsave import load, save
+from nibabel.processing import conform
+
+
+def _get_parser():
+    """Return command-line argument parser."""
+    p = argparse.ArgumentParser(description=__doc__)
+    p.add_argument("infile",
+                   help="Neuroimaging volume to conform.")
+    p.add_argument("outfile",
+                   help="Name of output file.")
+    p.add_argument("--out-shape", nargs=3, default=(256, 256, 256), type=int,
+                   help="Shape of the conformed output.")
+    p.add_argument("--voxel-size", nargs=3, default=(1, 1, 1), type=int,
+                   help="Voxel size in millimeters of the conformed output.")
+    p.add_argument("--orientation", default="RAS",
+                   help="Orientation of the conformed output.")
+    p.add_argument("-f", "--force", action="store_true",
+                   help="Overwrite existing output files.")
+    p.add_argument("-V", "--version", action="version", version="{} {}".format(p.prog, __version__))
+
+    return p
+
+
+def main(args=None):
+    """Main program function."""
+    parser = _get_parser()
+    opts = parser.parse_args(args)
+    from_img = load(opts.infile)
+
+    if not opts.force and Path(opts.outfile).exists():
+        raise FileExistsError("Output file exists: {}".format(opts.outfile))
+
+    out_img = conform(
+        from_img=from_img,
+        out_shape=opts.out_shape,
+        voxel_size=opts.voxel_size,
+        order=3,
+        cval=0.0,
+        orientation=opts.orientation)
+
+    save(out_img, opts.outfile)

nibabel/cmdline/tests/test_conform.py

@@ -0,0 +1,56 @@
+#!python
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+
+import unittest
+
+import pytest
+
+import nibabel as nib
+from nibabel.testing import test_data
+from nibabel.cmdline.conform import main
+from nibabel.optpkg import optional_package
+
+_, have_scipy, _ = optional_package('scipy.ndimage')
+needs_scipy = unittest.skipUnless(have_scipy, 'These tests need scipy')
+
+
+@needs_scipy
+def test_default(tmpdir):
+    infile = test_data(fname="anatomical.nii")
+    outfile = tmpdir / "output.nii.gz"
+    main([str(infile), str(outfile)])
+    assert outfile.isfile()
+    c = nib.load(outfile)
+    assert c.shape == (256, 256, 256)
+    assert c.header.get_zooms() == (1, 1, 1)
+    assert nib.orientations.aff2axcodes(c.affine) == ('R', 'A', 'S')
+
+    with pytest.raises(FileExistsError):
+        main([str(infile), str(outfile)])
+
+    main([str(infile), str(outfile), "--force"])
+    assert outfile.isfile()
+
+
+@needs_scipy
+def test_nondefault(tmpdir):
+    infile = test_data(fname="anatomical.nii")
+    outfile = tmpdir / "output.nii.gz"
+    out_shape = (100, 100, 150)
+    voxel_size = (1, 2, 4)
+    orientation = "LAS"
+    args = "{} {} --out-shape {} --voxel-size {} --orientation {}".format(
+        infile, outfile, " ".join(map(str, out_shape)), " ".join(map(str, voxel_size)), orientation)
+    main(args.split())
+    assert outfile.isfile()
+    c = nib.load(outfile)
+    assert c.shape == out_shape
+    assert c.header.get_zooms() == voxel_size
+    assert nib.orientations.aff2axcodes(c.affine) == tuple(orientation)

nibabel/processing.py

@@ -21,9 +21,10 @@
 from .optpkg import optional_package
 spnd, _, _ = optional_package('scipy.ndimage')
 
-from .affines import AffineError, to_matvec, from_matvec, append_diag
+from .affines import AffineError, to_matvec, from_matvec, append_diag, rescale_affine
 from .spaces import vox2out_vox
 from .nifti1 import Nifti1Image
+from .orientations import axcodes2ornt, io_orientation, ornt_transform
 from .imageclasses import spatial_axes_first
 
 SIGMA2FWHM = np.sqrt(8 * np.log(2))
@@ -310,3 +311,80 @@ def smooth_image(img,
                                    mode=mode,
                                    cval=cval)
     return out_class(sm_data, img.affine, img.header)
+
+
+def conform(from_img,
+            out_shape=(256, 256, 256),
+            voxel_size=(1.0, 1.0, 1.0),
+            order=3,
+            cval=0.0,
+            orientation='RAS',
+            out_class=None):
+    """ Resample image to ``out_shape`` with voxels of size ``voxel_size``.
+
+    Using the default arguments, this function is meant to replicate most parts
+    of FreeSurfer's ``mri_convert --conform`` command. Specifically, this
+    function:
+        - Resamples data to ``output_shape``
+        - Resamples voxel sizes to ``voxel_size``
+        - Reorients to RAS (``mri_convert --conform`` reorients to LIA)
+
+    Unlike ``mri_convert --conform``, this command does not:
+        - Transform data to range [0, 255]
+        - Cast to unsigned eight-bit integer
+
+    Parameters
+    ----------
+    from_img : object
+        Object having attributes ``dataobj``, ``affine``, ``header`` and
+        ``shape``. If `out_class` is not None, ``img.__class__`` should be able
+        to construct an image from data, affine and header.
+    out_shape : sequence, optional
+        The shape of the output volume. Default is (256, 256, 256).
+    voxel_size : sequence, optional
+        The size in millimeters of the voxels in the resampled output. Default
+        is 1mm isotropic.
+    order : int, optional
+        The order of the spline interpolation, default is 3.  The order has to
+        be in the range 0-5 (see ``scipy.ndimage.affine_transform``)
+    cval : scalar, optional
+        Value used for points outside the boundaries of the input if
+        ``mode='constant'``. Default is 0.0 (see
+        ``scipy.ndimage.affine_transform``)
+    orientation : str, optional
+        Orientation of output image. Default is "RAS".
+    out_class : None or SpatialImage class, optional
+        Class of output image.  If None, use ``from_img.__class__``.
+
+    Returns
+    -------
+    out_img : object
+        Image of instance specified by `out_class`, containing data output from
+        resampling `from_img` into axes aligned to the output space of
+        ``from_img.affine``
+    """
+    # Only support 3D images. This can be made more general in the future, once tests
+    # are written.
+    required_ndim = 3
+    if from_img.ndim != required_ndim:
+        raise ValueError("Only 3D images are supported.")
+    elif len(out_shape) != required_ndim:
+        raise ValueError("`out_shape` must have {} values".format(required_ndim))
+    elif len(voxel_size) != required_ndim:
+        raise ValueError("`voxel_size` must have {} values".format(required_ndim))
+
+    start_ornt = io_orientation(from_img.affine)
+    end_ornt = axcodes2ornt(orientation)
+    transform = ornt_transform(start_ornt, end_ornt)
+
+    # Reorient first to ensure shape matches expectations
+    reoriented = from_img.as_reoriented(transform)
+
+    out_aff = rescale_affine(reoriented.affine, reoriented.shape, voxel_size, out_shape)
+
+    # Resample input image.
+    out_img = resample_from_to(
+        from_img=from_img, to_vox_map=(out_shape, out_aff), order=order, mode="constant",
+        cval=cval, out_class=out_class)
+
+    return out_img

nibabel/tests/test_affines.py

@@ -7,7 +7,8 @@
 
 from ..eulerangles import euler2mat
 from ..affines import (AffineError, apply_affine, append_diag, to_matvec,
-                       from_matvec, dot_reduce, voxel_sizes, obliquity)
+                       from_matvec, dot_reduce, voxel_sizes, obliquity, rescale_affine)
+from ..orientations import aff2axcodes
 
 
 import pytest
@@ -192,3 +193,22 @@ def test_obliquity():
     assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0])
     assert_almost_equal(obliquity(oblique) * 180 / pi,
                         [0.0810285, 5.1569949, 5.1569376])
+
+
+def test_rescale_affine():
+    rng = np.random.RandomState(20200415)
+    orig_shape = rng.randint(low=20, high=512, size=(3,))
+    orig_aff = np.eye(4)
+    orig_aff[:3, :] = rng.normal(size=(3, 4))
+    orig_zooms = voxel_sizes(orig_aff)
+    orig_axcodes = aff2axcodes(orig_aff)
+    orig_centroid = apply_affine(orig_aff, (orig_shape - 1) // 2)
+    
+    for new_shape in (None, tuple(orig_shape), (256, 256, 256), (64, 64, 40)):
+        for new_zooms in ((1, 1, 1), (2, 2, 3), (0.5, 0.5, 0.5)):
+            new_aff = rescale_affine(orig_aff, orig_shape, new_zooms, new_shape)
+            assert aff2axcodes(new_aff) == orig_axcodes
+            if new_shape is None:
+                new_shape = tuple(orig_shape)
+            new_centroid = apply_affine(new_aff, (np.array(new_shape) - 1) // 2)
+            assert_almost_equal(new_centroid, orig_centroid)

nibabel/tests/test_processing.py

@@ -19,10 +19,11 @@
 
 import nibabel as nib
 from nibabel.processing import (sigma2fwhm, fwhm2sigma, adapt_affine,
-                                resample_from_to, resample_to_output, smooth_image)
+                                resample_from_to, resample_to_output, smooth_image,
+                                conform)
 from nibabel.nifti1 import Nifti1Image
 from nibabel.nifti2 import Nifti2Image
-from nibabel.orientations import flip_axis, inv_ornt_aff
+from nibabel.orientations import aff2axcodes, flip_axis, inv_ornt_aff
 from nibabel.affines import (AffineError, from_matvec, to_matvec, apply_affine,
                              voxel_sizes)
 from nibabel.eulerangles import euler2mat
@@ -420,3 +421,35 @@ def test_against_spm_resample():
     moved2output = resample_to_output(moved_anat, 4, order=1, cval=np.nan)
     spm2output = nib.load(pjoin(DATA_DIR, 'reoriented_anat_moved.nii'))
     assert_spm_resampling_close(moved_anat, moved2output, spm2output);
+
+
+@needs_scipy
+def test_conform():
+    anat = nib.load(pjoin(DATA_DIR, 'anatomical.nii'))
+
+    # Test with default arguments.
+    c = conform(anat)
+    assert c.shape == (256, 256, 256)
+    assert c.header.get_zooms() == (1, 1, 1)
+    assert c.dataobj.dtype.type == anat.dataobj.dtype.type
+    assert aff2axcodes(c.affine) == ('R', 'A', 'S')
+    assert isinstance(c, Nifti1Image)
+
+    # Test with non-default arguments.
+    c = conform(anat, out_shape=(100, 100, 200), voxel_size=(2, 2, 1.5),
+        orientation="LPI", out_class=Nifti2Image)
+    assert c.shape == (100, 100, 200)
+    assert c.header.get_zooms() == (2, 2, 1.5)
+    assert c.dataobj.dtype.type == anat.dataobj.dtype.type
+    assert aff2axcodes(c.affine) == ('L', 'P', 'I')
+    assert isinstance(c, Nifti2Image)
+
+    # TODO: support nD images in `conform` in the future, but for now, test that we get
+    # errors on non-3D images.
+    func = nib.load(pjoin(DATA_DIR, 'functional.nii'))
+    with pytest.raises(ValueError):
+        conform(func)
+    with pytest.raises(ValueError):
+        conform(anat, out_shape=(100, 100))
+    with pytest.raises(ValueError):
+        conform(anat, voxel_size=(2, 2))

setup.cfg

@@ -69,6 +69,7 @@ all =
 
 [options.entry_points]
 console_scripts =
+    nib-conform=nibabel.cmdline.conform:main
     nib-ls=nibabel.cmdline.ls:main
     nib-dicomfs=nibabel.cmdline.dicomfs:main
     nib-diff=nibabel.cmdline.diff:main