Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/qd_rigid.jl
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh
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)
tfmcoarse = Translation(best_shft) ∘ tfmcoarse0
tfmcoarse = tfmcoarse0 ∘ Translation(best_shft)
Copy link
Contributor

Choose a reason for hiding this comment

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

After thinking on it I'm less certain that the new ordering is correct. The mm value that we're using to choose a transform was computed by first rotating the image here and then finding the best translation of the rotated image here. Since tfmcoarse0 is the rotation then I think it should be applied first, i.e. Translation(best_shft) ∘ tfmcoarse0. Or am I missing something? Note that we redundantly compute the best shift again here (the redundancy is bad design on my part) but in both cases the rotation is applied first.

Copy link
Member

Choose a reason for hiding this comment

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

I agree this is confusing; and yes, the test fails (badly) without this. Here's how I now think about it: if m is the moving image and t1 is the transform that you rotate it with, you produce a new image m1. Mathematically, we can think of this as defining a function (if you think about images as being continuous rather than discrete) as

m1(x) = m(t1(x))

because in warp the transformation operates on the coordinates at which you evaluate the image. Now we take m1 as the new moving image and calculate a new transformation t2 (which happens to be a translation). At the end we get m2 which is defined

m2(x) = m1(t2(x)) = m(t1(t2(x)))

The key thing to note is that we have t1(t2(x)) = (t1 ∘ t2)(x). So it's the opposite of what I would have initially expected from the order in which these transformations were defined, which is probably why I didn't catch this when I reviewed #75.

Copy link
Contributor

Choose a reason for hiding this comment

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

Oh wow!! I'm surprised that I made it this long without reckoning with this misunderstanding... I think I inadvertently put a lot of effort into working around it. Thanks for that explanation, it's not intuitive but it makes perfect sense now!

return tfmcoarse, mm
end

Expand Down
40 changes: 38 additions & 2 deletions test/register_optimize.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
using StaticArrays, AffineTransforms, Interpolations, Base.Test
import BlockRegistration, RegisterOptimize
using RegisterCore, RegisterPenalty, RegisterDeformation, RegisterMismatch, RegisterFit
using Images, CoordinateTransformations, Rotations, RegisterOptimize
using Images, CoordinateTransformations, Rotations, RegisterOptimize, TestImages

using RegisterTestUtilities

Expand Down Expand Up @@ -334,7 +334,7 @@ b = AffineTransforms.transform(a, tformtranslate([2.0;0.0]) * tformrotate(pi/6))
tfm0 = tformtranslate([-2.0;0.0]) * tformrotate(-pi/6)
#note: maxshift must be GREATER than the true shift in order to find the true shift
tfm, mm = rotation_gridsearch(a, b, [11;11], [pi/6], [11])
@assert tfm.offset == tfm0.offset
@assert tfm.offset == tfm0.offset
@assert tfm.scalefwd == tfm0.scalefwd

## 3D
Expand Down Expand Up @@ -382,3 +382,39 @@ tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh=thres

@test mm < 1e-4
@test sum(abs.(vcat(tfm0.m[:], tfm0.v) - vcat(RotXYZ(tfm.m)[:], tfm.v))) < 0.1




# tests with standard images
@testset "tests with standard images" begin
img = testimage("cameraman");

tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now
img2 = warp(img,tfm)


inds = intersect.(indices(img), indices(img2))
img = img[inds...]
img2 = img2[inds...]

fixed = img

mxshift = (100,100) #make sure this isn't too small
mxrot = (0.5,)
minwidth_rot = fill(0.002, 3)
SD = diagm([pixelspacing(fixed)...])

moving = img2

tform, mm = qd_rigid(fixed, moving, mxshift, mxrot, minwidth_rot, SD, rtol=0, fvalue=0.01)

# imgw = warp(img2, tform)
# inds2 = intersect.(indices(img), indices(imgw))
# imshow(colorview(RGB,img[inds2...],imgw[inds2...],zeroarray))

testimg = tfm ∘ tform #should be the identity matrix
@test 0.99 < testimg.m[1,1] < 1.09 && 0.99 < testimg.m[2,2] < 1.09 && abs(testimg.m[1,2]) < 0.01 && abs(testimg.m[2,1]) < 0.01
@test abs(testimg.v[1]) < 2 && abs(testimg.v[2]) < 2

end #tests with standard images