From 084201622d6f20185625ab9ec082dc30bc89be7a Mon Sep 17 00:00:00 2001 From: Cody-G Date: Wed, 4 Jul 2018 17:52:54 -0500 Subject: [PATCH 1/3] Allow specifying initial transformations with all quaddirect-based registration methods. In some cases ordering of transform compositions was also incorrect, this fixes that. --- src/qd_rigid.jl | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/qd_rigid.jl b/src/qd_rigid.jl index 940ce48..a7dba7d 100644 --- a/src/qd_rigid.jl +++ b/src/qd_rigid.jl @@ -30,8 +30,8 @@ function rot(theta, img, SD=eye(ndims(img))) end #Finds the best shift aligning moving to fixed, possibly after an initial transformation `intial_tfm` -function best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=-1) - if initial_tfm != -1 +function best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=IdentityTransformation()) + if initial_tfm != IdentityTransformation() newmov = warp(moving, initial_tfm) inds1, inds2 = indices(fixed), indices(newmov) inds = ([intersect(x, y) for (x,y) in zip(inds1, inds2)]...) @@ -44,22 +44,21 @@ function best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, in end #rotation + shift, slow because it warps for every rotation and shift -function slow_mm(tfm, fixed, moving, thresh, SD; initial_tfm = -1) +function slow_mm(tfm, fixed, moving, thresh, SD; initial_tfm = IdentityTransformation()) tfm = tfmx(tfm, moving, SD) - if initial_tfm != -1 - tfm = tfm ∘ initial_tfm #compose with rotation already computed - end + tfm = initial_tfm ∘ tfm#compose with rotation already computed newmov = warp(moving, tfm) inds1, inds2 = indices(fixed), indices(newmov) inds = ([intersect(x, y) for (x,y) in zip(inds1, inds2)]...) - mm = mismatch0(view(fixed, inds...), view(newmov, inds...); normalization=:pixels) + mm = mismatch0(view(fixed, inds...), newmov[inds...]; normalization=:intensity) rslt = ratio(mm, thresh, Inf) return rslt end #rotation + shift, fast because it uses fourier method for shift -function fast_mm(theta, mxshift, fixed, moving, thresh, SD) +function fast_mm(theta, mxshift, fixed, moving, thresh, SD; initial_tfm = IdentityTransformation()) tfm = rot(theta, moving, SD) + tfm = initial_tfm ∘ tfm#compose with rotation already computed newmov = warp(moving, tfm) inds1, inds2 = indices(fixed), indices(newmov) inds = ([intersect(x, y) for (x,y) in zip(inds1, inds2)]...) @@ -67,20 +66,20 @@ 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))])), kwargs...) - f(x) = fast_mm(x, mxshift, fixed, moving, thresh, SD) +function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; initial_tfm = IdentityTransformation(), thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), kwargs...) + f(x) = fast_mm(x, mxshift, fixed, moving, thresh, SD; initial_tfm = initial_tfm) #note: if a trial rotation results in image overlap < thresh for all possible shifts then QuadDIRECT throws an error upper = [mxrot...] lower = -upper splits = ([[-x; 0.0; x] for x in mxrot]...) root_coarse, x0coarse = analyze(f, splits, lower, upper; maxevals=10^4, minwidth=minwidth_rot, print_interval=100, kwargs...) box_coarse = minimum(root_coarse) - tfmcoarse0 = rot(position(box_coarse, x0coarse), moving) + tfmcoarse0 = initial_tfm ∘ rot(position(box_coarse, x0coarse), moving) best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm = tfmcoarse0) tfmcoarse = tfmcoarse0 ∘ Translation(best_shft) return tfmcoarse, mm end -function qd_rigid_fine(fixed, moving, mxrot, minwidth, SD; initial_tfm = -1, thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), kwargs...) +function qd_rigid_fine(fixed, moving, mxrot, minwidth, SD; initial_tfm = IdentityTransformation(), thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), kwargs...) f2(x) = slow_mm(x, fixed, moving, thresh, SD; initial_tfm = initial_tfm) upper_shft = fill(1, ndims(fixed)) upper_rot = mxrot @@ -92,8 +91,7 @@ function qd_rigid_fine(fixed, moving, mxrot, minwidth, SD; initial_tfm = -1, thr minwidth = vcat(minwidth_shfts, minwidth_rots) root, x0 = analyze(f2, splits, lower, upper; maxevals=10^4, minwidth=minwidth, print_interval=100, kwargs...) box = minimum(root) - params = position(box, x0) - tfmfine = tfmx(position(box, x0), moving) + tfmfine = initial_tfm ∘ tfmx(position(box, x0), moving) return tfmfine, value(box) end @@ -109,12 +107,11 @@ 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))])), kwargs...) +function qd_rigid(fixed, moving, mxshift, mxrot, minwidth_rot, SD=eye(ndims(fixed)); thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), tfm0 = IdentityTransformation(), kwargs...) mxrot = [mxrot...] print("Running coarse step\n") - tfm_coarse, mm_coarse = qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; thresh = thresh, kwargs...) + tfm_coarse, mm_coarse = qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot, SD; initial_tfm = tfm0, thresh = thresh, 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, kwargs...) - final_tfm = tfm_fine ∘ tfm_coarse + final_tfm, mm_fine = qd_rigid_fine(fixed, moving, mxrot./10, minwidth_rot, SD; initial_tfm = tfm_coarse, thresh = thresh, kwargs...) return final_tfm, mm_fine end From 4e1f8bd407ffbefb649ee0809a82228fc0989592 Mon Sep 17 00:00:00 2001 From: Cody-G Date: Fri, 6 Jul 2018 15:56:14 -0500 Subject: [PATCH 2/3] Test for overflow with mismatch0 (test fails) --- test/register_mismatch.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/register_mismatch.jl b/test/register_mismatch.jl index d35e7ad..4287f84 100644 --- a/test/register_mismatch.jl +++ b/test/register_mismatch.jl @@ -89,6 +89,15 @@ for imsz in ((7,10), (6,5)) end end +# Test for denom overflow with mismatch0 +C16 = N0f16[0.6 0.1; 0.7 0.1] +D16 = N0f16[0.7 0.1; 0.6 0.1] +C = Float64.(C16) +D = Float64.(D16) +nd = RegisterMismatch.mismatch0(C16, D16; normalization=:intensity) +@test nd.denom ≈ sum((C.^2).+(D.^2)) +@test nd.num ≈ sum((C.-D).^2) + # Compare direct, global, and apertured C = rand(7,9) D = rand(7,9) From 27f8c2b670772e549df2422bedaefb6ede5828c9 Mon Sep 17 00:00:00 2001 From: Cody-G Date: Fri, 6 Jul 2018 15:56:44 -0500 Subject: [PATCH 3/3] Fix overflow bug with mismatch0 (test passes now) --- src/RegisterMismatchCommon.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/RegisterMismatchCommon.jl b/src/RegisterMismatchCommon.jl index 76593dc..346b436 100644 --- a/src/RegisterMismatchCommon.jl +++ b/src/RegisterMismatchCommon.jl @@ -102,11 +102,14 @@ nanval{T}(::Type{T}) = convert(Float32, NaN) """ function mismatch0{Tf,Tm,N}(fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm,N}; normalization = :intensity) size(fixed) == size(moving) || throw(DimensionMismatch("Size $(size(fixed)) of fixed is not equal to size $(size(moving)) of moving")) - num = denom = zero(promote_type(typeof((oneunit(Tf) - oneunit(Tm))^2), typeof(oneunit(Tf)^2+oneunit(Tm)^2))) + _mismatch0(zero(Float64), zero(Float64), fixed, moving; normalization=normalization) +end + +function _mismatch0{T,Tf,Tm,N}(num::T, denom::T, fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm,N}; normalization = :intensity) if normalization == :intensity for i in eachindex(fixed, moving) - vf = fixed[i] - vm = moving[i] + vf = T(fixed[i]) + vm = T(moving[i]) if isfinite(vf) && isfinite(vm) num += (vf-vm)^2 denom += vf^2 + vm^2 @@ -114,8 +117,8 @@ function mismatch0{Tf,Tm,N}(fixed::AbstractArray{Tf,N}, moving::AbstractArray{Tm end elseif normalization == :pixels for i in eachindex(fixed, moving) - vf = fixed[i] - vm = moving[i] + vf = T(fixed[i]) + vm = T(moving[i]) if isfinite(vf) && isfinite(vm) num += (vf-vm)^2 denom += 1