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

Enhance SLR with applying series of transformations and fix optimize bug for parameter missing in old scipy versions #470

Closed
wants to merge 13 commits into from
Closed

Enhance SLR with applying series of transformations and fix optimize bug for parameter missing in old scipy versions #470

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
8a0453c
ENH: allowing cascades of linear transforms e.g. first rigid then aff…
Garyfallidis Nov 19, 2014
67ff60f
TEST: checking different testing ideas for the cascade
Garyfallidis Nov 19, 2014
f114b16
Removed some empty lines
Garyfallidis Nov 20, 2014
ccb86ea
TEST: Completed test for applying a series of transformations
Garyfallidis Nov 20, 2014
0a05eab
BF: added default options for optimizer for older scipy versions
Garyfallidis Nov 20, 2014
d991c23
TEST: covering more versions of options
Garyfallidis Nov 20, 2014
a741106
Adding exception catch for 0.11 to 0.12 parameter missing in L-BFGS-B
Garyfallidis Nov 20, 2014
c8e36f1
Null commit to wake Travis
Garyfallidis Nov 20, 2014
6486b77
Using maxfun when maxiter is not available
Garyfallidis Nov 21, 2014
54a431b
BF: allowing for None check for iterations
Garyfallidis Nov 21, 2014
4b9be1b
Changing verbosity to identify where Travis stops
Garyfallidis Nov 21, 2014
507990b
Reduced to standard verbosity to identify where Travis stops
Garyfallidis Nov 21, 2014
fa68569
Test was too strict. Removing iterations asserts.
Garyfallidis Nov 21, 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
36 changes: 26 additions & 10 deletions dipy/align/streamlinear.py
View file
Open in desktop
Expand Up @@ -262,7 +262,7 @@ def __init__(self, metric=None, x0="rigid", method='L-BFGS-B',
self.options = options
self.evolution = evolution

def optimize(self, static, moving):
def optimize(self, static, moving, mat=None):
""" Find the minimum of the provided metric.

Parameters
Expand All @@ -272,6 +272,11 @@ def optimize(self, static, moving):
Reference or fixed set of streamlines.
moving : streamlines
Moving set of streamlines.
mat : array
Transformation (4, 4) matrix to start the registration. ``mat``
is applied to moving. Default value None which means that initial
transformation will be generated by shiftin the centers of moving
Copy link
Contributor

Choose a reason for hiding this comment

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

typo: shifting

Copy link
Contributor

Choose a reason for hiding this comment

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

This seems like a strange default. Isn't it better to have the default be np.eye(4)? That is, the default is not to pre-apply any transformation?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

No eye(4) would be a very bad strategy. The default idea in registration is that you first find the center of mass of the two objects that you want to register and then you move them to the origin. In that way you make sure that your rotations during the registration will make sense. Is that clear?

and static sets of streamlines to the origin.

Returns
-------
Expand All @@ -291,8 +296,19 @@ def optimize(self, static, moving):
if not np.all(np.array(list(map(len, moving))) == static[0].shape[0]):
raise ValueError('Static and moving streamlines ' + msg)

static_centered, static_shift = center_streamlines(static)
moving_centered, moving_shift = center_streamlines(moving)
if mat is None:
static_centered, static_shift = center_streamlines(static)
moving_centered, moving_shift = center_streamlines(moving)
static_mat = compose_matrix44([static_shift[0], static_shift[1],
static_shift[2], 0, 0, 0])

moving_mat = compose_matrix44([-moving_shift[0], -moving_shift[1],
-moving_shift[2], 0, 0, 0])
else:
static_centered = static
moving_centered = transform_streamlines(moving, mat)
static_mat = np.eye(4)
moving_mat = mat

self.metric.setup(static_centered, moving_centered)

Expand Down Expand Up @@ -324,12 +340,6 @@ def optimize(self, static, moving):

opt_mat = compose_matrix44(opt.xopt)

static_mat = compose_matrix44([static_shift[0], static_shift[1],
static_shift[2], 0, 0, 0])

moving_mat = compose_matrix44([-moving_shift[0], -moving_shift[1],
-moving_shift[2], 0, 0, 0])

mat = compose_transformations(moving_mat, opt_mat, static_mat)

mat_history = []
Expand Down Expand Up @@ -371,7 +381,13 @@ def _set_x0(self, x0):
if x0 not in [6, 7, 12]:
msg = 'Only 6, 7 and 12 are accepted as integers'
raise ValueError(msg)
return x0
else:
if x0 == 6:
return np.zeros(6)
if x0 == 7:
return np.array([0, 0, 0, 0, 0, 0, 1.])
if x0 == 12:
return np.array([0, 0, 0, 0, 0, 0, 1., 1., 1., 0, 0, 0])

raise ValueError('Wrong input')

Expand Down
56 changes: 41 additions & 15 deletions dipy/align/tests/test_streamlinear.py
View file
Open in desktop
@@ -1,5 +1,6 @@
import numpy as np
from numpy.testing import (run_module_suite,
assert_,
assert_equal,
assert_almost_equal,
assert_array_equal,
Expand Down Expand Up @@ -112,22 +113,14 @@ def test_rigid_partial_real_bundles():
moving_center, shift2 = center_streamlines(moving)

print(shift2)
#mat = compose_matrix44([0, 0, 0, 40, 0, 0])
#print(mat)
mat = compose_matrix(translate=np.array([0, 0, 0.]),
angles=np.deg2rad([40, 0, 0.]))
#print(mat)
moved = transform_streamlines(moving_center, mat)

srr = StreamlineLinearRegistration()

srm = srr.optimize(static_center, moved)
print(srm.fopt)

# bmd = BundleMinDistanceMetric()
# bmd.setup(static_center, moving)
# print(bmd.distance(np.array([0, 0, 0, 0., 0., 0])))

print(srm.iterations)
print(srm.funcs)

Expand All @@ -154,7 +147,6 @@ def test_rigid_partial_real_bundles():
vol2[i, j, k] = 1

overlap = np.sum(np.logical_and(vol, vol2)) / float(np.sum(vol2))
#print(overlap)

assert_equal(overlap * 100 > 40, True)

Expand All @@ -169,17 +161,12 @@ def test_stream_rigid():
moving = transform_streamlines(moving, mat)

srr = StreamlineLinearRegistration()

sr_params = srr.optimize(static, moving)

moved = transform_streamlines(moving, sr_params.matrix)

srr = StreamlineLinearRegistration(verbose=True)

srm = srr.optimize(static, moving)

moved2 = transform_streamlines(moving, srm.matrix)

moved3 = srm.transform(moving)

assert_array_equal(moved[0], moved2[0])
Expand Down Expand Up @@ -435,9 +422,11 @@ def test_x0_input():
StreamlineLinearRegistration(x0=x0)

for x0 in [8, 20, "Whatever", np.random.rand(20), np.random.rand(20, 3)]:

assert_raises(ValueError, StreamlineLinearRegistration, x0=x0)

x0 = np.random.rand(4, 3)
assert_raises(ValueError, StreamlineLinearRegistration, x0=x0)


def test_compose_decompose_matrix44():

Expand All @@ -453,6 +442,43 @@ def test_compose_decompose_matrix44():
assert_raises(ValueError, decompose_matrix44, mat, 20)


def test_cascade_of_optimizations():

cingulum_bundles = two_cingulum_bundles()

cb1 = cingulum_bundles[0]
cb1 = set_number_of_points(cb1, 20)

test_x0 = np.array([10, 4, 3, 0, 20, 10, 1., 1.3, 0.9, 0.1, 0.2, -0.2])

cb2 = transform_streamlines(cingulum_bundles[0],
compose_matrix44(test_x0))
cb2 = set_number_of_points(cb2, 20)

# first rigid
slr = StreamlineLinearRegistration(x0=6)
slm = slr.optimize(cb1, cb2)

# then similarity
slr2 = StreamlineLinearRegistration(x0=7)
slm2 = slr2.optimize(cb1, cb2, slm.matrix)

# then affine
slr3 = StreamlineLinearRegistration(x0=12, options={'maxiter': 400})
slm3 = slr3.optimize(cb1, cb2, slm2.matrix)

# all affine params at once
slr4 = StreamlineLinearRegistration(x0=12, options={'maxiter': 400})
slm4 = slr4.optimize(cb1, cb2)

assert_array_almost_equal(slm4.matrix, slm3.matrix, 2)
assert_(slm2.fopt < slm.fopt)
assert_(slm3.fopt < slm2.fopt)
assert_almost_equal(slm4.fopt, slm3.fopt, 4)
assert_(slm3.iterations < slm4.iterations)
assert_(slm3.funcs < slm4.funcs)


if __name__ == '__main__':

run_module_suite()
54 changes: 44 additions & 10 deletions dipy/core/optimize.py
View file
Open in desktop
Expand Up @@ -13,11 +13,8 @@
SCIPY_LESS_0_12 = LooseVersion(scipy.__version__) < '0.12'

if not SCIPY_LESS_0_12:

from scipy.optimize import minimize

else:

from scipy.optimize import fmin_l_bfgs_b, fmin_powell


Expand Down Expand Up @@ -128,22 +125,48 @@ def __init__(self, fun, x0, args=(), method='L-BFGS-B', jac=None,
self.tmp_files = []
self._evol_kx = None

_eps = np.finfo(float).eps

if SCIPY_LESS_0_12:

if method == 'L-BFGS-B':
default_options = {'maxcor': 10, 'ftol': 1e-7, 'gtol': 1e-5,
'eps': 1e-8, 'maxiter': 15000}

if jac is None:
approx_grad = True
else:
approx_grad = False

out = fmin_l_bfgs_b(fun, x0, args,
approx_grad=approx_grad,
bounds=bounds,
m=options['maxcor'],
factr=options['ftol'] / np.finfo(float).eps,
pgtol=options['gtol'],
epsilon=options['eps'])
if options is None:
options = default_options

if options is not None:
for key in options:
default_options[key] = options[key]
options = default_options

try:
out = fmin_l_bfgs_b(fun, x0, args,
approx_grad=approx_grad,
bounds=bounds,
m=options['maxcor'],
factr=options['ftol']/_eps,
pgtol=options['gtol'],
epsilon=options['eps'],
maxiter=options['maxiter'])
except TypeError:

msg = 'In current version of Scipy `maxiter` parameter is'
msg += ' not available for L-BFGS-B.'
print(msg)
out = fmin_l_bfgs_b(fun, x0, args,
approx_grad=approx_grad,
bounds=bounds,
m=options['maxcor'],
factr=options['ftol']/_eps,
pgtol=options['gtol'],
epsilon=options['eps'])

res = {'x': out[0], 'fun': out[1], 'nfev': out[2]['funcalls']}
try:
Expand All @@ -153,6 +176,17 @@ def __init__(self, fun, x0, args=(), method='L-BFGS-B', jac=None,

elif method == 'Powell':

default_options = {'xtol': 0.0001, 'ftol': 0.0001,
'maxiter': None}

if options is None:
options = default_options

if options is not None:
for key in options:
default_options[key] = options[key]
options = default_options

out = fmin_powell(fun, x0, args,
xtol=options['xtol'],
ftol=options['ftol'],
Expand Down
49 changes: 44 additions & 5 deletions dipy/core/tests/test_optimize.py
View file
Open in desktop
Expand Up @@ -24,17 +24,14 @@ def func2(x):
opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]), method='Powell')

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))

assert_almost_equal(opt.fopt, 0)

opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]), method='L-BFGS-B',
options={'maxcor': 10, 'ftol': 1e-7,
'gtol': 1e-5, 'eps': 1e-8})

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))

assert_almost_equal(opt.fopt, 0)

assert_equal(opt.evolution, None)

opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]), method='L-BFGS-B',
Expand All @@ -43,7 +40,6 @@ def func2(x):
evolution=False)

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))

assert_almost_equal(opt.fopt, 0)

opt.print_summary()
Expand Down Expand Up @@ -84,7 +80,6 @@ def func2(x):
'gtol': 1e-5, 'eps': 1e-8})

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))

assert_almost_equal(opt.fopt, 0)

print(opt.nit)
Expand All @@ -102,6 +97,50 @@ def func2(x):

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))

opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]),
method='L-BFGS-B',
options={'maxcor': 10, 'eps': 1e-8})

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
assert_almost_equal(opt.fopt, 0)

print(opt.nit)
print(opt.fopt)
print(opt.nfev)

opt = Optimizer(fun=func, x0=np.array([1., 1., 1.]),
method='L-BFGS-B',
options=None)

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0]))
assert_almost_equal(opt.fopt, 0)

print(opt.nit)
print(opt.fopt)
print(opt.nfev)

opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
method='Powell',
options={'maxiter': 1e6},
evolution=True)

print(opt.nit)
print(opt.fopt)
print(opt.nfev)

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))

opt = Optimizer(fun=func2, x0=np.array([1., 1., 1., 5.]),
method='Powell',
options={'maxiter': 1e6},
evolution=True)

print(opt.nit)
print(opt.fopt)
print(opt.nfev)

assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.]))


if __name__ == '__main__':

Expand Down