Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix diffeomorphic registration test failures #489

Merged
merged 12 commits into from Dec 8, 2014
+486 −301
Merged

Fix diffeomorphic registration test failures #489

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
f0321dd
NF: write __config__.py with compile variables
matthew-brett Dec 4, 2014
f80a478
TEST: compare against different energy profile with SSE2
omarocegueda Nov 15, 2014
271fc95
Export build flags to __config__.py
omarocegueda Nov 29, 2014
c096da5
TEST: use exported build flags to select oracle energy profile
omarocegueda Nov 29, 2014
6c5726f
TEST: skip regression tests when SSE2 is not supported
omarocegueda Nov 29, 2014
74afdff
BF: apparently, 'has_key(...)' was removed from dict in python3
omarocegueda Nov 29, 2014
a10b3f9
Generate the config file earlier
omarocegueda Nov 30, 2014
0438941
HACK: generate __config__ module as a cython extension to make sure i…
omarocegueda Dec 1, 2014
4c76179
Remove __config__.pyx, will be auto-generated
omarocegueda Dec 1, 2014
54543eb
Using simpler __setup__.py, suggested by Matthew
omarocegueda Dec 5, 2014
8eac874
Remove DipyConfig and build_flags classes
omarocegueda Dec 6, 2014
24dc623
PEP8: remove extra empty lines and set maximum line length to 80 char…
omarocegueda Dec 8, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Expand Up @@ -24,3 +24,4 @@ dist/
dipy.egg-info/
nibabel*egg/
pyx-stamps
__config__.py
3 changes: 3 additions & 0 deletions dipy/align/imwarp.py
View file
Open in desktop
Expand Up @@ -1452,8 +1452,11 @@ def _iterate(self):
#set zero displacements at the boundary
fw_step[0, ...] = 0
fw_step[:, 0, ...] = 0
fw_step[-1, ...] = 0
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need to set to zero all the boundaries. I missed the other side!

fw_step[:, -1, ...] = 0
if(self.dim == 3):
fw_step[:, :, 0, ...] = 0
fw_step[:, :, -1, ...] = 0

#Normalize the forward step
nrm = np.sqrt(np.sum((fw_step/current_disp_spacing)**2, -1)).max()
Expand Down
6 changes: 5 additions & 1 deletion dipy/align/tests/test_crosscorr.py
View file
Open in desktop
Expand Up @@ -11,6 +11,8 @@ def test_cc_factors_2d():
"""
a = np.array(range(20*20), dtype=floating).reshape(20,20)
b = np.array(range(20*20)[::-1], dtype=floating).reshape(20,20)
a /= a.max()
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is to keep the values in a manageable range (not too large) so we don't get different results in the first decimals when the processor is using 64 bits or 80 bits

b /= b.max()
for radius in [0, 1, 3, 6]:
factors = np.asarray(cc.precompute_cc_factors_2d(a,b,radius))
expected = np.asarray(cc.precompute_cc_factors_2d_test(a,b,radius))
Expand All @@ -25,10 +27,12 @@ def test_cc_factors_3d():
"""
a = np.array(range(20*20*20), dtype=floating).reshape(20,20,20)
b = np.array(range(20*20*20)[::-1], dtype=floating).reshape(20,20,20)
a /= a.max()
b /= b.max()
for radius in [0, 1, 3, 6]:
factors = np.asarray(cc.precompute_cc_factors_3d(a,b,radius))
expected = np.asarray(cc.precompute_cc_factors_3d_test(a,b,radius))
assert_array_almost_equal(factors, expected)
assert_array_almost_equal(factors, expected, decimal=5)


def test_compute_cc_steps_2d():
Expand Down
437 changes: 240 additions & 197 deletions dipy/align/tests/test_imwarp.py
View file
Open in desktop

Large diffs are not rendered by default.

154 changes: 111 additions & 43 deletions dipy/align/tests/test_vector_fields.py
View file
Open in desktop
Expand Up @@ -240,24 +240,40 @@ def test_interpolate_scalar_2d():
target_shape = (sz, sz)
image = np.ndarray(target_shape, dtype=floating)
image[...] = np.random.randint(0, 10, np.size(image)).reshape(target_shape)
#Select some coordinates to interpolate at

extended_image = np.zeros((sz+2, sz+2), dtype=floating)
extended_image[1:sz+1, 1:sz+1] = image[...]

#Select some coordinates inside the image to interpolate at
nsamples = 200
locations = np.random.ranf(2 * nsamples).reshape((nsamples,2)) * (sz+2) - 1.0
extended_locations = locations + 1.0 # shift coordinates one voxel

#Call the implementation under test
interp, inside = vfu.interpolate_scalar_2d(image, locations)

#Call the reference implementation
expected = map_coordinates(image, locations.transpose(), order=1)
expected = map_coordinates(extended_image, extended_locations.transpose(), order=1)

assert_array_almost_equal(expected, interp)

