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 .