Cjg/rigid upgrades

[Cody-G]

This adds a couple of options for Ipopt-based rigid registration and also introduces QuadDIRECT-based rigid registration.

[Cody-G]

Note that I relied on CoordinateTransformations rather than the AffineTransforms package. This is a departure from our other rigid registration code, but it seems like the right move because AffineTransforms seems deprecated. At some point I suppose we can switch everything over to CoordinateTransformations.

[Cody-G]

The test failure is only on nightly.

[timholy]

Really fantastic work, thanks for tackling this!

Merge at will.

[timholy]

Do you know whether it ever takes anything close to this number of evaluations? I guess the termination criteria could make it happen?

[Cody-G]

In my cases so far it usually takes a few hundred evaluations, but I've seen it take 1-2k before. I've never seen it go to 10k, but yes the criteria could allow for that to happen. fval or ftol could help, I'll experiment with those as I continue to run registrations.

[timholy]

Should this really be filling with 0?

[Cody-G]

Oh that does seem risky. I think I'll go with Inf instead of the default NaN?

[timholy]

Also, FYI I just turned on precompilation for QuadDIRECT; I'm surprised its absence didn't cause trouble here (I thought you couldn't use a non-precompiled module inside a precompiled one...). Also, QuadDIRECT changed quite a lot since you first submitted this, in particular the incremental Gaussian elimination only got added three days afterwards. Might be worth experimenting a bit to make sure you're still happy with it, and of course report any problems you discover.

[Cody-G]

Also, QuadDIRECT changed quite a lot since you first submitted this, in particular the incremental Gaussian elimination only got added three days afterwards. Might be worth experimenting a bit to make sure you're still happy with it, and of course report any problems you discover.

The recent updates did change the behavior of the algorithm. In particular this test now fails consistently. It doesn't fail by much. Here's the output from my laptop:

Final minimum (268 evaluations): Box0.0006401215324105769@[0.0031369, 0.0707756, 0.226949, -0.00349066, -0.0013941, 0.00118257]
Test Failed
  Expression: mm < 0.0001
   Evaluated: 0.0006401215324105769 < 0.0001

It seems that the newer version completes fewer iterations before stopping (If you check the last good CI run it took greater than 1k iterations). I may be able to fix this by decreasing ftol but would like to get your opinion first.

[timholy]

Yep, I have the same guess as you. As I'm sure you understand, any change to the algorithm that decreases the iteration count needed for equivalent registration performance is a good thing. However, changes in the algorithm's behavior may also trigger the rtol criterion at different stages of optimization. In this case, I'd bet it's faster to find some local minimum (the improved quasi-Newton optimization is much more reliable in finding local minima), and can then go through several iterations where the minimum is no longer decreasing but it hasn't finished enough "exploration" to realize that there's a better minimum out there. That's why an fvalue criterion, when available, is so nice.

"How do you know you're done?" seems like a much harder problem for global optimization than for local optimization.

[Cody-G]

Interesting, I just tried decreasing rtol all the way to 1e-14 (vs the default of 1e-7) and it still fails that test quite often, even when it does roughly as many iterations as before. Here's an example test result:

Final minimum (1409 evaluations): Box0.0006754915092378871@[-0.00595834, 0.0606882, 0.21179, -0.00349066, -0.00120664, 0.00197094]
Test Failed
  Expression: mm < 0.0001
   Evaluated: 0.0006754915092378871 < 0.0001

Maybe this is a more difficult problem than average "real world" data because the test is registering a set of random pixels that probably have more local minima than structured data, but I'm still a bit troubled by the fact that the test was passing more easily before the recent QuadDIRECT changes.

[Cody-G]

I should also mention that 'rtol' wasn't getting passed properly to the "coarse" registration function, but fixing that doesn't change the results.

[timholy]

The changes fixed some pretty clear bugs, for example other problems that were taking 10x more trials before than they are now. But it seems like it didn't lead to improvement on all problems. I spent some time looking into it, and I noticed some other things that weren't really being handled correctly. If you incorporate timholy/QuadDIRECT.jl#4 I think you'll see it reliably converges.

BTW, I know you got this from the demo/images.jl, but you can probably just use analyze, there's no real need for a second call to analyze! with more iterations. I just did that for further testing purposes. I ran your code with the following diff:

diff --git a/src/qd_rigid.jl b/src/qd_rigid.jl
index 9b69ae1..4d983bd 100644
--- a/src/qd_rigid.jl
+++ b/src/qd_rigid.jl
@@ -67,13 +67,12 @@ function fast_mm(theta, mxshift, fixed, moving, thresh, SD)
     return mm
 end
 
-function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), rtol=1e-7)
+function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), rtol=1e-7, kwargs...)
     f(x) = fast_mm(x, mxshift, fixed, moving, thresh, SD)
     upper = [mxrot...]
     lower = -upper
     splits = ([[-x; 0.0; x] for x in mxrot]...)
