resample.py 10.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
#!/usr/bin/env python
#
# resample.py - The resample functino
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module defines the :func:`resample` function, which can be used
to resample an :class:`.Image` object to a different resolution.

10
11
12
The :func:`resampleToPixdims` and :func:`resampleToReference` functions
are convenience wrappers around :func:`resample`.

13
The :func:`applySmoothing` function is a sub-function of :func:`resample`.
14
15
16
"""


17
18
import numpy                as np
import scipy.ndimage        as ndimage
19

20
import fsl.transform.affine as affine
21
import fsl.utils.deprecated as deprecated
22
23


24
25
26
27
28
29
30
31
32
def resampleToPixdims(image, newPixdims, **kwargs):
    """Resample ``image`` so that it has the specified voxel dimensions.

    This is a wrapper around :func:`resample` - refer to its documenttion
    for details on the other arguments and the return values.

    :arg image:   :class:`.Image` to resample
    :arg pixdims: New voxel dimensions to resample ``image`` to.
    """
Paul McCarthy's avatar
Paul McCarthy committed
33
34
35
36
    newPixdims = np.array(newPixdims)
    oldShape   = np.array(image.shape)
    oldPixdims = np.array(image.pixdim)
    newShape   = oldShape * (oldPixdims / newPixdims)
37
38
39
    return resample(image, newShape, **kwargs)


40
def resampleToReference(image, reference, matrix=None, **kwargs):
41
42
43
44
45
    """Resample ``image`` into the space of the ``reference``.

    This is a wrapper around :func:`resample` - refer to its documenttion
    for details on the other arguments and the return values.

46
47
48
    When resampling to a reference image, resampling will only be applied
    along the spatial (first three) dimensions.

49
    :arg image:     :class:`.Image` to resample
50
    :arg matrix:    Optional world-to-world affine alignment matrix
51
52
53
54
    :arg reference: :class:`.Nifti` defining the space to resample ``image``
                    into
    """

55
56
57
    oldShape = list(image.shape)
    newShape = list(reference.shape[:3])

58
59
60
    if matrix is None:
        matrix = np.eye(4)

61
62
63
64
65
66
67
    # If the input image is >3D, pad the
    # new shape so that we only resample
    # along the first 3 dimensions.
    if len(newShape) < len(oldShape):
        newShape = newShape + oldShape[len(newShape):]

    # Align the two images together
68
69
70
    # via their vox-to-world affines,
    # and the world-to-world affine
    # if provided
71
72
73
    matrix = affine.concat(image.worldToVoxMat,
                           affine.invert(matrix),
                           reference.voxToWorldMat)
74
75
76
77
78
79
80
81
82
83
84
85
86

    # If the input image is >3D, we
    # have to adjust the resampling
    # matrix to take into account the
    # additional dimensions (see scipy.
    # ndimage.affine_transform)
    if len(newShape) > 3:
        rotmat  = matrix[:3, :3]
        offsets = matrix[:3,  3]
        matrix  = np.eye(len(newShape) + 1)
        matrix[:3, :3] = rotmat
        matrix[:3, -1] = offsets

87
    kwargs['mode']     = kwargs.get('mode', 'constant')
88
89
90
    kwargs['newShape'] = newShape
    kwargs['matrix']   = matrix

91
92
93
94
95
96
    data = resample(image, **kwargs)[0]

    # The image is now in the same space
    # as the reference, so it inherits
    # ref's voxel-to-world affine
    return data, reference.voxToWorldMat
97
98


99
100
101
102
103
104
def resample(image,
             newShape,
             sliceobj=None,
             dtype=None,
             order=1,
             smooth=True,
Paul McCarthy's avatar
Paul McCarthy committed
105
             origin=None,
106
             matrix=None,
Paul McCarthy's avatar
Paul McCarthy committed
107
             mode=None,
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
             cval=0):
    """Returns a copy of the data in the ``image``, resampled to the specified
    ``newShape``.

    The space that the image is resampled into can be defined in one of the
    following ways, in decreasing order of precedence:

      1. If a ``matrix`` is provided, it is applied to the voxel coordinates
         when retrieving values from the ``image``

      2. Otherwise the image is simply scaled according to the ratio calculated
         by ``image.shape / newShape``. In this case the ``origin`` argument
         may be used to adjust the alignemnt of the original and resampled
         voxel grids.

    See the ``scipy.ndimage.affine_transform`` function for more details,
    particularly on the ``order``, ``matrix``, ``mode`` and
    ``cval`` arguments.

