diff --git a/nibabel/affines.py b/nibabel/affines.py index c2b2a3b1d0..87bec4be07 100644 --- a/nibabel/affines.py +++ b/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) diff --git a/nibabel/cmdline/conform.py b/nibabel/cmdline/conform.py new file mode 100644 index 0000000000..65b4ccc388 --- /dev/null +++ b/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) diff --git a/nibabel/cmdline/tests/test_conform.py b/nibabel/cmdline/tests/test_conform.py new file mode 100644 index 0000000000..fd29cbf5a2 --- /dev/null +++ b/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) diff --git a/nibabel/processing.py b/nibabel/processing.py index 0c5f921d87..b3bd83d706 100644 --- a/nibabel/processing.py +++ b/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 diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index 6fd2f59fab..8013fdc2eb 100644 --- a/nibabel/tests/test_affines.py +++ b/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) diff --git a/nibabel/tests/test_processing.py b/nibabel/tests/test_processing.py index 1e9e94091e..64139bab43 100644 --- a/nibabel/tests/test_processing.py +++ b/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)) diff --git a/setup.cfg b/setup.cfg index 3d26facc08..824a8bb5c3 100644 --- a/setup.cfg +++ b/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