Merge pull request #1 from effigies/add/fs-conform

RF: Update conformation to reorient, rescale and resample

Author: kaczmarj

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/processing.py

@@ -21,7 +21,7 @@
 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
@@ -373,19 +373,18 @@ def conform(from_img,
     elif len(voxel_size) != required_ndim:
         raise ValueError("`voxel_size` must have {} values".format(required_ndim))
 
-    # Create template image to which input is resampled.
-    tmpl_hdr = from_img.header_class().from_header(from_img.header)
-    tmpl_hdr.set_data_shape(out_shape)
-    tmpl_hdr.set_zooms(voxel_size)
-    tmpl = from_img.__class__(np.empty(out_shape), affine=np.eye(4), header=tmpl_hdr)
+    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=tmpl, order=order, mode="constant",
+        from_img=from_img, to_vox_map=(out_shape, out_aff), order=order, mode="constant",
         cval=cval, out_class=out_class)
 
-    # Reorient to desired orientation.
-    start_ornt = io_orientation(out_img.affine)
-    end_ornt = axcodes2ornt(orientation)
-    transform = ornt_transform(start_ornt, end_ornt)
-    return out_img.as_reoriented(transform)
+    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)