test_image_resample.py 10 KB
Newer Older
1
2
3
4
5
6
7
8
#!/usr/bin/env python


import itertools as it
import numpy     as np

import pytest

9
10
import scipy.ndimage       as ndimage

11
import fsl.data.image           as     fslimage
Paul McCarthy's avatar
Paul McCarthy committed
12
import fsl.transform.affine     as     affine
13
14
15
16
import fsl.utils.image.resample as     resample

from . import make_random_image

17
def random_affine():
Paul McCarthy's avatar
Paul McCarthy committed
18
    return affine.compose(
19
20
21
22
23
        0.25   + 4.75      * np.random.random(3),
        -50    + 100       * np.random.random(3),
        -np.pi + 2 * np.pi * np.random.random(3))


24
25
26

def test_resample(seed):

Paul McCarthy's avatar
Paul McCarthy committed
27
28
    # Random base image shapes
    for i in range(25):
29

Paul McCarthy's avatar
Paul McCarthy committed
30
31
        shape = np.random.randint(5, 50, 3)
        img = fslimage.Image(make_random_image(dims=shape))
32

Paul McCarthy's avatar
Paul McCarthy committed
33
34
35
36
37
        # bad shape
        with pytest.raises(ValueError):
            resample.resample(img, (10, 10))
        with pytest.raises(ValueError):
            resample.resample(img, (10, 10, 10, 10))
38

Paul McCarthy's avatar
Paul McCarthy committed
39
40
41
42
        # resampling to the same shape should be a no-op
        samei, samex = resample.resample(img, shape)
        assert np.all(samei == img[:])
        assert np.all(samex == img.voxToWorldMat)
43

Paul McCarthy's avatar
Paul McCarthy committed
44
45
        # Random resampled image shapes
        for j in range(10):
46

Paul McCarthy's avatar
Paul McCarthy committed
47
48
            rshape        = np.random.randint(5, 50, 3)
            resampled, xf = resample.resample(img, rshape, order=0)
49

Paul McCarthy's avatar
Paul McCarthy committed
50
            assert tuple(resampled.shape) == tuple(rshape)
51

Paul McCarthy's avatar
Paul McCarthy committed
52
53
54
55
56
57
58
59
            # We used nearest neighbour interp, so the
            # values in the resampled image should match
            # corresponding values in the original. Let's
            # check some whynot.
            restestcoords = np.array([
                np.random.randint(0, rshape[0], 100),
                np.random.randint(0, rshape[1], 100),
                np.random.randint(0, rshape[2], 100)]).T
60

Paul McCarthy's avatar
Paul McCarthy committed
61
62
            resx,  resy,  resz  = restestcoords.T
            resvals  = resampled[resx, resy, resz]
63

Paul McCarthy's avatar
Paul McCarthy committed
64
            res2orig = affine.concat(img.worldToVoxMat, xf)
65

Paul McCarthy's avatar
Paul McCarthy committed
66
            origtestcoords = affine.transform(restestcoords, res2orig)
67

Paul McCarthy's avatar
Paul McCarthy committed
68
69
70
71
72
73
74
75
76
            # remove any coordinates which are out of
            # bounds in the original image space, or
            # are right on a voxel boundary (where the
            # nn interp could have gone either way), or
            # have value == 0 in the resampled space.
            out = ((origtestcoords < 0)            |
                   (origtestcoords >= shape - 0.5) |
                   (np.isclose(np.modf(origtestcoords)[0], 0.5)))
            out = np.any(out, axis=1) | (resvals == 0)
77

78
            origtestcoords = np.array(origtestcoords.round(), dtype=int)
79

Paul McCarthy's avatar
Paul McCarthy committed
80
81
            origtestcoords = origtestcoords[~out, :]
            restestcoords  = restestcoords[ ~out, :]
82

Paul McCarthy's avatar
Paul McCarthy committed
83
84
            resx,  resy,  resz  = restestcoords.T
            origx, origy, origz = origtestcoords.T
85

Paul McCarthy's avatar
Paul McCarthy committed
86
87
            origvals = img[:][origx, origy, origz]
            resvals  = resampled[resx, resy, resz]
88

Paul McCarthy's avatar
Paul McCarthy committed
89
            assert np.all(np.isclose(resvals, origvals))
90
91


Paul McCarthy's avatar
Paul McCarthy committed
92
def test_resample_4d(seed):
93

Paul McCarthy's avatar
Paul McCarthy committed
94
95
96
97
98
    # resample one volume
    img = fslimage.Image(make_random_image(dims=(10, 10, 10, 10)))
    slc = (slice(None), slice(None), slice(None), 3)
    resampled = resample.resample(img, img.shape[:3], slc)[0]
    assert np.all(resampled == img[..., 3])
99

Paul McCarthy's avatar
Paul McCarthy committed
100
101
102
    # resample up
    resampled = resample.resample(img, (15, 15, 15), slc)[0]
    assert tuple(resampled.shape) == (15, 15, 15)
103