-    root0, x0coarse = analyze(f, splits, lower, upper; maxevals=10, minwidth=minwidth_rot)
-    root_coarse = analyze!(deepcopy(root0), f, x0coarse, splits, lower, upper; maxevals=10000, rtol=rtol, minwidth=minwidth_rot, print_interval=100)
+    root_coarse, x0coarse = analyze(f, splits, lower, upper; maxevals=10^4, minwidth=minwidth_rot, rtol=rtol, print_interval=100, kwargs...)
     box_coarse = minimum(root_coarse)
     tfmcoarse0 = rot(position(box_coarse, x0coarse), moving)
     best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm = tfmcoarse0)
@@ -81,7 +80,7 @@ function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh
     return tfmcoarse, mm
 end
 
-function qd_rigid_fine(fixed, moving, mxrot, minwidth, SD; initial_tfm = -1, thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), rtol=1e-7)
+function qd_rigid_fine(fixed, moving, mxrot, minwidth, SD; initial_tfm = -1, thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), rtol=1e-7, kwargs...)
     f2(x) = slow_mm(x, fixed, moving, thresh, SD; initial_tfm = initial_tfm)
     upper_shft = fill(1, ndims(fixed))
     upper_rot = mxrot
@@ -91,8 +90,7 @@ function qd_rigid_fine(fixed, moving, mxrot, minwidth, SD; initial_tfm = -1, thr
     minwidth_shfts = fill(0.01, ndims(fixed))
     minwidth_rots = ndims(fixed) == 2 ? [0.0001;] : fill(0.0001, 3) #assume 2 or 3 dimensional input
     minwidth = vcat(minwidth_shfts, minwidth_rots)
-    root0, x0 = analyze(f2, splits, lower, upper; maxevals=10, minwidth=minwidth)
-    root = analyze!(deepcopy(root0), f2, x0, splits, lower, upper; maxevals=10000, rtol=rtol, minwidth=minwidth, print_interval=100)
+    root, x0 = analyze(f2, splits, lower, upper; maxevals=10^4, minwidth=minwidth, print_interval=100, rtol=rtol, kwargs...)
     box = minimum(root)
     params = position(box, x0)
     tfmfine = tfmx(position(box, x0), moving)
@@ -109,12 +107,12 @@ Use `SD` if your axes are not uniformly sampled, for example `SD = diagm(voxelsp
 is a vector encoding the spacing along all axes of the image. `thresh` enforces a certain amount of sum-of-squared-intensity overlap between
 the two images; with non-zero `thresh`, it is not permissible to "align" the images by shifting one entirely out of the way of the other.
 """
-function qd_rigid(fixed, moving, mxshift, mxrot, minwidth_rot, SD=eye(ndims(fixed)); thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), rtol=1e-7)
+function qd_rigid(fixed, moving, mxshift, mxrot, minwidth_rot, SD=eye(ndims(fixed)); thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), rtol=1e-7, kwargs...)
     mxrot = [mxrot...]
     print("Running coarse step\n")
-    tfm_coarse, mm_coarse = qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh = thresh, rtol=1e-5)
+    tfm_coarse, mm_coarse = qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh = thresh, rtol=1e-5, kwargs...)
     print("Running fine step\n")
-    tfm_fine, mm_fine = qd_rigid_fine(fixed, moving, mxrot./10, minwidth_rot, SD; initial_tfm = tfm_coarse, thresh = thresh, rtol=rtol)
+    tfm_fine, mm_fine = qd_rigid_fine(fixed, moving, mxrot./10, minwidth_rot, SD; initial_tfm = tfm_coarse, thresh = thresh, rtol=rtol, kwargs...)
     final_tfm = tfm_fine ∘ tfm_coarse
     return final_tfm, mm_fine
 end

Aside from consolidating the analyze calls, this allows one to pass through additional keyword arguments (e.g., I can use fvalue when I'm confident about the ground truth).

[timholy]

The main problem with the new version is, as usual, it often does a lot of additional iterations after one might (retrospectively) say it converged.

[Cody-G]

timholy/QuadDIRECT.jl#4 did the trick, thanks! I've updated the PR to pass kwargs as you suggested. Agreed that it's still easy to do more iterations than necessary, but at least it's finding a low minimum -- lower than even before the changes! I'll merge after Travis runs again.