Add qd_translate and qd_affine methods for QuadDIRECT registration #81
Conversation
…ion. Also add and reorganize tests.
There was a problem hiding this comment.
This is looking really good. Thanks for enhancing the functionality so dramatically.
Two thoughts:
- I suspect @ChantalJuntao has another pull request coming. If it's close, we might want to merge that first (since this is a pretty dramatic reorganization and that's tough for newcomers to cope with from the standpoint of git etc)
- I suspect it might be worthwhile poking a bit more at why the affine transformation test fails. I am sure some of the limitations are in QuadDIRECT (and that won't be fixed until I have its replacement ready), but some may also be in subtle issues here. Before merging I might want to play with this a bit and see if I can gain more insight.
| error() | ||
| end | ||
| #################### Utilities ########################## | ||
| function inds_intersect(img1, img2) |
There was a problem hiding this comment.
You can also write this more compactly as intersect.(indices(A), indices(B)). Might not need a function for that? But I'm not against having one.
| end | ||
| newmov = warp(moving, tfm) | ||
| inds = inds_intersect(newmov, fixed) | ||
| return newmov[inds...], view(fixed, inds...) |
There was a problem hiding this comment.
I noticed this before, but didn't ask. Now I'm curious, why use newmov[inds...] but view(fixed, inds...)?
There was a problem hiding this comment.
I wanted to use views for both, but I remember I was getting an error. Now I can't remember what the error was (it's been a while) but I'll try changing it back and see.
There was a problem hiding this comment.
I found the error again and opened issue #83
Sorry for not reporting it sooner, I think I was in a hurry at the time I encountered it and then forgot.
There was a problem hiding this comment.
Thanks for opening the issue and finding a minimal reproducer. I'll try to get that one figured out myself.
| error() | ||
| function warp_and_intersect(moving, fixed, tfm) | ||
| if tfm == IdentityTransformation() | ||
| return moving, fixed |
There was a problem hiding this comment.
For full consistency here, before returning we should check that indices(moving) == indices(fixed); if that's not true, we should still intersect and subset.
Or, if having "extra pixels" in one of the images might be useful, then we shouldn't trim at all. Either way, it might be best to make only the warp part conditional on tfm and keep the various branches more parallel.
There was a problem hiding this comment.
I take your first point, and as for the extra pixels I can imagine them being useful sometimes. I can think of at least one case where that would require other changes: we're using mismatch0 in some places and that doesn't work when images aren't of equal size. I'll see how easy it looks, and if it's not so easy maybe I'll take the "always intersect" approach for this PR if that sounds good to you.
There was a problem hiding this comment.
I still like the idea of keeping the "extra pixels" around, but I think that will be easier to do after #83 is addressed so I'll continue intersecting for now.
| mms = mismatch(fixed, moving, mxshift; normalization=normalization) | ||
| best_i = indmin_mismatch(mms, thresh) | ||
| return best_i.I, ratio(mms[best_i], 0.0, Inf) | ||
| end | ||
|
|
||
| #returns new minbounds and maxbounds with range sizes change by fac | ||
| function scalebounds(minb, maxb, fac::Float64) |
There was a problem hiding this comment.
fac::Real? Float64 seems unnecessarily specific.
| #(Includes rotations, scaling, shear, etc) | ||
| #The space is centered on the identity matrix | ||
| function default_linmap_bounds(img::AbstractArray{T,N}) where {T, N} | ||
| d=0.05; nd=0.05; #not sure yet whether off-diagonals should be treated differently |
There was a problem hiding this comment.
I suspect these constants should be a keyword argument. Otherwise "reasonable default" might depend on how much movement one gets with this prep.
There was a problem hiding this comment.
Makes sense! Initially I wanted to try to phrase these in more understandable terms, letting the user specify maximum rotation, scaling, shear, etc but I found it difficult to combine those constraints. For now I guess I'll go with diagonal and off-diagonal kwargs.
| return ratio(mm, thresh, Inf) | ||
| end | ||
|
|
||
| #sets splits based on lower and upper bounds |
There was a problem hiding this comment.
👍 . I should bake this into QuadDIRECT (or rather, its successor) when I get to it.
| `kwargs...` can also include any other keyword argument that can be passed to `QuadDIRECT.analyze`. | ||
| It's recommended that you pass your own stopping criteria when possible (i.e. `rtol`, `atol`, and/or `fvalue`). | ||
|
|
||
| If you have a good initial guess at the solution, pass it with the `tfm0` kwarg to jump-start the search. |
| Therefore if you try to use the transform to warp an image of a size other than `size(moving)` then the result will not be as expected. | ||
| If you want to keep the rotation centered you must call `recenter(tfm, newctr)` where `newctr` is the displacement of the new center from the old center. | ||
| """ | ||
| function qd_affine(fixed, moving, mxshift, linmins, linmaxs, SD=eye(ndims(fixed)); thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), initial_tfm = IdentityTransformation(), kwargs...) |
There was a problem hiding this comment.
It's worth noting that affine doesn't technically require SD, since an affine transformation can accommodate the compensating transformation. I suspect we should still use it as a hint, however, because the intializations may be more "rigid-like" with it.
There was a problem hiding this comment.
I see your point, and as you say I think it's a good hint to have so that we can make the most of the user's intuitions about real-world bounds on movement.
| # a different set of errors than the tests above) | ||
| @testset "QuadDIRECT tests with standard images" begin | ||
| img = testimage("cameraman"); | ||
| img = Float64.(img) |
There was a problem hiding this comment.
Ah It used to be before I fixed the overflow in mismatch0, I can change it back.
There was a problem hiding this comment.
Wow so this was more complicated than I thought. After removing the Float64 conversion some tests failed to reach the fvalue criteria that we had in place. At first this seemed like a problem because those tests were passing with the same criteria before this PR. But after a lot of checks I became convinced that it was just a lucky accident that the tests in question passed with the 8-bit cameraman image before. What I realized is that with low precision sometimes the "coarse" step can achieve a lower objective value than the fine step. I'm pretty sure this is how that can happen:
The coarse step of qd_rigid involves two warps (a rotation and then a pixel-aligned shift) before computing the mismatch while the fine step only does a single warp (combined rotation and shift). Remember that neither method can reduce the mismatch to 0 because the "moving" image was created by interpolating the fixed image. It seems that with low-precision images sometimes the algorithm can find a pair of transforms in the coarse step that exhibit favorable rounding errors and achieve a lower mismatch than the single combined transform. In the cases where I've seen this happen the coarse step is already very close to the best answer so the fine step doesn't need to do anything. But it can make setting stopping criteria very confusing. After that rabbit hole my feeling is that it's simplest to just convert the images to higher precision before doing the registration.
There was a problem hiding this comment.
Hmm. The tests should be deterministic, though? There isn't any randomness I can see, and I don't think analyze has anything non-deterministic about it. (I could be forgetting, however.)
I wonder if we might be seeing overflow in warp? (i.e., Interpolations?) So that very subtly different ways of computing things yield quite different results?
There was a problem hiding this comment.
It is deterministic. If I run this PR without converting to Float64 I always get the same result. The reason it's inconsistent with Chantal's PR is that I've since changed some things about how the function works. The change may be due to my correcting the order of more of the compositions, though I haven't tracked it down to a single line. But I do know that the mismatch value shown in the log of her Travis run is lower than expected even when applying the perfect inverse transform to moving, consistent with my idea that the test passed with 8-bit images only by a happy accident.
There was a problem hiding this comment.
And all of this could also be consistent with overflow within warp, good idea!
There was a problem hiding this comment.
Okay, I got it! The problem is with warp, but it's not overflow. Since we can't represent NaN with fixed point numbers warp imputes zeros when sampling outside of the image. This can change the normalizing term in the mismatch calculation dramatically, especially for small images. Here's an example:
using BlockRegistration, RegisterMismatch, FixedPointNumbers, CoordinateTransformations
a = fill(N0f16(0.5), 5,5)
b = Float64.(a)
tfm = Translation(0.5,0.5)
wa, fa = RegisterOptimize.warp_and_intersect(a, a, tfm)
wb, fb = RegisterOptimize.warp_and_intersect(b, b, tfm)
@show ratio(mismatch0(wa, fa), 0) #nonzero
@show ratio(mismatch0(wb, fb), 0) #zeroThere may be a nice Julian way to fix this. Maybe by defining an efficient element type that wraps a fixed point number and can handle missing values? Or maybe the new missing value from Julia 0.7? I haven't played with that yet but my guess is that performance may suffer. Maybe easiest would be just to never try to warp an array of fixed point numbers, forcing conversion beforehand?
There was a problem hiding this comment.
Great analysis! Yes, this makes it much clearer, and a good example of how important NaN-poisoning is. (Note that most transformation frameworks, e.g., ANTs and matlab, fill with zeros.) Now I have no objection to coercing it. Maybe better it's to use float.(img), though; the latter will even work for RGB, I think, whereas Float64.(img) will not.
I'm filing an issue (#82) to remind myself to go through BlockRegistration more generally and see whether there are other places we should be doing this.
There was a problem hiding this comment.
Great, and I think I will use float(img) rather than float.(img) because the former is smart enough not to allocate a new array when it already has element type Float64
| inds = intersect.(indices(img), indices(img2)) | ||
| fixed = img[inds...] | ||
| moving = img2[inds...] | ||
| return fixed, moving |
There was a problem hiding this comment.
weird indentation here and below
|
Continued in #85 |
Also added new tests and refactored existing ones. You'll notice that some new
qd_affinetests are commented out because they incorporate randomness and they fail sometimes. Theqd_translateandqd_rigidtests with randomness pass reliably so I kept those. We may be able to switch the affine ones back on in the future with some newer version of QuadDIRECT.