Skip to content

Cjg/initial tfm #80

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
timholy merged 3 commits into from
Jul 6, 2018
Merged

Cjg/initial tfm #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0842016
Allow specifying initial transformations with all quaddirect-based re…
Cody-G Jul 4, 2018
4e1f8bd
Test for overflow with mismatch0 (test fails)
Cody-G Jul 6, 2018
27f8c2b
Fix overflow bug with mismatch0 (test passes now)
Cody-G Jul 6, 2018
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
13 changes: 8 additions & 5 deletions src/RegisterMismatchCommon.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,20 +102,23 @@ 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
end
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
Expand Down
33 changes: 15 additions & 18 deletions src/qd_rigid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)]...)
Expand All @@ -44,43 +44,42 @@ 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())
Copy link
Copy Markdown
Member

@timholy timholy Jul 6, 2018

Choose a reason for hiding this comment

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

Even though it's effectively the same thing I like this better.

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)]...)
bshft, mm = best_shift(view(fixed, inds...), newmov[inds...], mxshift, thresh; normalization=:intensity)
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
Expand All @@ -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

Expand All @@ -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
9 changes: 9 additions & 0 deletions test/register_mismatch.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down