1 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-  2 # vi: set ft=python sts=4 ts=4 sw=4 et:  3 33 """ Utility routines for working with points and affine transforms  4 """  5 33 import numpy as np  6 7 33 from functools import reduce  8 9 10 33 class AffineError(ValueError):  11  """ Errors in calculating or using affines """  12  # Inherits from ValueError to keep compatibility with ValueError previously  13  # raised in append_diag  14 33  pass  15 16 17 33 def apply_affine(aff, pts):  18  """ Apply affine matrix aff to points pts  19 20  Returns result of application of aff to the *right* of pts. The  21  coordinate dimension of pts should be the last.  22 23  For the 3D case, aff will be shape (4,4) and pts will have final axis  24  length 3 - maybe it will just be N by 3. The return value is the  25  transformed points, in this case::  26 27  res = np.dot(aff[:3,:3], pts.T) + aff[:3,3:4]  28  transformed_pts = res.T  29 30  This routine is more general than 3D, in that aff can have any shape  31  (N,N), and pts can have any shape, as long as the last dimension is for  32  the coordinates, and is therefore length N-1.  33 34  Parameters  35  ----------  36  aff : (N, N) array-like  37  Homogenous affine, for 3D points, will be 4 by 4. Contrary to first  38  appearance, the affine will be applied on the left of pts.  39  pts : (..., N-1) array-like  40  Points, where the last dimension contains the coordinates of each  41  point. For 3D, the last dimension will be length 3.  42 43  Returns  44  -------  45  transformed_pts : (..., N-1) array  46  transformed points  47 48  Examples  49  --------  50  >>> aff = np.array([[0,2,0,10],[3,0,0,11],[0,0,4,12],[0,0,0,1]])  51  >>> pts = np.array([[1,2,3],[2,3,4],[4,5,6],[6,7,8]])  52  >>> apply_affine(aff, pts) #doctest: +ELLIPSIS  53  array([[14, 14, 24],  54  [16, 17, 28],  55  [20, 23, 36],  56  [24, 29, 44]]...)  57 58  Just to show that in the simple 3D case, it is equivalent to:  59 60  >>> (np.dot(aff[:3,:3], pts.T) + aff[:3,3:4]).T #doctest: +ELLIPSIS  61  array([[14, 14, 24],  62  [16, 17, 28],  63  [20, 23, 36],  64  [24, 29, 44]]...)  65 66  But pts can be a more complicated shape:  67 68  >>> pts = pts.reshape((2,2,3))  69  >>> apply_affine(aff, pts) #doctest: +ELLIPSIS  70  array([[[14, 14, 24],  71  [16, 17, 28]],  72   73  [[20, 23, 36],  74  [24, 29, 44]]]...)  75  """  76 33  aff = np.asarray(aff)  77 33  pts = np.asarray(pts)  78 33  shape = pts.shape  79 33  pts = pts.reshape((-1, shape[-1]))  80  # rzs == rotations, zooms, shears  81 33  rzs = aff[:-1, :-1]  82 33  trans = aff[:-1, -1]  83 33  res = np.dot(pts, rzs.T) + trans[None, :]  84 33  return res.reshape(shape)  85 86 87 33 def to_matvec(transform):  88  """Split a transform into its matrix and vector components.  89 90  The tranformation must be represented in homogeneous coordinates and is  91  split into its rotation matrix and translation vector components.  92 93  Parameters  94  ----------  95  transform : array-like  96  NxM transform matrix in homogeneous coordinates representing an affine  97  transformation from an (N-1)-dimensional space to an (M-1)-dimensional  98  space. An example is a 4x4 transform representing rotations and  99  translations in 3 dimensions. A 4x3 matrix can represent a  100  2-dimensional plane embedded in 3 dimensional space.  101 102  Returns  103  -------  104  matrix : (N-1, M-1) array  105  Matrix component of transform  106  vector : (M-1,) array  107  Vector compoent of transform  108 109  See Also  110  --------  111  from_matvec  112 113  Examples  114  --------  115  >>> aff = np.diag([2, 3, 4, 1])  116  >>> aff[:3,3] = [9, 10, 11]  117  >>> to_matvec(aff)  118  (array([[2, 0, 0],  119  [0, 3, 0],  120  [0, 0, 4]]), array([ 9, 10, 11]))  121  """  122 33  transform = np.asarray(transform)  123 33  ndimin = transform.shape[0] - 1  124 33  ndimout = transform.shape[1] - 1  125 33  matrix = transform[0:ndimin, 0:ndimout]  126 33  vector = transform[0:ndimin, ndimout]  127 33  return matrix, vector  128 129 130 33 def from_matvec(matrix, vector=None):  131  """ Combine a matrix and vector into an homogeneous affine  132 133  Combine a rotation / scaling / shearing matrix and translation vector into  134  a transform in homogeneous coordinates.  135 136  Parameters  137  ----------  138  matrix : array-like  139  An NxM array representing the the linear part of the transform.  140  A transform from an M-dimensional space to an N-dimensional space.  141  vector : None or array-like, optional  142  None or an (N,) array representing the translation. None corresponds to  143  an (N,) array of zeros.  144 145  Returns  146  -------  147  xform : array  148  An (N+1, M+1) homogenous transform matrix.  149 150  See Also  151  --------  152  to_matvec  153 154  Examples  155  --------  156  >>> from_matvec(np.diag([2, 3, 4]), [9, 10, 11])  157  array([[ 2, 0, 0, 9],  158  [ 0, 3, 0, 10],  159  [ 0, 0, 4, 11],  160  [ 0, 0, 0, 1]])  161 162  The vector argument is optional:  163 164  >>> from_matvec(np.diag([2, 3, 4]))  165  array([[2, 0, 0, 0],  166  [0, 3, 0, 0],  167  [0, 0, 4, 0],  168  [0, 0, 0, 1]])  169  """  170 33  matrix = np.asarray(matrix)  171 33  nin, nout = matrix.shape  172 33  t = np.zeros((nin + 1, nout + 1), matrix.dtype)  173 33  t[0:nin, 0:nout] = matrix  174 33  t[nin, nout] = 1.  175 33  if vector is not None:  176 33  t[0:nin, nout] = vector  177 33  return t  178 179 180 33 def append_diag(aff, steps, starts=()):  181  """ Add diagonal elements steps and translations starts to affine  182 183  Typical use is in expanding 4x4 affines to larger dimensions. Nipy is the  184  main consumer because it uses NxM affines, whereas we generally only use  185  4x4 affines; the routine is here for convenience.  186 187  Parameters  188  ----------  189  aff : 2D array  190  N by M affine matrix  191  steps : scalar or sequence  192  diagonal elements to append.  193  starts : scalar or sequence  194  elements to append to last column of aff, representing translations  195  corresponding to the steps. If empty, expands to a vector of zeros  196  of the same length as steps  197 198  Returns  199  -------  200  aff_plus : 2D array  201  Now P by Q where L = len(steps) and P == N+L, Q=N+L  202 203  Examples  204  --------  205  >>> aff = np.eye(4)  206  >>> aff[:3,:3] = np.arange(9).reshape((3,3))  207  >>> append_diag(aff, [9, 10], [99,100])  208  array([[ 0., 1., 2., 0., 0., 0.],  209  [ 3., 4., 5., 0., 0., 0.],  210  [ 6., 7., 8., 0., 0., 0.],  211  [ 0., 0., 0., 9., 0., 99.],  212  [ 0., 0., 0., 0., 10., 100.],  213  [ 0., 0., 0., 0., 0., 1.]])  214  """  215 33  aff = np.asarray(aff)  216 33  steps = np.atleast_1d(steps)  217 33  starts = np.atleast_1d(starts)  218 33  n_steps = len(steps)  219 33  if len(starts) == 0:  220 33  starts = np.zeros(n_steps, dtype=steps.dtype)  221 33  elif len(starts) != n_steps:  222 33  raise AffineError('Steps should have same length as starts')  223 33  old_n_out, old_n_in = aff.shape[0] - 1, aff.shape[1] - 1  224  # make new affine  225 33  aff_plus = np.zeros((old_n_out + n_steps + 1,  226  old_n_in + n_steps + 1), dtype=aff.dtype)  227  # Get stuff from old affine  228 33  aff_plus[:old_n_out, :old_n_in] = aff[:old_n_out, :old_n_in]  229 33  aff_plus[:old_n_out, -1] = aff[:old_n_out, -1]  230  # Add new diagonal elements  231 33  for i, el in enumerate(steps):  232 33  aff_plus[old_n_out + i, old_n_in + i] = el  233  # Add translations for new affine, plus last 1  234 33  aff_plus[old_n_out:, -1] = list(starts) + [1]  235 33  return aff_plus  236 237 238 33 def dot_reduce(*args):  239  r""" Apply numpy dot product function from right to left on arrays  240 241  For passed arrays :math:A, B, C, ... Z returns :math:A \dot B \dot C ...  242  \dot Z where "." is the numpy array dot product.  243 244  Parameters  245  ----------  246  \*\*args : arrays  247  Arrays that can be passed to numpy dot function  248 249  Returns  250  -------  251  dot_product : array  252  If there are N arguments, result of arg[0].dot(arg[1].dot(arg[2].dot  253  ... arg[N-2].dot(arg[N-1])))...  254  """  255 33  return reduce(lambda x, y: np.dot(y, x), args[::-1])  256 257 258 33 def voxel_sizes(affine):  259  r""" Return voxel size for each input axis given affine  260 261  The affine is the mapping between array (voxel) coordinates and mm  262  (world) coordinates.  263 264  The voxel size for the first voxel (array) axis is the distance moved in  265  world coordinates when moving one unit along the first voxel (array) axis.  266  This is the distance between the world coordinate of voxel (0, 0, 0) and  267  the world coordinate of voxel (1, 0, 0). The world coordinate vector of  268  voxel coordinate vector (0, 0, 0) is given by v0 = affine.dot((0, 0, 0,  269  1)[:3]. The world coordinate vector of voxel vector (1, 0, 0) is  270  v1_ax1 = affine.dot((1, 0, 0, 1))[:3]. The final 1 in the voxel  271  vectors and the [:3] at the end are because the affine works on  272  homogenous coodinates. The translations part of the affine is trans =  273  affine[:3, 3], and the rotations, zooms and shearing part of the affine  274  is rzs = affine[:3, :3]. Because of the final 1 in the input voxel  275  vector, v0 == rzs.dot((0, 0, 0)) + trans, and v1_ax1 == rzs.dot((1,  276  0, 0)) + trans, and the difference vector is rzs.dot((0, 0, 0)) -  277  rzs.dot((1, 0, 0)) == rzs.dot((1, 0, 0)) == rzs[:, 0]. The distance  278  vectors in world coordinates between (0, 0, 0) and (1, 0, 0), (0, 1, 0),  279  (0, 0, 1) are given by rzs.dot(np.eye(3)) = rzs. The voxel sizes are  280  the Euclidean lengths of the distance vectors. So, the voxel sizes are  281  the Euclidean lengths of the columns of the affine (excluding the last row  282  and column of the affine).  283 284  Parameters  285  ----------  286  affine : 2D array-like  287  Affine transformation array. Usually shape (4, 4), but can be any 2D  288  array.  289 290  Returns  291  -------  292  vox_sizes : 1D array  293  Voxel sizes for each input axis of affine. Usually 1D array length 3,  294  but in general has length (N-1) where input affine is shape (M, N).  295  """  296 33  top_left = affine[:-1, :-1]  297 33  return np.sqrt(np.sum(top_left ** 2, axis=0))  298 299 300 33 def obliquity(affine):  301  r"""  302  Estimate the *obliquity* an affine's axes represent.  303 304  The term *obliquity* is defined here as the rotation of those axes with  305  respect to the cardinal axes.  306  This implementation is inspired by AFNI's implementation  307  _.  308  For further details about *obliquity*, check AFNI's documentation  309  _.  310 311  Parameters  312  ----------  313  affine : 2D array-like  314  Affine transformation array. Usually shape (4, 4), but can be any 2D  315  array.  316 317  Returns  318  -------  319  angles : 1D array-like  320  The *obliquity* of each axis with respect to the cardinal axes, in radians.  321 322  """  323 33  vs = voxel_sizes(affine)  324 33  best_cosines = np.abs(affine[:-1, :-1] / vs).max(axis=1)  325 33  return np.arccos(best_cosines)  326 327 328 33 def rescale_affine(affine, shape, zooms, new_shape=None):  329  """ Return a new affine matrix with updated voxel sizes (zooms)  330 331  This function preserves the rotations and shears of the original  332  affine, as well as the RAS location of the central voxel of the  333  image.  334 335  Parameters  336  ----------  337  affine : (N, N) array-like  338  NxN transform matrix in homogeneous coordinates representing an affine  339  transformation from an (N-1)-dimensional space to an (N-1)-dimensional  340  space. An example is a 4x4 transform representing rotations and  341  translations in 3 dimensions.  342  shape : (N-1,) array-like  343  The extent of the (N-1) dimensions of the original space  344  zooms : (N-1,) array-like  345  The size of voxels of the output affine  346  new_shape : (N-1,) array-like, optional  347  The extent of the (N-1) dimensions of the space described by the  348  new affine. If None, use shape.  349 350  Returns  351  -------  352  affine : (N, N) array  353  A new affine transform with the specified voxel sizes  354 355  """  356 33  shape = np.array(shape, copy=False)  357 33  new_shape = np.array(new_shape if new_shape is not None else shape)  358 359 33  s = voxel_sizes(affine)  360 33  rzs_out = affine[:3, :3] * zooms / s  361 362  # Using xyz = A @ ijk, determine translation  363 33  centroid = apply_affine(affine, (shape - 1) // 2)  364 33  t_out = centroid - rzs_out @ ((new_shape - 1) // 2)  365 33  return from_matvec(rzs_out, t_out) 

Read our documentation on viewing source code .