#Test the 'inside' flag
for i in range(nsamples):
if (locations[i, 0]<0 or locations[i, 0]>(sz-1)) or (locations[i, 1]<0 or locations[i, 1]>(sz-1)):
assert_equal(inside[i], 0)
else:
assert_equal(inside[i], 1)
#Test interpolation stability along the boundary
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added this test to verify that small differences in the interpolation coordinates near the boundary produce small differences in the interpolation result (i.e. the new interpolation is more stable along the boundary)

epsilon = 5e-8
for k in range(2):
for offset in [0, sz-1]:
delta = ((np.random.ranf(nsamples) * 2) -1) * epsilon
locations[:, k] = delta + offset
locations[:, (k+1)%2] = np.random.ranf(nsamples) * (sz-1)
interp, inside = vfu.interpolate_scalar_2d(image, locations)

locations[:, k] = offset
expected = map_coordinates(image, locations.transpose(), order=1)
assert_array_almost_equal(expected, interp)
if offset == 0:
expected_flag = np.array(delta>=0, dtype = np.int32)
else:
expected_flag = np.array(delta<=0, dtype = np.int32)
assert_array_almost_equal(expected_flag, inside)


def test_interpolate_scalar_nn_2d():
Expand Down Expand Up @@ -320,26 +336,42 @@ def test_interpolate_scalar_3d():
target_shape = (sz, sz, sz)
image = np.ndarray(target_shape, dtype=floating)
image[...] = np.random.randint(0, 10, np.size(image)).reshape(target_shape)
#Select some coordinates to interpolate at
nsamples = 200

extended_image = np.zeros((sz+2, sz+2, sz+2), dtype=floating)
extended_image[1:sz+1, 1:sz+1, 1:sz+1] = image[...]

#Select some coordinates inside the image to interpolate at
nsamples = 800
locations = np.random.ranf(3 * nsamples).reshape((nsamples,3)) * (sz+2) - 1.0
extended_locations = locations + 1.0 # shift coordinates one voxel

#Call the implementation under test
interp, inside = vfu.interpolate_scalar_3d(image, locations)

#Call the reference implementation
expected = map_coordinates(image, locations.transpose(), order=1)
expected = map_coordinates(extended_image, extended_locations.transpose(), order=1)

assert_array_almost_equal(expected, interp)

#Test the 'inside' flag
for i in range(nsamples):
expected_inside = 1
for axis in range(3):
if (locations[i, axis]<0 or locations[i, axis]>(sz-1)):
expected_inside = 0
break
assert_equal(inside[i], expected_inside)
#Test interpolation stability along the boundary
epsilon = 5e-8
for k in range(3):
for offset in [0, sz-1]:
delta = ((np.random.ranf(nsamples) * 2) -1) * epsilon
locations[:, k] = delta + offset
locations[:, (k+1)%3] = np.random.ranf(nsamples) * (sz-1)
locations[:, (k+2)%3] = np.random.ranf(nsamples) * (sz-1)
interp, inside = vfu.interpolate_scalar_3d(image, locations)

locations[:, k] = offset
expected = map_coordinates(image, locations.transpose(), order=1)
assert_array_almost_equal(expected, interp)

if offset == 0:
expected_flag = np.array(delta>=0, dtype = np.int32)
else:
expected_flag = np.array(delta<=0, dtype = np.int32)
assert_array_almost_equal(expected_flag, inside)


def test_interpolate_vector_3d():
Expand All @@ -348,28 +380,44 @@ def test_interpolate_vector_3d():
target_shape = (sz, sz, sz)
field = np.ndarray(target_shape+(3,), dtype=floating)
field[...] = np.random.randint(0, 10, np.size(field)).reshape(target_shape+(3,))

extended_field = np.zeros((sz+2, sz+2, sz+2, 3), dtype=floating)
extended_field[1:sz+1, 1:sz+1, 1:sz+1] = field
#Select some coordinates to interpolate at
nsamples = 200
nsamples = 800
locations = np.random.ranf(3 * nsamples).reshape((nsamples,3)) * (sz+2) - 1.0
extended_locations = locations + 1

#Call the implementation under test
interp, inside = vfu.interpolate_vector_3d(field, locations)

#Call the reference implementation
expected = np.zeros_like(interp)
for i in range(3):
expected[...,i] = map_coordinates(field[...,i], locations.transpose(), order=1)
expected[...,i] = map_coordinates(extended_field[...,i], extended_locations.transpose(), order=1)

assert_array_almost_equal(expected, interp)

#Test the 'inside' flag
for i in range(nsamples):
expected_inside = 1
for axis in range(3):
if (locations[i, axis]<0 or locations[i, axis]>(sz-1)):
expected_inside = 0
break
assert_equal(inside[i], expected_inside)
#Test interpolation stability along the boundary
epsilon = 5e-8
for k in range(3):
for offset in [0, sz-1]:
delta = ((np.random.ranf(nsamples) * 2) -1) * epsilon
locations[:, k] = delta + offset
locations[:, (k+1)%3] = np.random.ranf(nsamples) * (sz-1)
locations[:, (k+2)%3] = np.random.ranf(nsamples) * (sz-1)
interp, inside = vfu.interpolate_vector_3d(field, locations)

