|
| 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 | +""" Image processing functions for: |
| 10 | +
|
| 11 | +* smoothing |
| 12 | +* resampling |
| 13 | +* converting sd to and from FWHM |
| 14 | +
|
| 15 | +Smoothing and resampling routines need scipy |
| 16 | +""" |
| 17 | +from __future__ import print_function, division, absolute_import |
| 18 | + |
| 19 | +import numpy as np |
| 20 | +import numpy.linalg as npl |
| 21 | + |
| 22 | +from .optpkg import optional_package |
| 23 | +spnd, _, _ = optional_package('scipy.ndimage') |
| 24 | + |
| 25 | +from .affines import AffineError, to_matvec, from_matvec, append_diag |
| 26 | +from .spaces import vox2out_vox |
| 27 | +from .nifti1 import Nifti1Image |
| 28 | + |
| 29 | +SIGMA2FWHM = np.sqrt(8 * np.log(2)) |
| 30 | + |
| 31 | + |
| 32 | +def fwhm2sigma(fwhm): |
| 33 | + """ Convert a FWHM value to sigma in a Gaussian kernel. |
| 34 | +
|
| 35 | + Parameters |
| 36 | + ---------- |
| 37 | + fwhm : array-like |
| 38 | + FWHM value or values |
| 39 | +
|
| 40 | + Returns |
| 41 | + ------- |
| 42 | + sigma : array or float |
| 43 | + sigma values corresponding to `fwhm` values |
| 44 | +
|
| 45 | + Examples |
| 46 | + -------- |
| 47 | + >>> sigma = fwhm2sigma(6) |
| 48 | + >>> sigmae = fwhm2sigma([6, 7, 8]) |
| 49 | + >>> sigma == sigmae[0] |
| 50 | + True |
| 51 | + """ |
| 52 | + return np.asarray(fwhm) / SIGMA2FWHM |
| 53 | + |
| 54 | + |
| 55 | +def sigma2fwhm(sigma): |
| 56 | + """ Convert a sigma in a Gaussian kernel to a FWHM value |
| 57 | +
|
| 58 | + Parameters |
| 59 | + ---------- |
| 60 | + sigma : array-like |
| 61 | + sigma value or values |
| 62 | +
|
| 63 | + Returns |
| 64 | + ------- |
| 65 | + fwhm : array or float |
| 66 | + fwhm values corresponding to `sigma` values |
| 67 | +
|
| 68 | + Examples |
| 69 | + -------- |
| 70 | + >>> fwhm = sigma2fwhm(3) |
| 71 | + >>> fwhms = sigma2fwhm([3, 4, 5]) |
| 72 | + >>> fwhm == fwhms[0] |
| 73 | + True |
| 74 | + """ |
| 75 | + return np.asarray(sigma) * SIGMA2FWHM |
| 76 | + |
| 77 | + |
| 78 | +def adapt_affine(affine, n_dim): |
| 79 | + """ Adapt input / output dimensions of spatial `affine` for `n_dims` |
| 80 | +
|
| 81 | + Adapts a spatial (4, 4) affine that is being applied to an image with fewer |
| 82 | + than 3 spatial dimensions, or more than 3 dimensions. If there are more |
| 83 | + than three dimensions, assume an identity transformation for these |
| 84 | + dimensions. |
| 85 | +
|
| 86 | + Parameters |
| 87 | + ---------- |
| 88 | + affine : array-like |
| 89 | + affine transform. Usually shape (4, 4). For what follows ``N, M = |
| 90 | + affine.shape`` |
| 91 | + n_dims : int |
| 92 | + Number of dimensions of underlying array, and therefore number of input |
| 93 | + dimensions for affine. |
| 94 | +
|
| 95 | + Returns |
| 96 | + ------- |
| 97 | + adapted : shape (M, n_dims+1) array |
| 98 | + Affine array adapted to number of input dimensions. Columns of the |
| 99 | + affine corresponding to missing input dimensions have been dropped, |
| 100 | + columns corresponding to extra input dimensions have an extra identity |
| 101 | + column added |
| 102 | + """ |
| 103 | + affine = np.asarray(affine) |
| 104 | + rzs, trans = to_matvec(affine) |
| 105 | + # For missing input dimensions, drop columns in rzs |
| 106 | + rzs = rzs[:, :n_dim] |
| 107 | + adapted = from_matvec(rzs, trans) |
| 108 | + n_extra_columns = n_dim - adapted.shape[1] + 1 |
| 109 | + if n_extra_columns > 0: |
| 110 | + adapted = append_diag(adapted, np.ones((n_extra_columns,))) |
| 111 | + return adapted |
| 112 | + |
| 113 | + |
| 114 | +def resample_from_to(from_img, |
| 115 | + to_vox_map, |
| 116 | + order=3, |
| 117 | + mode='constant', |
| 118 | + cval=0., |
| 119 | + out_class=Nifti1Image): |
| 120 | + """ Resample image `from_img` to mapped voxel space `to_vox_map` |
| 121 | +
|
| 122 | + Resample using N-d spline interpolation. |
| 123 | +
|
| 124 | + Parameters |
| 125 | + ---------- |
| 126 | + from_img : object |
| 127 | + Object having attributes ``dataobj``, ``affine``, ``header``. If |
| 128 | + `out_class` is not None, ``img.__class__`` should be able to construct |
| 129 | + an image from data, affine and header. |
| 130 | + to_vox_map : image object or length 2 sequence |
| 131 | + If object, has attributes ``shape`` giving input voxel shape, and |
| 132 | + ``affine`` giving mapping of input voxels to output space. If length 2 |
| 133 | + sequence, elements are (shape, affine) with same meaning as above. The |
| 134 | + affine is a (4, 4) array-like. |
| 135 | + order : int, optional |
| 136 | + The order of the spline interpolation, default is 3. The order has to |
| 137 | + be in the range 0-5 (see ``scipy.ndimage.affine_transform``) |
| 138 | + mode : str, optional |
| 139 | + Points outside the boundaries of the input are filled according |
| 140 | + to the given mode ('constant', 'nearest', 'reflect' or 'wrap'). |
| 141 | + Default is 'constant' (see ``scipy.ndimage.affine_transform``) |
| 142 | + cval : scalar, optional |
| 143 | + Value used for points outside the boundaries of the input if |
| 144 | + ``mode='constant'``. Default is 0.0 (see |
| 145 | + ``scipy.ndimage.affine_transform``) |
| 146 | + out_class : None or SpatialImage class, optional |
| 147 | + Class of output image. If None, use ``from_img.__class__``. |
| 148 | +
|
| 149 | + Returns |
| 150 | + ------- |
| 151 | + out_img : object |
| 152 | + Image of instance specified by `out_class`, containing data output from |
| 153 | + resampling `from_img` into axes aligned to the output space of |
| 154 | + ``from_img.affine`` |
| 155 | + """ |
| 156 | + try: |
| 157 | + to_shape, to_affine = to_vox_map.shape, to_vox_map.affine |
| 158 | + except AttributeError: |
| 159 | + to_shape, to_affine = to_vox_map |
| 160 | + a_to_affine = adapt_affine(to_affine, len(to_shape)) |
| 161 | + if out_class is None: |
| 162 | + out_class = from_img.__class__ |
| 163 | + from_n_dim = len(from_img.shape) |
| 164 | + if from_n_dim < 3: |
| 165 | + raise AffineError('from_img must be at least 3D') |
| 166 | + a_from_affine = adapt_affine(from_img.affine, from_n_dim) |
| 167 | + to_vox2from_vox = npl.inv(a_from_affine).dot(a_to_affine) |
| 168 | + rzs, trans = to_matvec(to_vox2from_vox) |
| 169 | + data = spnd.affine_transform(from_img.dataobj, |
| 170 | + rzs, |
| 171 | + trans, |
| 172 | + to_shape, |
| 173 | + order=order, |
| 174 | + mode=mode, |
| 175 | + cval=cval) |
| 176 | + return out_class(data, to_affine, from_img.header) |
| 177 | + |
| 178 | + |
| 179 | +def resample_to_output(in_img, |
| 180 | + voxel_sizes=None, |
| 181 | + order=3, |
| 182 | + mode='constant', |
| 183 | + cval=0., |
| 184 | + out_class=Nifti1Image): |
| 185 | + """ Resample image `in_img` to output voxel axes (world space) |
| 186 | +
|
| 187 | + Parameters |
| 188 | + ---------- |
| 189 | + in_img : object |
| 190 | + Object having attributes ``dataobj``, ``affine``, ``header``. If |
| 191 | + `out_class` is not None, ``img.__class__`` should be able to construct |
| 192 | + an image from data, affine and header. |
| 193 | + voxel_sizes : None or sequence |
| 194 | + Gives the diagonal entries of ``out_img.affine` (except the trailing 1 |
| 195 | + for the homogenous coordinates) (``out_img.affine == |
| 196 | + np.diag(voxel_sizes + [1])``). If None, return identity |
| 197 | + `out_img.affine`. If scalar, interpret as vector ``[voxel_sizes] * |
| 198 | + len(in_img.shape)``. |
| 199 | + order : int, optional |
| 200 | + The order of the spline interpolation, default is 3. The order has to |
| 201 | + be in the range 0-5 (see ``scipy.ndimage.affine_transform``). |
| 202 | + mode : str, optional |
| 203 | + Points outside the boundaries of the input are filled according to the |
| 204 | + given mode ('constant', 'nearest', 'reflect' or 'wrap'). Default is |
| 205 | + 'constant' (see ``scipy.ndimage.affine_transform``). |
| 206 | + cval : scalar, optional |
| 207 | + Value used for points outside the boundaries of the input if |
| 208 | + ``mode='constant'``. Default is 0.0 (see |
| 209 | + ``scipy.ndimage.affine_transform``). |
| 210 | + out_class : None or SpatialImage class, optional |
| 211 | + Class of output image. If None, use ``in_img.__class__``. |
| 212 | +
|
| 213 | + Returns |
| 214 | + ------- |
| 215 | + out_img : object |
| 216 | + Image of instance specified by `out_class`, containing data output from |
| 217 | + resampling `in_img` into axes aligned to the output space of |
| 218 | + ``in_img.affine`` |
| 219 | + """ |
| 220 | + if out_class is None: |
| 221 | + out_class = in_img.__class__ |
| 222 | + in_shape = in_img.shape |
| 223 | + n_dim = len(in_shape) |
| 224 | + if voxel_sizes is not None: |
| 225 | + voxel_sizes = np.asarray(voxel_sizes) |
| 226 | + if voxel_sizes.ndim == 0: # Scalar |
| 227 | + voxel_sizes = np.repeat(voxel_sizes, n_dim) |
| 228 | + # Allow 2D images by promoting to 3D. We might want to see what a slice |
| 229 | + # looks like when resampled into world coordinates |
| 230 | + if n_dim < 3: # Expand image to 3D, make voxel sizes match |
| 231 | + new_shape = in_shape + (1,) * (3 - n_dim) |
| 232 | + data = in_img.get_data().reshape(new_shape) # 2D data should be small |
| 233 | + in_img = out_class(data, in_img.affine, in_img.header) |
| 234 | + if voxel_sizes is not None and len(voxel_sizes) == n_dim: |
| 235 | + # Need to pad out voxel sizes to match new image dimensions |
| 236 | + voxel_sizes = tuple(voxel_sizes) + (1,) * (3 - n_dim) |
| 237 | + out_vox_map = vox2out_vox((in_img.shape, in_img.affine), voxel_sizes) |
| 238 | + return resample_from_to(in_img, out_vox_map, order, mode, cval, out_class) |
| 239 | + |
| 240 | + |
| 241 | +def smooth_image(img, |
| 242 | + fwhm, |
| 243 | + mode='nearest', |
| 244 | + cval=0., |
| 245 | + out_class=Nifti1Image): |
| 246 | + """ Smooth image `img` along voxel axes by FWHM `fwhm` millimeters |
| 247 | +
|
| 248 | + Parameters |
| 249 | + ---------- |
| 250 | + img : object |
| 251 | + Object having attributes ``dataobj``, ``affine``, ``header``. If |
| 252 | + `out_class` is not None, ``img.__class__`` should be able to construct |
| 253 | + an image from data, affine and header. |
| 254 | + fwhm : scalar or length 3 sequence |
| 255 | + FWHM *in mm* over which to smooth. The smoothing applies to the voxel |
| 256 | + axes, not to the output axes, but is in millimeters. The function |
| 257 | + adjusts the FWHM to voxels using the voxel sizes calculated from the |
| 258 | + affine. A scalar implies the same smoothing across the spatial |
| 259 | + dimensions of the image, but 0 smoothing over any further dimensions |
| 260 | + such as time. A vector should be the same length as the number of |
| 261 | + image dimensions. |
| 262 | + mode : str, optional |
| 263 | + Points outside the boundaries of the input are filled according |
| 264 | + to the given mode ('constant', 'nearest', 'reflect' or 'wrap'). |
| 265 | + Default is 'nearest'. This is different from the default for |
| 266 | + ``scipy.ndimage.affine_transform``, which is 'constant'. 'nearest' |
| 267 | + might be a better choice when smoothing to the edge of an image where |
| 268 | + there is still strong brain signal, otherwise this signal will get |
| 269 | + blurred towards zero. |
| 270 | + cval : scalar, optional |
| 271 | + Value used for points outside the boundaries of the input if |
| 272 | + ``mode='constant'``. Default is 0.0 (see |
| 273 | + ``scipy.ndimage.affine_transform``). |
| 274 | + out_class : None or SpatialImage class, optional |
| 275 | + Class of output image. If None, use ``img.__class__``. |
| 276 | +
|
| 277 | + Returns |
| 278 | + ------- |
| 279 | + smoothed_img : object |
| 280 | + Image of instance specified by `out_class`, containing data output from |
| 281 | + smoothing `img` data by given FWHM kernel. |
| 282 | + """ |
| 283 | + if out_class is None: |
| 284 | + out_class = img.__class__ |
| 285 | + n_dim = len(img.shape) |
| 286 | + # TODO: make sure time axis is last |
| 287 | + # Pad out fwhm from scalar, adding 0 for fourth etc (time etc) dimensions |
| 288 | + fwhm = np.asarray(fwhm) |
| 289 | + if fwhm.size == 1: |
| 290 | + fwhm_scalar = fwhm |
| 291 | + fwhm = np.zeros((n_dim,)) |
| 292 | + fwhm[:3] = fwhm_scalar |
| 293 | + # Voxel sizes |
| 294 | + RZS = img.affine[:-1, :n_dim] |
| 295 | + vox = np.sqrt(np.sum(RZS ** 2, 0)) |
| 296 | + # Smoothing in terms of voxels |
| 297 | + vox_fwhm = fwhm / vox |
| 298 | + vox_sd = fwhm2sigma(vox_fwhm) |
| 299 | + # Do the smoothing |
| 300 | + sm_data = spnd.gaussian_filter(img.dataobj, |
| 301 | + vox_sd, |
| 302 | + mode=mode, |
| 303 | + cval=cval) |
| 304 | + return out_class(sm_data, img.affine, img.header) |
0 commit comments