Paul McCarthy's avatar
Paul McCarthy committed
104
105
106
    # resample down
    resampled = resample.resample(img, (5, 5, 5), slc)[0]
    assert tuple(resampled.shape) == (5, 5, 5)
107

Paul McCarthy's avatar
Paul McCarthy committed
108
109
110
    # resample the entire image
    resampled = resample.resample(img, (15, 15, 15, 10), None)[0]
    assert tuple(resampled.shape) == (15, 15, 15, 10)
111

Paul McCarthy's avatar
Paul McCarthy committed
112
113
    resampled = resample.resample(img, (5, 5, 5, 10), None)[0]
    assert tuple(resampled.shape) == (5, 5, 5, 10)
114

Paul McCarthy's avatar
Paul McCarthy committed
115
116
117
    # resample along the fourth dim
    resampled = resample.resample(img, (15, 15, 15, 15), None)[0]
    assert tuple(resampled.shape) == (15, 15, 15, 15)
118

Paul McCarthy's avatar
Paul McCarthy committed
119
120
    resampled = resample.resample(img, (5, 5, 5, 15), None)[0]
    assert tuple(resampled.shape) == (5, 5, 5, 15)
121
122


Paul McCarthy's avatar
Paul McCarthy committed
123
def test_resample_origin(seed):
124

Paul McCarthy's avatar
Paul McCarthy committed
125
126
127
128
129
130
131
132
    img = fslimage.Image(make_random_image(dims=(10, 10, 10)))

    # with origin='corner', image
    # bounding boxes should match
    for i in range(25):
        shape = np.random.randint(5, 50, 3)
        res = resample.resample(img, shape, origin='corner')
        res = fslimage.Image(res[0], xform=res[1])
Paul McCarthy's avatar
Paul McCarthy committed
133
134
        imgb = affine.axisBounds(img.shape, img.voxToWorldMat)
        resb = affine.axisBounds(res.shape, res.voxToWorldMat)
Paul McCarthy's avatar
Paul McCarthy committed
135
136
137
138
139
140
141
142
143
144
        assert np.all(np.isclose(imgb, resb, rtol=1e-5, atol=1e-5))

    # with origin='centre' image
    # bounding boxes should be offset
    # by (size_resampled - size_orig) / 2
    for i in range(25):
        shape = np.random.randint(5, 50, 3)
        res = resample.resample(img, shape, origin='centre')
        res = fslimage.Image(res[0], xform=res[1])
        off = (np.array(img.shape) / np.array(res.shape) - 1) / 2
Paul McCarthy's avatar
Paul McCarthy committed
145
146
        imgb = np.array(affine.axisBounds(img.shape, img.voxToWorldMat))
        resb = np.array(affine.axisBounds(res.shape, res.voxToWorldMat))
Paul McCarthy's avatar
Paul McCarthy committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
        assert np.all(np.isclose(imgb, resb + off, rtol=1e-5, atol=1e-5))

    # with origin='corner', using
    # linear interp, when we down-
    # sample an image to a shape
    # that divides evenly into the
    # original shape, a downsampled
    # voxel should equal the average
    # of the original voxels inside
    # it.
    res = resample.resample(img, (5, 5, 5), smooth=False, origin='corner')[0]
    for x, y, z in it.product(range(5), range(5), range(5)):
        block = img[x * 2: x * 2 + 2,
                    y * 2: y * 2 + 2,
                    z * 2: z * 2 + 2]
        assert np.isclose(res[x, y, z], block.mean())
163
164


Paul McCarthy's avatar
Paul McCarthy committed
165
def test_resampleToPixdims():
166

Paul McCarthy's avatar
Paul McCarthy committed
167
    img          = fslimage.Image(make_random_image(dims=(10, 10, 10)))
Paul McCarthy's avatar
Paul McCarthy committed
168
    imglo, imghi = affine.axisBounds(img.shape, img.voxToWorldMat)
169
170
    oldpix       = np.array(img.pixdim, dtype=float)
    oldshape     = np.array(img.shape,  dtype=float)
171

Paul McCarthy's avatar
Paul McCarthy committed
172
    for i, origin in it.product(range(25), ('centre', 'corner')):
173

Paul McCarthy's avatar
Paul McCarthy committed
174
175
176
        # random pixdims in the range 0.1 - 5.0
        newpix   = 0.1 + 4.9 * np.random.random(3)
        expshape = np.round(oldshape * (oldpix / newpix))
177

Paul McCarthy's avatar
Paul McCarthy committed
178
179
        res = resample.resampleToPixdims(img, newpix, origin=origin)
        res = fslimage.Image(res[0], xform=res[1])
Paul McCarthy's avatar
Paul McCarthy committed
180
        reslo, reshi = affine.axisBounds(res.shape, res.voxToWorldMat)
Paul McCarthy's avatar
Paul McCarthy committed
181
182
        resfov       = reshi - reslo
        expfov       = newpix * res.shape
183

