1
# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
2
# vi: set ft=python sts=4 ts=4 sw=4 et:
3
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
4
#
5
#   See COPYING file distributed along with the NiBabel package for the
6
#   copyright and license terms.
7
#
8
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
9 33
""" Utilities for calculating and applying affine orientations """
10

11

12 33
import numpy as np
13 33
import numpy.linalg as npl
14

15 33
from .deprecated import deprecate_with_version
16

17

18 33
class OrientationError(Exception):
19 33
    pass
20

21

22 33
def io_orientation(affine, tol=None):
23
    """ Orientation of input axes in terms of output axes for `affine`
24

25
    Valid for an affine transformation from ``p`` dimensions to ``q``
26
    dimensions (``affine.shape == (q + 1, p + 1)``).
27

28
    The calculated orientations can be used to transform associated
29
    arrays to best match the output orientations. If ``p`` > ``q``, then
30
    some of the output axes should be considered dropped in this
31
    orientation.
32

33
    Parameters
34
    ----------
35
    affine : (q+1, p+1) ndarray-like
36
       Transformation affine from ``p`` inputs to ``q`` outputs.  Usually this
37
       will be a shape (4,4) matrix, transforming 3 inputs to 3 outputs, but
38
       the code also handles the more general case
39
    tol : {None, float}, optional
40
       threshold below which SVD values of the affine are considered zero. If
41
       `tol` is None, and ``S`` is an array with singular values for `affine`,
42
       and ``eps`` is the epsilon value for datatype of ``S``, then `tol` set
43
       to ``S.max() * max((q, p)) * eps``
44

45
    Returns
46
    -------
47
    orientations : (p, 2) ndarray
48
       one row per input axis, where the first value in each row is the closest
49
       corresponding output axis. The second value in each row is 1 if the
50
       input axis is in the same direction as the corresponding output axis and
51
       -1 if it is in the opposite direction.  If a row is [np.nan, np.nan],
52
       which can happen when p > q, then this row should be considered dropped.
53
    """
54 33
    affine = np.asarray(affine)
55 33
    q, p = affine.shape[0] - 1, affine.shape[1] - 1
56
    # extract the underlying rotation, zoom, shear matrix
57 33
    RZS = affine[:q, :p]
58 33
    zooms = np.sqrt(np.sum(RZS * RZS, axis=0))
59
    # Zooms can be zero, in which case all elements in the column are zero, and
60
    # we can leave them as they are
61 33
    zooms[zooms == 0] = 1
62 33
    RS = RZS / zooms
63
    # Transform below is polar decomposition, returning the closest
64
    # shearless matrix R to RS
65 33
    P, S, Qs = npl.svd(RS, full_matrices=False)
66
    # Threshold the singular values to determine the rank.
67 33
    if tol is None:
68 33
        tol = S.max() * max(RS.shape) * np.finfo(S.dtype).eps
69 33
    keep = (S > tol)
70 33
    R = np.dot(P[:, keep], Qs[keep])
71
    # the matrix R is such that np.dot(R,R.T) is projection onto the
72
    # columns of P[:,keep] and np.dot(R.T,R) is projection onto the rows
73
    # of Qs[keep].  R (== np.dot(R, np.eye(p))) gives rotation of the
74
    # unit input vectors to output coordinates.  Therefore, the row
75
    # index of abs max R[:,N], is the output axis changing most as input
76
    # axis N changes.  In case there are ties, we choose the axes
77
    # iteratively, removing used axes from consideration as we go
78 33
    ornt = np.ones((p, 2), dtype=np.int8) * np.nan
79 33
    for in_ax in range(p):
80 33
        col = R[:, in_ax]
81 33
        if not np.allclose(col, 0):
82 33
            out_ax = np.argmax(np.abs(col))
83 33
            ornt[in_ax, 0] = out_ax
84 33
            assert col[out_ax] != 0
85 33
            if col[out_ax] < 0:
86 33
                ornt[in_ax, 1] = -1
87
            else:
88 33
                ornt[in_ax, 1] = 1
89
            # remove the identified axis from further consideration, by
90
            # zeroing out the corresponding row in R
91 33
            R[out_ax, :] = 0
92 33
    return ornt
93

94