127
128
129
130
131
132
    .. note:: If a custom resampling ``matrix`` is specified, the adjusted
              voxel-to-world afffine cannot be calculated by this function,
              so ``None`` will be returned instead.

    :arg image:    :class:`.Image` object to resample

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
    :arg newShape: Desired shape. May containg floating point values, in which
                   case the resampled image will have shape
                   ``round(newShape)``, but the voxel sizes will have scales
                   ``self.shape / newShape`` (unless ``matrix`` is specified).

    :arg sliceobj: Slice into this ``Image``. If ``None``, the whole
                   image is resampled, and it is assumed that it has the
                   same number of dimensions as  ``newShape``. A
                   :exc:`ValueError` is raised if this is not the case.

    :arg dtype:    ``numpy`` data type of the resampled data. If ``None``,
                   the :meth:`dtype` of this ``Image`` is used.

    :arg order:    Spline interpolation order, passed through to the
                   ``scipy.ndimage.affine_transform`` function - ``0``
                   corresponds to nearest neighbour interpolation, ``1``
                   (the default) to linear interpolation, and ``3`` to
                   cubic interpolation.

    :arg smooth:   If ``True`` (the default), the data is smoothed before
                   being resampled, but only along axes which are being
                   down-sampled (i.e. where ``newShape[i] < self.shape[i]``).

    :arg origin:   ``'centre'`` (the default) or ``'corner'``. ``'centre'``
                   resamples the image such that the centre of the corner
                   voxels of this image and the resampled data are
                   aligned. ``'corner'`` resamples the image such that
                   the corner of the corner voxels are aligned (and
                   therefore the voxel grids are aligned).
                   Ignored if ``offset`` or ``matrix`` is specified.

    :arg matrix:   Arbitrary affine transformation matrix to apply to the
                   voxel coordinates of ``image`` when resampling.

    :arg mode:     How to handle regions which are outside of the image FOV.
                   Defaults to `''nearest'``.

    :arg cval:     Constant value to use when ``mode='constant'``.

    :returns: A tuple containing:

               - A ``numpy`` array of shape ``newShape``, containing
                 an interpolated copy of the data in this ``Image``.

               - A ``numpy`` array of shape ``(4, 4)``, containing the
                 adjusted voxel-to-world transformation for the spatial
179
180
                 dimensions of the resampled data, or ``None`` if a ``matrix``
                 was provided.
181
182
183
184
    """

    if sliceobj is None:     sliceobj = slice(None)
    if dtype    is None:     dtype    = image.dtype
Paul McCarthy's avatar
Paul McCarthy committed
185
    if origin   is None:     origin   = 'centre'
Paul McCarthy's avatar
Paul McCarthy committed
186
    if mode     is None:     mode     = 'nearest'
187
188
    if origin   == 'center': origin   = 'centre'

189
190
    ownMatrix = matrix is None

191
192
193
194
195
196
197
198
199
200
201
202
203
204
    if origin not in ('centre', 'corner'):
        raise ValueError('Invalid value for origin: {}'.format(origin))

    data = np.array(image[sliceobj], dtype=dtype, copy=False)

    if len(data.shape) != len(newShape):
        raise ValueError('Data dimensions do not match new shape: '
                         'len({}) != len({})'.format(data.shape, newShape))

    # If matrix not provided, calculate
    # a scaling/offset matrix from the
    # old/new shape ratio and the origin
    # setting.
    if matrix is None:
205
        matrix = affine.rescale(data.shape, newShape, origin)
206

207
    # identity matrix? the image
208
    # doesn't need to be resampled
209
    if np.all(np.isclose(matrix, np.eye(len(newShape) + 1))):
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
        return data, image.voxToWorldMat

    newShape = np.array(np.round(newShape), dtype=np.int)

    # Apply smoothing if requested,
    # and if not using nn interp
    if order > 0 and smooth:
        data = applySmoothing(data, matrix, newShape)

    # Do the resample thing
    data = ndimage.affine_transform(data,
                                    matrix,
                                    output_shape=newShape,
                                    order=order,
                                    mode=mode,
                                    cval=cval)

    # Construct an affine transform which
    # puts the resampled image into the
    # same world coordinate system as this
230
231
232
    # image. We may be working with >3D data,
    # so here we discard the non-spatial
    # parts of the resampling matrix
233
    if matrix.shape != (4, 4):
234
235
236
237
238
239
        rotmat         = matrix[:3, :3]
        offsets        = matrix[:3, -1]
        matrix         = np.eye(4)
        matrix[:3, :3] = rotmat
        matrix[:3, -1] = offsets

240
241
242
243
244
245
246
247
248
    # We can only adjust the v2w affine if
    # the input space and resampling space
    # are aligned (e.g. if we're just
    # resampling to different dimensions).
    # We can't assume this when a custom
    # matrix is specified, so we just give
    # up and return None.
    if ownMatrix: matrix = affine.concat(image.voxToWorldMat, matrix)
    else:         matrix = None
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270

    return data, matrix


def applySmoothing(data, matrix, newShape):
    """Called by the :func:`resample` function.

    If interpolating and smoothing, we apply a gaussian filter along axes with
    a resampling ratio greater than 1.1. We do this so that interpolation has
    an effect when down-sampling to a resolution where the voxel centres are
    aligned (as otherwise any interpolation regime will be equivalent to
    nearest neighbour). This more-or-less mimics the behaviour of FLIRT.

    See the ``scipy.ndimage.gaussian_filter`` function for more details.

    :arg data:     Data to be smoothed.
    :arg matrix:   Affine matrix to be used during resampling. The voxel
                   scaling factors are extracted from this.
    :arg newShape: Shape the data is to be resampled into.
    :returns:      A smoothed copy of ``data``.
    """

271
    ratio = affine.decompose(matrix[:3, :3])[0]
272
273
274
275
276
277
278
279
280
281
282
283
284
285

    if len(newShape) > 3:
        ratio = np.concatenate((
            ratio,
            [float(o) / float(s)
             for o, s in zip(data.shape[3:], newShape[3:])]))

    sigma                = np.array(ratio)
    sigma[ratio <  1.1]  = 0
    sigma[ratio >= 1.1] *= 0.425

    return ndimage.gaussian_filter(data, sigma)


286
287
@deprecated.deprecated('2.9.0', '3.0.0',
                       'Use fsl.transform.affine.rescale instead')
288
def calculateMatrix(oldShape, newShape, origin):
289
290
291
    """Deprecated - use :func:`.affine.rescale` instead. """
    xform = affine.rescale(oldShape, newShape, origin)
    if np.all(np.isclose(xform, np.eye(len(oldShape) + 1))):
292
        return None
293
    return xform[:-1, :]