Paul McCarthy's avatar
Paul McCarthy committed
184
185
186
        assert np.all(np.isclose(res.shape,  expshape))
        assert np.all(np.isclose(res.pixdim, newpix))
        assert np.all(np.isclose(resfov, expfov, rtol=1e-2, atol=1e-2))
187

Paul McCarthy's avatar
Paul McCarthy committed
188
189
190
191
        if origin == 'corner':
            assert np.all(np.isclose(imglo, reslo))
            assert np.all(np.isclose(reshi, reslo + expfov,
                                     rtol=1e-2, atol=1e-2))
192
193
194



195
def test_resampleToReference1():
Paul McCarthy's avatar
Paul McCarthy committed
196
197
198
199
200
201
202

    # Basic test - output has same
    # dimensions/space as reference
    for i in range(25):

        ishape = np.random.randint(5, 50, 3)
        rshape = np.random.randint(5, 50, 3)
203
204
        iv2w   = random_affine()
        rv2w   = random_affine()
Paul McCarthy's avatar
Paul McCarthy committed
205
206
207
208
209
210
211
        img    = fslimage.Image(make_random_image(dims=ishape, xform=iv2w))
        ref    = fslimage.Image(make_random_image(dims=rshape, xform=rv2w))
        res    = resample.resampleToReference(img, ref)
        res    = fslimage.Image(res[0], header=ref.header)

        assert res.sameSpace(ref)

212
213
214

def test_resampleToReference2():

Paul McCarthy's avatar
Paul McCarthy committed
215
216
217
    # More specific test - output
    # data gets transformed correctly
    # into reference space
218
    img          = np.zeros((5, 5, 5), dtype=float)
Paul McCarthy's avatar
Paul McCarthy committed
219
220
221
    img[1, 1, 1] = 1
    img          = fslimage.Image(img)

Paul McCarthy's avatar
Paul McCarthy committed
222
    refv2w = affine.scaleOffsetXform([1, 1, 1], [-1, -1, -1])
223
    ref    = np.zeros((5, 5, 5), dtype=float)
Paul McCarthy's avatar
Paul McCarthy committed
224
225
226
    ref    = fslimage.Image(ref, xform=refv2w)
    res    = resample.resampleToReference(img, ref, order=0)

227
    exp          = np.zeros((5, 5, 5), dtype=float)
Paul McCarthy's avatar
Paul McCarthy committed
228
229
230
    exp[2, 2, 2] = 1

    assert np.all(np.isclose(res[0], exp))
231
232
233
234
235
236
237


def test_resampleToReference3():

    # Test resampling image to ref
    # with mismatched dimensions
    imgdata = np.random.randint(0, 65536, (5, 5, 5))
Paul McCarthy's avatar
Paul McCarthy committed
238
    img     = fslimage.Image(imgdata, xform=affine.scaleOffsetXform(
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
        (2, 2, 2), (0.5, 0.5, 0.5)))

    # reference/expected data when
    # resampled to ref (using nn interp).
    # Same as image, upsampled by a
    # factor of 2
    refdata = np.repeat(np.repeat(np.repeat(imgdata, 2, 0), 2, 1), 2, 2)
    refdata = np.array([refdata] * 8).transpose((1, 2, 3, 0))
    ref     = fslimage.Image(refdata)

    # We should be able to use a 4D reference
    resampled, xform = resample.resampleToReference(img, ref, order=0, mode='nearest')
    assert np.all(resampled == ref.data[..., 0])

    # If resampling a 4D image with a 3D reference,
    # the fourth dimension should be passed through
    resampled, xform = resample.resampleToReference(ref, img, order=0, mode='nearest')
    exp = np.array([imgdata] * 8).transpose((1, 2, 3, 0))
    assert np.all(resampled == exp)

    # When resampling 4D to 4D, only the
    # first 3 dimensions should be resampled
    imgdata = np.array([imgdata] * 15).transpose((1, 2, 3, 0))
    img     = fslimage.Image(imgdata, xform=img.voxToWorldMat)
    exp     = np.array([refdata[..., 0]] * 15).transpose((1, 2, 3, 0))
    resampled, xform = resample.resampleToReference(img, ref, order=0, mode='nearest')
    assert np.all(resampled == exp)
266
267
268
269
270
271
272


def test_resampleToReference4():

    # the image and ref are out of
    # alignment, but this affine
    # will bring them into alignment
Paul McCarthy's avatar
Paul McCarthy committed
273
    img2ref = affine.scaleOffsetXform([2, 2, 2], [10, 10, 10])
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290

    imgdata = np.random.randint(0, 65536, (5, 5, 5))
    refdata = np.zeros((5, 5, 5))
    img     = fslimage.Image(imgdata)
    ref     = fslimage.Image(refdata, xform=img2ref)

    # Without the affine, the image
    # will be out of the FOV of the
    # reference
    resampled, xform = resample.resampleToReference(img, ref)
    assert np.all(resampled == 0)

    # But applying the affine will
    # cause them to overlap
    # perfectly in world coordinates
    resampled, xform = resample.resampleToReference(img, ref, matrix=img2ref)
    assert np.all(resampled == imgdata)