add 3d rigid registration test

[Cody-G]

Added a 3 dimensional test of rigid registration. The test currently fails, highlighting a problem with optimization of 3d transformations. I'm trying to track down the problem now, but wanted to get this out in case it takes me a while. The Ipopt optimization of the affine transformation seems to be stuck evaluating the same point over and over.

[timholy]

I think we'll need to add ImageMagick as a test dependency---this didn't even get as far as the problem you were wanting to illustrate.

My guess: it's been snookered by the fact that the function isn't continuous at the edges, when pixels drop in and out. You could put a cap on the number of iterations, see the Ipopt options here.

Originally I used Optim's Nelder-Mead, which doesn't rely as much on continuity (it doesn't use derivatives). But I found that Nelder-Mead was very unreliable in terms of finding a good solution.

[Cody-G]

Strange, I didn't add anything new--just reused the cameraman image that you already used for the 2d test. I think I was prompted interactively to install ImageMagick when I first ran the test. I suppose it would be better to add it to the REQUIRE in the test folder?

I think the failure is not even as interesting as being snookered. The only point that ever gets evaluated is the initial point, so it returns the identity transformation. As far as I can tell the gradient always evaluates to zero, so I'm wondering if this is a problem with ForwardDiff. It would be consistent with your idea that the objective is discontinuous at pixel boundaries, except that for the 3d test I have done nothing but copy the cameraman image into three dimensions, and rotate it about the same axis. So from an optimization standpoint it should be the same problem.

I have been reading to make sure I understand the machinery of MathProgBase and the Ipopt interface to see if there could be a problem there.

I should have mentioned that you also get this error when registering a 3d image:

ERROR: LoadError: MethodError: `_rotation3` has no method matching _rotation3(::Array{ForwardDiff.GradientNumber{6,Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}},1}, ::ForwardDiff.GradientNumber{6,Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}})

I fixed that by just sending to rotation3 the Float64 representation of the gradient number, i.e. by calling real() on it. But maybe that's not valid? Strangely, 2d registration never has this issue because p2rigid is always passed an array of Float64s

[timholy]

Strange, I didn't add anything new--just reused the cameraman image that you already used for the 2d test. I think I was prompted interactively to install ImageMagick when I first ran the test.

Images got updated to use FileIO, so this is new behavior.

I suppose it would be better to add it to the REQUIRE in the test folder?

Yes, that was what I meant.

I fixed that by just sending to rotation3 the Float64 representation of the gradient number, i.e. by calling real() on it.

That's almost certainly why you are getting a gradient of 0. Try for the 2d image and you'll see that it works there. You probably need to loosen the signature of _rotation3, see timholy/AffineTransforms.jl@60b4f36 for inspiration.

[timholy]

I can't dig into your bug now, but for your reference:
https://en.wikipedia.org/wiki/Dual_number#Differentiation
http://julialang.org/blog/2015/10/auto-diff-in-julia/

Definitely one of the coolest things I've learned by coming to Julia.

[timholy]

Strangely, 2d registration never has this issue because p2rigid is always passed an array of Float64s

Are you sure? I have a hard time believing that given these definitions.

[Cody-G]

You probably need to loosen the signature of _rotation3, see timholy/AffineTransforms.jl@60b4f36 for inspiration.

I tried relaxing the signature as you suggested (after having removed the call to real()) and I get this error:

ERROR: LoadError: InexactError()
 in trunc at ./float.jl:362
 [inlined code] from cartesian.jl:33
 in _transform! at /home/cody/git/juliapackages_new/v0.4/AffineTransforms/src/tformedarrays.jl:32
 [inlined code] from /home/cody/git/juliapackages_new/v0.4/AffineTransforms/src/tformedarrays.jl:95
 in transform! at /home/cody/git/juliapackages_new/v0.4/AffineTransforms/src/tformedarrays.jl:127
 in transform at /home/cody/git/juliapackages_new/v0.4/AffineTransforms/src/tformedarrays.jl:68
 in transform at /home/cody/git/juliapackages_new/v0.4/AffineTransforms/src/tformedarrays.jl:72
 in call at /home/cody/git/juliapackages_new/v0.4/BlockRegistration/src/RegisterOptimize.jl:173
 in _calc_gradient at /home/cody/git/juliapackages_new/v0.4/ForwardDiff/src/api/gradient.jl:86
 in g at /home/cody/git/juliapackages_new/v0.4/ForwardDiff/src/api/gradient.jl:54
 [inlined code] from show.jl:127
 in eval_grad_f at /home/cody/git/juliapackages_new/v0.4/BlockRegistration/src/RegisterOptimize.jl:213
 in eval_grad_f_cb at /home/cody/git/juliapackages_new/v0.4/Ipopt/src/IpoptSolverInterface.jl:293
 in eval_grad_f_wrapper at /home/cody/git/juliapackages_new/v0.4/Ipopt/src/Ipopt.jl:113
while loading /home/cody/git/juliapackages_new/v0.4/BlockRegistrationScheduler/test/rigid.jl, in expression starting on line 38

The line number is wrong, but the error is inside a cartesian loop inside a generated function. I'm not sure how to diagnose it. Actually I find macros hard to debug in general, so if you have any tips I would appreciate them! Maybe the computation gets corrupted by an integer being mixed with a floating point operation? I looked at the types passed to _transform! in both the 2d and 3d case, and they seem to be the same except for the extra dimension.

[Cody-G]

Definitely one of the coolest things I've learned by coming to Julia.

Yes, so cool!!

[Cody-G]

Are you sure? I have a hard time believing that given these definitions.