95 33
def ornt_transform(start_ornt, end_ornt):
96
    """Return the orientation that transforms from `start_ornt` to `end_ornt`.
97

98
    Parameters
99
    ----------
100
    start_ornt : (n,2) orientation array
101
        Initial orientation.
102

103
    end_ornt : (n,2) orientation array
104
        Final orientation.
105

106
    Returns
107
    -------
108
    orientations : (p, 2) ndarray
109
       The orientation that will transform the `start_ornt` to the `end_ornt`.
110
    """
111 33
    start_ornt = np.asarray(start_ornt)
112 33
    end_ornt = np.asarray(end_ornt)
113 33
    if start_ornt.shape != end_ornt.shape:
114 33
        raise ValueError("The orientations must have the same shape")
115 33
    if start_ornt.shape[1] != 2:
116 33
        raise ValueError(f"Invalid shape for an orientation: {start_ornt.shape}")
117 33
    result = np.empty_like(start_ornt)
118 33
    for end_in_idx, (end_out_idx, end_flip) in enumerate(end_ornt):
119 33
        for start_in_idx, (start_out_idx, start_flip) in enumerate(start_ornt):
120 33
            if end_out_idx == start_out_idx:
121 33
                if start_flip == end_flip:
122 33
                    flip = 1
123
                else:
124 33
                    flip = -1
125 33
                result[start_in_idx, :] = [end_in_idx, flip]
126 33
                break
127
        else:
128 33
            raise ValueError("Unable to find out axis %d in start_ornt" %
129
                             end_out_idx)
130 33
    return result
131

132

133 33
def apply_orientation(arr, ornt):
134
    """ Apply transformations implied by `ornt` to the first
135
    n axes of the array `arr`
136

137
    Parameters
138
    ----------
139
    arr : array-like of data with ndim >= n
140
    ornt : (n,2) orientation array
141
       orientation transform. ``ornt[N,1]` is flip of axis N of the
142
       array implied by `shape`, where 1 means no flip and -1 means
143
       flip.  For example, if ``N==0`` and ``ornt[0,1] == -1``, and
144
       there's an array ``arr`` of shape `shape`, the flip would
145
       correspond to the effect of ``np.flipud(arr)``.  ``ornt[:,0]`` is
146
       the transpose that needs to be done to the implied array, as in
147
       ``arr.transpose(ornt[:,0])``
148

149
    Returns
150
    -------
151
    t_arr : ndarray
152
       data array `arr` transformed according to ornt
153
    """
154 33
    t_arr = np.asarray(arr)
155 33
    ornt = np.asarray(ornt)
156 33
    n = ornt.shape[0]
157 33
    if t_arr.ndim < n:
158 33
        raise OrientationError('Data array has fewer dimensions than '
159
                               'orientation')
160
    # no coordinates can be dropped for applying the orientations
161 33
    if np.any(np.isnan(ornt[:, 0])):
162 33
        raise OrientationError('Cannot drop coordinates when '
163
                               'applying orientation to data')
164
    # apply ornt transformations
165 33
    for ax, flip in enumerate(ornt[:, 1]):
166 33
        if flip == -1:
167 33
            t_arr = np.flip(t_arr, axis=ax)
168 33
    full_transpose = np.arange(t_arr.ndim)
169
    # ornt indicates the transpose that has occurred - we reverse it
170 33
    full_transpose[:n] = np.argsort(ornt[:, 0])
171 33
    t_arr = t_arr.transpose(full_transpose)
172 33
    return t_arr
173

174