locations[:, k] = offset
for i in range(3):
expected[...,i] = map_coordinates(field[...,i], locations.transpose(), order=1)
assert_array_almost_equal(expected, interp)

if offset == 0:
expected_flag = np.array(delta>=0, dtype = np.int32)
else:
expected_flag = np.array(delta<=0, dtype = np.int32)
assert_array_almost_equal(expected_flag, inside)


def test_interpolate_vector_2d():
Expand All @@ -378,28 +426,43 @@ def test_interpolate_vector_2d():
target_shape = (sz, sz)
field = np.ndarray(target_shape+(2,), dtype=floating)
field[...] = np.random.randint(0, 10, np.size(field)).reshape(target_shape+(2,))
extended_field = np.zeros((sz+2, sz+2, 2), dtype=floating)
extended_field[1:sz+1, 1:sz+1] = field
#Select some coordinates to interpolate at
nsamples = 200
locations = np.random.ranf(2 * nsamples).reshape((nsamples,2)) * (sz+2) - 1.0
extended_locations = locations + 1

#Call the implementation under test
interp, inside = vfu.interpolate_vector_2d(field, locations)

#Call the reference implementation
expected = np.zeros_like(interp)
for i in range(2):
expected[...,i] = map_coordinates(field[...,i], locations.transpose(), order=1)
expected[...,i] = map_coordinates(extended_field[...,i], extended_locations.transpose(), order=1)

assert_array_almost_equal(expected, interp)

#Test the 'inside' flag
for i in range(nsamples):
expected_inside = 1
for axis in range(2):
if (locations[i, axis]<0 or locations[i, axis]>(sz-1)):
expected_inside = 0
break
assert_equal(inside[i], expected_inside)
#Test interpolation stability along the boundary
epsilon = 5e-8
for k in range(2):
for offset in [0, sz-1]:
delta = ((np.random.ranf(nsamples) * 2) -1) * epsilon
locations[:, k] = delta + offset
locations[:, (k+1)%2] = np.random.ranf(nsamples) * (sz-1)
interp, inside = vfu.interpolate_vector_2d(field, locations)

locations[:, k] = offset
for i in range(2):
expected[...,i] = map_coordinates(field[...,i], locations.transpose(), order=1)
assert_array_almost_equal(expected, interp)

if offset == 0:
expected_flag = np.array(delta>=0, dtype = np.int32)
else:
expected_flag = np.array(delta<=0, dtype = np.int32)
assert_array_almost_equal(expected_flag, inside)



def test_warping_2d():
Expand Down Expand Up @@ -466,12 +529,14 @@ def test_warping_2d():
# Reorient the displacement field according to its grid-to-space transform
dcopy = np.copy(d)
vfu.reorient_vector_field_2d(dcopy, field_affine)
extended_dcopy = np.zeros((nr+2, nc+2, 2), dtype=floating)
extended_dcopy[1:nr+1, 1:nc+1, :] = dcopy

# Compute the warping coordinates (see warp_2d documentation)
Y = np.apply_along_axis(A.dot, 0, X)[0:2,...]
Z = np.zeros_like(X)
Z[0,...] = map_coordinates(dcopy[...,0], Y, order=1)
Z[1,...] = map_coordinates(dcopy[...,1], Y, order=1)
Z[0,...] = map_coordinates(extended_dcopy[...,0], Y + 1, order=1)
Z[1,...] = map_coordinates(extended_dcopy[...,1], Y + 1, order=1)
Z[2,...] = 0
Z = np.apply_along_axis(C.dot, 0, Z)[0:2,...]
T = np.apply_along_axis(B.dot, 0, X)[0:2,...]
Expand Down Expand Up @@ -557,12 +622,15 @@ def test_warping_3d():
dcopy = np.copy(d)
vfu.reorient_vector_field_3d(dcopy, field_affine)

extended_dcopy = np.zeros((ns+2, nr+2, nc+2, 3), dtype=floating)
extended_dcopy[1:ns+1, 1:nr+1, 1:nc+1, :] = dcopy

# Compute the warping coordinates (see warp_2d documentation)
Y = np.apply_along_axis(A.dot, 0, X)[0:3,...]
Z = np.zeros_like(X)
Z[0,...] = map_coordinates(dcopy[...,0], Y, order=1)
Z[1,...] = map_coordinates(dcopy[...,1], Y, order=1)
Z[2,...] = map_coordinates(dcopy[...,2], Y, order=1)
Z[0,...] = map_coordinates(extended_dcopy[...,0], Y + 1, order=1)
Z[1,...] = map_coordinates(extended_dcopy[...,1], Y + 1, order=1)
Z[2,...] = map_coordinates(extended_dcopy[...,2], Y + 1, order=1)
Z[3,...] = 0
Z = np.apply_along_axis(C.dot, 0, Z)[0:3,...]
T = np.apply_along_axis(B.dot, 0, X)[0:3,...]
Expand Down