Skip to content
Snippets Groups Projects
Commit 5f67d58c authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

TEST: Expand resample tests

parent 4e0acb87
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
tests/test_image_resample.py
+ 168
− 133
View file @ 5f67d58c
......@@ -2,7 +2,6 @@
import itertools as it
import os.path as op
import numpy as np
import pytest
......@@ -10,181 +9,217 @@ import pytest
import fsl.data.image as fslimage
import fsl.utils.transform as transform
import fsl.utils.image.resample as resample
from fsl.utils.tempdir import tempdir
from . import make_random_image
def test_resample(seed):
with tempdir() as td:
# Random base image shapes
for i in range(25):
fname = op.join(td, 'test.nii')
shape = np.random.randint(5, 50, 3)
img = fslimage.Image(make_random_image(dims=shape))
# Random base image shapes
for i in range(25):
# bad shape
with pytest.raises(ValueError):
resample.resample(img, (10, 10))
with pytest.raises(ValueError):
resample.resample(img, (10, 10, 10, 10))
shape = np.random.randint(5, 50, 3)
make_random_image(fname, shape)
img = fslimage.Image(fname, mmap=False)
# 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)
# bad shape
with pytest.raises(ValueError):
resample.resample(img, (10, 10))
with pytest.raises(ValueError):
resample.resample(img, (10, 10, 10, 10))
# Random resampled image shapes
for j in range(10):
# 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)
rshape = np.random.randint(5, 50, 3)
resampled, xf = resample.resample(img, rshape, order=0)
# Random resampled image shapes
for j in range(10):
assert tuple(resampled.shape) == tuple(rshape)
rshape = np.random.randint(5, 50, 3)
resampled, xf = resample.resample(img, rshape, order=0)
# 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
img.save('base.nii.gz')
fslimage.Image(resampled, xform=xf,
mmap=False).save('res.nii.gz')
resx, resy, resz = restestcoords.T
resvals = resampled[resx, resy, resz]
assert tuple(resampled.shape) == tuple(rshape)
res2orig = transform.concat(img.worldToVoxMat, xf)
# 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
origtestcoords = transform.transform(restestcoords, res2orig)
resx, resy, resz = restestcoords.T
resvals = resampled[resx, resy, resz]
# 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)
res2orig = transform.concat(img.worldToVoxMat, xf)
origtestcoords = np.array(origtestcoords.round(), dtype=np.int)
origtestcoords = transform.transform(restestcoords, res2orig)
origtestcoords = origtestcoords[~out, :]
restestcoords = restestcoords[ ~out, :]
# 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)
resx, resy, resz = restestcoords.T
origx, origy, origz = origtestcoords.T
origtestcoords = np.array(origtestcoords.round(), dtype=np.int)
origvals = img[:][origx, origy, origz]
resvals = resampled[resx, resy, resz]
origtestcoords = origtestcoords[~out, :]
restestcoords = restestcoords[ ~out, :]
assert np.all(np.isclose(resvals, origvals))
resx, resy, resz = restestcoords.T
origx, origy, origz = origtestcoords.T
origvals = img[:][origx, origy, origz]
resvals = resampled[resx, resy, resz]
def test_resample_4d(seed):
assert np.all(np.isclose(resvals, origvals))
# 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])
del img
img = None
# resample up
resampled = resample.resample(img, (15, 15, 15), slc)[0]
assert tuple(resampled.shape) == (15, 15, 15)
# resample down
resampled = resample.resample(img, (5, 5, 5), slc)[0]
assert tuple(resampled.shape) == (5, 5, 5)
def test_resample_4d(seed):
# resample the entire image
resampled = resample.resample(img, (15, 15, 15, 10), None)[0]
assert tuple(resampled.shape) == (15, 15, 15, 10)
fname = 'test.nii.gz'
resampled = resample.resample(img, (5, 5, 5, 10), None)[0]
assert tuple(resampled.shape) == (5, 5, 5, 10)
with tempdir():
# resample along the fourth dim
resampled = resample.resample(img, (15, 15, 15, 15), None)[0]
assert tuple(resampled.shape) == (15, 15, 15, 15)
make_random_image(fname, (10, 10, 10, 10))
resampled = resample.resample(img, (5, 5, 5, 15), None)[0]
assert tuple(resampled.shape) == (5, 5, 5, 15)
# resample one volume
img = fslimage.Image(fname)
slc = (slice(None), slice(None), slice(None), 3)
resampled = resample.resample(img, img.shape[:3], slc)[0]
assert np.all(resampled == img[..., 3])
# resample up
resampled = resample.resample(img, (15, 15, 15), slc)[0]
assert tuple(resampled.shape) == (15, 15, 15)
def test_resample_origin(seed):
# resample down
resampled = resample.resample(img, (5, 5, 5), slc)[0]
assert tuple(resampled.shape) == (5, 5, 5)
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])
imgb = transform.axisBounds(img.shape, img.voxToWorldMat)
resb = transform.axisBounds(res.shape, res.voxToWorldMat)
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
imgb = np.array(transform.axisBounds(img.shape, img.voxToWorldMat))
resb = np.array(transform.axisBounds(res.shape, res.voxToWorldMat))
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())
# resample the entire image
resampled = resample.resample(img, (15, 15, 15, 10), None)[0]
assert tuple(resampled.shape) == (15, 15, 15, 10)
resampled = resample.resample(img, (5, 5, 5, 10), None)[0]
assert tuple(resampled.shape) == (5, 5, 5, 10)
def test_resampleToPixdims():
# resample along the fourth dim
resampled = resample.resample(img, (15, 15, 15, 15), None)[0]
assert tuple(resampled.shape) == (15, 15, 15, 15)
img = fslimage.Image(make_random_image(dims=(10, 10, 10)))
imglo, imghi = transform.axisBounds(img.shape, img.voxToWorldMat)
oldpix = np.array(img.pixdim, dtype=np.float)
oldshape = np.array(img.shape, dtype=np.float)
resampled = resample.resample(img, (5, 5, 5, 15), None)[0]
assert tuple(resampled.shape) == (5, 5, 5, 15)
for i, origin in it.product(range(25), ('centre', 'corner')):
del img
del resampled
img = None
resampled = None
# random pixdims in the range 0.1 - 5.0
newpix = 0.1 + 4.9 * np.random.random(3)
expshape = np.round(oldshape * (oldpix / newpix))
res = resample.resampleToPixdims(img, newpix, origin=origin)
res = fslimage.Image(res[0], xform=res[1])
reslo, reshi = transform.axisBounds(res.shape, res.voxToWorldMat)
resfov = reshi - reslo
expfov = newpix * res.shape
def test_resample_origin(seed):
with tempdir() as td:
fname = op.join(td, 'test.nii')
make_random_image(fname, (10, 10, 10))
img = fslimage.Image(fname)
# 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])
imgb = transform.axisBounds(img.shape, img.voxToWorldMat)
resb = transform.axisBounds(res.shape, res.voxToWorldMat)
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
imgb = np.array(transform.axisBounds(img.shape, img.voxToWorldMat))
resb = np.array(transform.axisBounds(res.shape, res.voxToWorldMat))
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())
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))
if origin == 'corner':
assert np.all(np.isclose(imglo, reslo))
assert np.all(np.isclose(reshi, reslo + expfov,
rtol=1e-2, atol=1e-2))
def test_resampleToPixdims():
pass
def test_resampleToReference():
pass
def random_v2w():
return transform.compose(
0.25 + 4.75 * np.random.random(3),
-50 + 100 * np.random.random(3),
-np.pi + 2 * np.pi * np.random.random(3))
# 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)
iv2w = random_v2w()
rv2w = random_v2w()
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)
# More specific test - output
# data gets transformed correctly
# into reference space
img = np.zeros((5, 5, 5), dtype=np.float)
img[1, 1, 1] = 1
img = fslimage.Image(img)
refv2w = transform.scaleOffsetXform([1, 1, 1], [-1, -1, -1])
ref = np.zeros((5, 5, 5), dtype=np.float)
ref = fslimage.Image(ref, xform=refv2w)
res = resample.resampleToReference(img, ref, order=0)
exp = np.zeros((5, 5, 5), dtype=np.float)
exp[2, 2, 2] = 1
assert np.all(np.isclose(res[0], exp))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment