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
    <BLANKLINE>
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
    <https://github.com/afni/afni/blob/b6a9f7a21c1f3231ff09efbd861f8975ad48e525/src/thd_coords.c#L660-L698>`_.
308
    For further details about *obliquity*, check `AFNI's documentation
309
    <https://sscc.nimh.nih.gov/sscc/dglen/Obliquity>`_.
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 .

Loading