/usr/share/pyshared/nibabel/orientations.py is in python-nibabel 1.3.0-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 | # 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.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
''' Utilities for calculating and applying affine orientations '''
import numpy as np
import numpy.linalg as npl
class OrientationError(Exception):
pass
def io_orientation(affine, tol=None):
''' Orientation of input axes in terms of output axes for `affine`
Valid for an affine transformation from ``p`` dimensions to ``q``
dimensions (``affine.shape == (q + 1, p + 1)``).
The calculated orientations can be used to transform associated
arrays to best match the output orientations. If ``p`` > ``q``, then
some of the output axes should be considered dropped in this
orientation.
Parameters
----------
affine : (q+1, p+1) ndarray-like
Transformation affine from ``p`` inputs to ``q`` outputs. Usually this
will be a shape (4,4) matrix, transforming 3 inputs to 3 outputs, but the
code also handles the more general case
tol : {None, float}, optional
threshold below which SVD values of the affine are considered zero. If
`tol` is None, and ``S`` is an array with singular values for `affine`,
and ``eps`` is the epsilon value for datatype of ``S``, then `tol` set to
``S.max() * max((q, p)) * eps``
Returns
-------
orientations : (p, 2) ndarray
one row per input axis, where the first value in each row is the closest
corresponding output axis. The second value in each row is 1 if the input
axis is in the same direction as the corresponding output axis and -1 if
it is in the opposite direction. If a row is [np.nan, np.nan], which can
happen when p > q, then this row should be considered dropped.
'''
affine = np.asarray(affine)
q, p = affine.shape[0]-1, affine.shape[1]-1
# extract the underlying rotation, zoom, shear matrix
RZS = affine[:q, :p]
zooms = np.sqrt(np.sum(RZS * RZS, axis=0))
# Zooms can be zero, in which case all elements in the column are zero, and
# we can leave them as they are
zooms[zooms == 0] = 1
RS = RZS / zooms
# Transform below is polar decomposition, returning the closest
# shearless matrix R to RS
P, S, Qs = npl.svd(RS)
# Threshold the singular values to determine the rank.
if tol is None:
tol = S.max() * max(RS.shape) * np.finfo(S.dtype).eps
keep = (S > tol)
R = np.dot(P[:, keep], Qs[keep])
# the matrix R is such that np.dot(R,R.T) is projection onto the
# columns of P[:,keep] and np.dot(R.T,R) is projection onto the rows
# of Qs[keep]. R (== np.dot(R, np.eye(p))) gives rotation of the
# unit input vectors to output coordinates. Therefore, the row
# index of abs max R[:,N], is the output axis changing most as input
# axis N changes. In case there are ties, we choose the axes
# iteratively, removing used axes from consideration as we go
ornt = np.ones((p, 2), dtype=np.int8) * np.nan
for in_ax in range(p):
col = R[:, in_ax]
if not np.allclose(col, 0):
out_ax = np.argmax(np.abs(col))
ornt[in_ax, 0] = out_ax
assert col[out_ax] != 0
if col[out_ax] < 0:
ornt[in_ax, 1] = -1
else:
ornt[in_ax, 1] = 1
# remove the identified axis from further consideration, by
# zeroing out the corresponding row in R
R[out_ax, :] = 0
return ornt
def apply_orientation(arr, ornt):
''' Apply transformations implied by `ornt` to the first
n axes of the array `arr`
Parameters
----------
arr : array-like of data with ndim >= n
ornt : (n,2) orientation array
orientation transform. ``ornt[N,1]` is flip of axis N of the
array implied by `shape`, where 1 means no flip and -1 means
flip. For example, if ``N==0`` and ``ornt[0,1] == -1``, and
there's an array ``arr`` of shape `shape`, the flip would
correspond to the effect of ``np.flipud(arr)``. ``ornt[:,0]`` is
the transpose that needs to be done to the implied array, as in
``arr.transpose(ornt[:,0])``
Returns
-------
t_arr : ndarray
data array `arr` transformed according to ornt
'''
t_arr = np.asarray(arr)
ornt = np.asarray(ornt)
n = ornt.shape[0]
if t_arr.ndim < n:
raise OrientationError('Data array has fewer dimensions than '
'orientation')
# no coordinates can be dropped for applying the orientations
if np.any(np.isnan(ornt[:, 0])):
raise OrientationError('Cannot drop coordinates when '
'applying orientation to data')
# apply ornt transformations
for ax, flip in enumerate(ornt[:, 1]):
if flip == -1:
t_arr = flip_axis(t_arr, axis=ax)
full_transpose = np.arange(t_arr.ndim)
# ornt indicates the transpose that has occurred - we reverse it
full_transpose[:n] = np.argsort(ornt[:, 0])
t_arr = t_arr.transpose(full_transpose)
return t_arr
def inv_ornt_aff(ornt, shape):
''' Affine transform reversing transforms implied in `ornt`
Imagine you have an array ``arr`` of shape `shape`, and you apply the
transforms implied by `ornt` (more below), to get ``tarr``.
``tarr`` may have a different shape ``shape_prime``. This routine
returns the affine that will take a array coordinate for ``tarr``
and give you the corresponding array coordinate in ``arr``.
Parameters
----------
ornt : (p, 2) ndarray
orientation transform. ``ornt[P, 1]` is flip of axis N of the array
implied by `shape`, where 1 means no flip and -1 means flip. For
example, if ``P==0`` and ``ornt[0, 1] == -1``, and there's an array
``arr`` of shape `shape`, the flip would correspond to the effect of
``np.flipud(arr)``. ``ornt[:,0]`` gives us the (reverse of the)
transpose that has been done to ``arr``. If there are any NaNs in
`ornt`, we raise an ``OrientationError`` (see notes)
shape : length p sequence
shape of array you may transform with `ornt`
Returns
-------
transform_affine : (p + 1, p + 1) ndarray
An array ``arr`` (shape `shape`) might be transformed according to
`ornt`, resulting in a transformed array ``tarr``. `transformed_affine`
is the transform that takes you from array coordinates in ``tarr`` to
array coordinates in ``arr``.
Notes
-----
If a row in `ornt` contains NaN, this means that the input row does not
influence the output space, and is thus effectively dropped from the output
space. In that case one ``tarr`` coordinate maps to many ``arr``
coordinates, we can't invert the transform, and we raise an error
'''
ornt = np.asarray(ornt)
if np.any(np.isnan(ornt)):
raise OrientationError("We cannot invert orientation transform")
p = ornt.shape[0]
shape = np.array(shape)[:p]
# ornt implies a flip, followed by a transpose. We need the affine
# that inverts these. Thus we need the affine that first undoes the
# effect of the transpose, then undoes the effects of the flip.
# ornt indicates the transpose that has occurred to get the current
# ordering, relative to canonical, so we just use that.
# undo_reorder is a row permutatation matrix
undo_reorder = np.eye(p + 1)[list(ornt[:, 0]) + [p], :]
undo_flip = np.diag(list(ornt[:, 1]) + [1.0])
center_trans = -(shape - 1) / 2.0
undo_flip[:p, p] = (ornt[:, 1] * center_trans) - center_trans
return np.dot(undo_flip, undo_reorder)
@np.deprecate_with_doc("Please use inv_ornt_aff instead")
def orientation_affine(ornt, shape):
return inv_ornt_aff(ornt, shape)
def flip_axis(arr, axis=0):
''' Flip contents of `axis` in array `arr`
``flip_axis`` is the same transform as ``np.flipud``, but for any
axis. For example ``flip_axis(arr, axis=0)`` is the same transform
as ``np.flipud(arr)``, and ``flip_axis(arr, axis=1)`` is the same
transform as ``np.fliplr(arr)``
Parameters
----------
arr : array-like
axis : int, optional
axis to flip. Default `axis` == 0
Returns
-------
farr : array
Array with axis `axis` flipped
Examples
--------
>>> a = np.arange(6).reshape((2,3))
>>> a
array([[0, 1, 2],
[3, 4, 5]])
>>> flip_axis(a, axis=0)
array([[3, 4, 5],
[0, 1, 2]])
>>> flip_axis(a, axis=1)
array([[2, 1, 0],
[5, 4, 3]])
'''
arr = np.asanyarray(arr)
arr = arr.swapaxes(0, axis)
arr = np.flipud(arr)
return arr.swapaxes(axis, 0)
def ornt2axcodes(ornt, labels=None):
""" Convert orientation `ornt` to labels for axis directions
Parameters
----------
ornt : (N,2) array-like
orientation array - see io_orientation docstring
labels : optional, None or sequence of (2,) sequences
(2,) sequences are labels for (beginning, end) of output axis. That is,
if the first row in `ornt` is ``[1, 1]``, and the second (2,) sequence
in `labels` is ('back', 'front') then the first returned axis code will
be ``'front'``. If the first row in `ornt` had been ``[1, -1]`` then
the first returned value would have been ``'back'``. If None,
equivalent to ``(('L','R'),('P','A'),('I','S'))`` - that is - RAS axes.
Returns
-------
axcodes : (N,) tuple
labels for positive end of voxel axes. Dropped axes get a label of
None.
Examples
--------
>>> ornt2axcodes([[1, 1],[0,-1],[2,1]], (('L','R'),('B','F'),('D','U')))
('F', 'L', 'U')
"""
if labels is None:
labels = zip('LPI', 'RAS')
axcodes = []
for axno, direction in np.asarray(ornt):
if np.isnan(axno):
axcodes.append(None)
continue
axint = int(np.round(axno))
if axint != axno:
raise ValueError('Non integer axis number %f' % axno)
elif direction == 1:
axcode = labels[axint][1]
elif direction == -1:
axcode = labels[axint][0]
else:
raise ValueError('Direction should be -1 or 1')
axcodes.append(axcode)
return tuple(axcodes)
def aff2axcodes(aff, labels=None, tol=None):
""" axis direction codes for affine `aff`
Parameters
----------
aff : (N,M) array-like
affine transformation matrix
labels : optional, None or sequence of (2,) sequences
Labels for negative and positive ends of output axes of `aff`. See
docstring for ``ornt2axcodes`` for more detail
tol : None or float
Tolerance for SVD of affine - see ``io_orientation`` for more detail.
Returns
-------
axcodes : (N,) tuple
labels for positive end of voxel axes. Dropped axes get a label of
None.
Examples
--------
>>> aff = [[0,1,0,10],[-1,0,0,20],[0,0,1,30],[0,0,0,1]]
>>> aff2axcodes(aff, (('L','R'),('B','F'),('D','U')))
('B', 'R', 'U')
"""
ornt = io_orientation(aff, tol)
return ornt2axcodes(ornt, labels)
|