175 33
def inv_ornt_aff(ornt, shape):
176
    """ Affine transform reversing transforms implied in `ornt`
177

178
    Imagine you have an array ``arr`` of shape `shape`, and you apply the
179
    transforms implied by `ornt` (more below), to get ``tarr``.
180
    ``tarr`` may have a different shape ``shape_prime``.  This routine
181
    returns the affine that will take a array coordinate for ``tarr``
182
    and give you the corresponding array coordinate in ``arr``.
183

184
    Parameters
185
    ----------
186
    ornt : (p, 2) ndarray
187
       orientation transform. ``ornt[P, 1]` is flip of axis N of the array
188
       implied by `shape`, where 1 means no flip and -1 means flip.  For
189
       example, if ``P==0`` and ``ornt[0, 1] == -1``, and there's an array
190
       ``arr`` of shape `shape`, the flip would correspond to the effect of
191
       ``np.flipud(arr)``.  ``ornt[:,0]`` gives us the (reverse of the)
192
       transpose that has been done to ``arr``.  If there are any NaNs in
193
       `ornt`, we raise an ``OrientationError`` (see notes)
194
    shape : length p sequence
195
       shape of array you may transform with `ornt`
196

197
    Returns
198
    -------
199
    transform_affine : (p + 1, p + 1) ndarray
200
       An array ``arr`` (shape `shape`) might be transformed according to
201
       `ornt`, resulting in a transformed array ``tarr``.  `transformed_affine`
202
       is the transform that takes you from array coordinates in ``tarr`` to
203
       array coordinates in ``arr``.
204

205
    Notes
206
    -----
207
    If a row in `ornt` contains NaN, this means that the input row does not
208
    influence the output space, and is thus effectively dropped from the output
209
    space.  In that case one ``tarr`` coordinate maps to many ``arr``
210
    coordinates, we can't invert the transform, and we raise an error
211
    """
212 33
    ornt = np.asarray(ornt)
213 33
    if np.any(np.isnan(ornt)):
214 33
        raise OrientationError("We cannot invert orientation transform")
215 33
    p = ornt.shape[0]
216 33
    shape = np.array(shape)[:p]
217
    # ornt implies a flip, followed by a transpose.   We need the affine
218
    # that inverts these.  Thus we need the affine that first undoes the
219
    # effect of the transpose, then undoes the effects of the flip.
220
    # ornt indicates the transpose that has occurred to get the current
221
    # ordering, relative to canonical, so we just use that.
222
    # undo_reorder is a row permutatation matrix
223 33
    axis_transpose = [int(v) for v in ornt[:, 0]]
224 33
    undo_reorder = np.eye(p + 1)[axis_transpose + [p], :]
225 33
    undo_flip = np.diag(list(ornt[:, 1]) + [1.0])
226 33
    center_trans = -(shape - 1) / 2.0
227 33
    undo_flip[:p, p] = (ornt[:, 1] * center_trans) - center_trans
228 33
    return np.dot(undo_flip, undo_reorder)
229

230

231 33
@deprecate_with_version('orientation_affine deprecated. '
232
                        'Please use inv_ornt_aff instead.',
233
                        '3.0',
234
                        '4.0')
235 5
def orientation_affine(ornt, shape):
236 33
    return inv_ornt_aff(ornt, shape)
237

238

239 33
@deprecate_with_version('flip_axis is deprecated. '
240
                        'Please use numpy.flip instead.',
241
                        '3.2',
242
                        '5.0')
243 33
def flip_axis(arr, axis=0):
244
    """ Flip contents of `axis` in array `arr`
245

246
    Equivalent to ``np.flip(arr, axis)``.
247

248
    Parameters
249
    ----------
250
    arr : array-like
251
    axis : int, optional
252
       axis to flip.  Default `axis` == 0
253

254
    Returns
255
    -------
256
    farr : array
257
       Array with axis `axis` flipped
258
    """
259 33
    return np.flip(arr, axis)
260

261

262 33
def ornt2axcodes(ornt, labels=None):
263
    """ Convert orientation `ornt` to labels for axis directions
264

265
    Parameters
266
    ----------
267
    ornt : (N,2) array-like
268
        orientation array - see io_orientation docstring
269
    labels : optional, None or sequence of (2,) sequences
270
        (2,) sequences are labels for (beginning, end) of output axis.  That
271
        is, if the first row in `ornt` is ``[1, 1]``, and the second (2,)
272
        sequence in `labels` is ('back', 'front') then the first returned axis
273
        code will be ``'front'``.  If the first row in `ornt` had been
274
        ``[1, -1]`` then the first returned value would have been ``'back'``.
275
        If None, equivalent to ``(('L','R'),('P','A'),('I','S'))`` - that is -
276
        RAS axes.
277

278
    Returns
279
    -------
280
    axcodes : (N,) tuple
281
        labels for positive end of voxel axes.  Dropped axes get a label of
282
        None.
283

284
    Examples
285
    --------
286
    >>> ornt2axcodes([[1, 1],[0,-1],[2,1]], (('L','R'),('B','F'),('D','U')))
287
    ('F', 'L', 'U')
288
    """