I don't immediately see what you mean. The p2rigid function doesn't specialize on arguments. It looks like when the objective function (which is a RigidOpt) is evaluated by the solver, the solver calls your RigidOpt.call(), which is also not specialized on x. The p2rigid must get passed whatever type the solver is passing it under the hood. Are you saying that it's hard to believe that the solver would be inconsistent in what it passes?

I was just noticing that if you show p inside of p2rigid, sometimes it's an Array{GradientNumber} and other times it's an Array{Float64} But now that I have looked at it more, I don't think it's an important observation.

[timholy]

The line number is wrong, but the error is inside a cartesian loop inside a generated function. I'm not sure how to diagnose it. Actually I find macros hard to debug in general, so if you have any tips I would appreciate them!

Yes, macros still make it a pain. Better than before, but the line numbering still has issues.

Putting (as you must have done) @show src dest at the top of the script is a good start. In cases where I'm desperate, I do this:

• create variables in the REPL with the same name and type as those in the function
• copy/paste the body of the function into

macroexpand(quote 
    # function body goes here
end)

and evaluate it.

• copy/paste the output to a julia file (editing as needed to make it work).

You now have a version of the function customized for your types that doesn't involve any macros.

I have a private "package" called MacroExpandJL that combines the last two steps and fixes some of the things you'd otherwise have to edit. Let me know if you want it.

The p2rigid must get passed whatever type the solver is passing it under the hood.

The solver will always pass a vector of Float64s, but when the solver asks for the gradient, my code passes a vector of GradientNumbers. (Specifically, the generated g function from ForwardDiff.gradient does that, see https://github.com/JuliaDiff/ForwardDiff.jl/blob/2f505104d30a752847033a9f3f2720f6afad9cfc/src/api/gradient.jl#L77.)

Are you saying that it's hard to believe that the solver would be inconsistent in what it passes?

I mean that I find it hard to believe your statement that p2rigid is always passed an array of Float64s, and that it never gets passed GradientNumbers. I'm basically sure it does; that's why I had to make that adjustment to AffineTransforms.

[Cody-G]

Thanks, I will try that tip for working with macros from now on. And maybe I'll also ask for you package sometime. I've nearly located the problem now. Try this:

using ForwardDiff
#I took the value of axis directly from the code when a transform is generated about the [0;0;0] axis
axis = ForwardDiff.GradientNumber{6,Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}[ForwardDiff.GradientNumber{6,Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}(-0.0,ForwardDiff.Partials{Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}((1.0,0.0,0.0,0.0,0.0,0.0))),ForwardDiff.GradientNumber{6,Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}(-0.0,ForwardDiff.Partials{Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}((0.0,1.0,0.0,0.0,0.0,0.0))),ForwardDiff.GradientNumber{6,Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}(-0.0,ForwardDiff.Partials{Float64,Tuple{Float64,Float64,Float64,Float64,Float64,Float64}}((0.0,0.0,1.0,0.0,0.0,0.0)))]
norm(axis) #should give you NaNs, even though norm([0;0;0]) does not

When you write out the norm calculation instead of using the library function, you see that square root is being applied to a negative 0.0, resulting in a NaN. I suppose negative 0.0 is useful for automatic differentiation, but it's causing some problem here?

[Cody-G]

I think I understand now. ForwardDiff is correct when it calculates NaN for the gradient at the initial point. That's because the derivative of the square root function is not defined at zero. All of this is because the angle of rotation is encoded in the norm of the axis vector, which means the gradient is undefined at the starting point. I confirmed this by providing a tiny rotation for the initial guess. The optimization runs correctly in that case. So what do you think is the best way to fix this? I think if we start with even a very small rotation in the wrong direction then we could get a much worse result. Would it be better to pass to Ipopt all elements of the rotation matrix?

[timholy]

Good progress!

Would it be better to pass to Ipopt all elements of the rotation matrix?

That would indeed solve this particular problem. It introduces another, however; in 3d, a general matrix has 9 parameters, but a rotation matrix has only 3. Assuming you really want it to be a rotation and not a general affine transformation, you'd need to set up the constraints (see http://mathprogbasejl.readthedocs.org/en/latest/nlp.html) to enforce the rotation. Probably an easier approach would be to pass 4 parameters to https://github.com/timholy/AffineTransforms.jl/blob/eeee407b3df35bfc94caab123a60af71ec414993/src/transform.jl#L199 and add a single constraint that sumabs2(uaxis) = 1. This is equivalent to the quaternion representation of rotations.

Note there is nothing stopping us from using a general affine transformation, and that mathematically and code-wise this problem is actually easier than imposing rigidity constraints. However, if the actual underlying deformation really is a rotation, allowing the flexibility of an affine transformation may make it harder to find the best optimum. A greedy algorithm will start going downhill in whatever direction makes it initially most productive, even if that starts by adding skew or other factors that ultimately steer it towards the wrong solution. If the optimization problem were convex this wouldn't matter (except to make the process slower), but image registration definitely isn't, and the existence of multiple local minima is a genuine problem in practice. Enforcing the constraint of it being a rigid transformation definitely reduces the complexity of the parameter space. Given that the fish's brain probably doesn't distort much when the fish moves, I suspect the rigidity assumption makes sense. (If there is a little bit of skew, it's likely that the best approach would still be to start by fitting a rigid transformation, and then doing a second optimization with a general affine transformation.)

However, an approach that requires no new code might be to start with 1e-5*randn(3). Let's say your longest axis is 2000 pixels; from the center, that's a radius of 1000 pixels. The angle that causes a single-pixel shift at the edge of the image is asin(1/1000) ≈ 0.001, so anything smaller than that will have a pretty negligible effect on the alignment.

[timholy]

Would also be great to get this merged.