289 33
    if labels is None:
290 33
        labels = list(zip('LPI', 'RAS'))
291 33
    axcodes = []
292 33
    for axno, direction in np.asarray(ornt):
293 33
        if np.isnan(axno):
294 33
            axcodes.append(None)
295 33
            continue
296 33
        axint = int(np.round(axno))
297 33
        if axint != axno:
298 33
            raise ValueError(f'Non integer axis number {axno:f}')
299 33
        elif direction == 1:
300 33
            axcode = labels[axint][1]
301 33
        elif direction == -1:
302 33
            axcode = labels[axint][0]
303
        else:
304 33
            raise ValueError('Direction should be -1 or 1')
305 33
        axcodes.append(axcode)
306 33
    return tuple(axcodes)
307

308

309 33
def axcodes2ornt(axcodes, labels=None):
310
    """ Convert axis codes `axcodes` to an orientation
311

312
    Parameters
313
    ----------
314
    axcodes : (N,) tuple
315
        axis codes - see ornt2axcodes docstring
316
    labels : optional, None or sequence of (2,) sequences
317
        (2,) sequences are labels for (beginning, end) of output axis.  That
318
        is, if the first element in `axcodes` is ``front``, and the second
319
        (2,) sequence in `labels` is ('back', 'front') then the first
320
        row of `ornt` will be ``[1, 1]``. If None, equivalent to
321
        ``(('L','R'),('P','A'),('I','S'))`` - that is - RAS axes.
322

323
    Returns
324
    -------
325
    ornt : (N,2) array-like
326
        orientation array - see io_orientation docstring
327

328
    Examples
329
    --------
330
    >>> axcodes2ornt(('F', 'L', 'U'), (('L','R'),('B','F'),('D','U')))
331
    array([[ 1.,  1.],
332
           [ 0., -1.],
333
           [ 2.,  1.]])
334
    """
335 33
    labels = list(zip('LPI', 'RAS')) if labels is None else labels
336 33
    allowed_labels = sum([list(L) for L in labels], []) + [None]
337 33
    if len(allowed_labels) != len(set(allowed_labels)):
338 33
        raise ValueError(f'Duplicate labels in {allowed_labels}')
339 33
    if not set(axcodes).issubset(allowed_labels):
340 33
        raise ValueError(f'Not all axis codes {list(axcodes)} in label set {allowed_labels}')
341 33
    n_axes = len(axcodes)
342 33
    ornt = np.ones((n_axes, 2), dtype=np.int8) * np.nan
343 33
    for code_idx, code in enumerate(axcodes):
344 33
        for label_idx, codes in enumerate(labels):
345 33
            if code is None:
346 33
                continue
347 33
            if code in codes:
348 33
                if code == codes[0]:
349 33
                    ornt[code_idx, :] = [label_idx, -1]
350
                else:
351 33
                    ornt[code_idx, :] = [label_idx, 1]
352 33
                break
353 33
    return ornt
354

355

356 33
def aff2axcodes(aff, labels=None, tol=None):
357
    """ axis direction codes for affine `aff`
358

359
    Parameters
360
    ----------
361
    aff : (N,M) array-like
362
        affine transformation matrix
363
    labels : optional, None or sequence of (2,) sequences
364
        Labels for negative and positive ends of output axes of `aff`.  See
365
        docstring for ``ornt2axcodes`` for more detail
366
    tol : None or float
367
        Tolerance for SVD of affine - see ``io_orientation`` for more detail.
368

369
    Returns
370
    -------
371
    axcodes : (N,) tuple
372
        labels for positive end of voxel axes.  Dropped axes get a label of
373
        None.
374

375
    Examples
376
    --------
377
    >>> aff = [[0,1,0,10],[-1,0,0,20],[0,0,1,30],[0,0,0,1]]
378
    >>> aff2axcodes(aff, (('L','R'),('B','F'),('D','U')))
379
    ('B', 'R', 'U')
380
    """
381 33
    ornt = io_orientation(aff, tol)
382 33
    return ornt2axcodes(ornt, labels)

Read our documentation on viewing source code .

Loading