From 332c761489bbb72e9e459d88a429d3b8a147b478 Mon Sep 17 00:00:00 2001 From: Dae Woo Kim Date: Tue, 26 May 2026 14:49:31 -0500 Subject: [PATCH] Apply runic formatting Co-Authored-By: Claude Opus 4.7 --- src/RegisterQD.jl | 14 ++-- src/affine.jl | 120 +++++++++++++++++++--------------- src/gridsearch.jl | 26 ++++---- src/rigid.jl | 88 ++++++++++++++----------- src/translations.jl | 44 +++++++------ src/util.jl | 77 +++++++++++----------- test/gpu_tests.jl | 88 ++++++++++++------------- test/gridsearch.jl | 16 ++--- test/initial_tfm.jl | 140 +++++++++++++++++++-------------------- test/qd_random.jl | 156 ++++++++++++++++++++++---------------------- test/qd_standard.jl | 70 ++++++++++---------- test/util.jl | 34 +++++----- 12 files changed, 452 insertions(+), 421 deletions(-) diff --git a/src/RegisterQD.jl b/src/RegisterQD.jl index 90ee40d..3900329 100644 --- a/src/RegisterQD.jl +++ b/src/RegisterQD.jl @@ -21,12 +21,12 @@ include("affine.jl") include("gridsearch.jl") export qd_translate, - qd_rigid, - qd_affine, - arrayscale, - grid_rotations, - rotation_gridsearch, - getSD, - qsmooth + qd_rigid, + qd_affine, + arrayscale, + grid_rotations, + rotation_gridsearch, + getSD, + qsmooth end # module diff --git a/src/affine.jl b/src/affine.jl index 871513f..1e9a954 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -9,12 +9,12 @@ If `initial_tfm` is supplied, it will be composed with the one from `mat`. See also [`RegisterQD.aff`](@ref), which includes a translational component. """ -function linmap(mat, img::AbstractArray{T,N}, initial_tfm=IdentityTransformation()) where {T,N} +function linmap(mat, img::AbstractArray{T, N}, initial_tfm = IdentityTransformation()) where {T, N} # img isn't used *except* to get the dimensionality in an inferrable way, so that # the (static) size of the array in `tfm` is known. mat = [mat...] - mat = reshape(mat, N,N) - lm = LinearMap(SMatrix{N,N}(mat)) + mat = reshape(mat, N, N) + lm = LinearMap(SMatrix{N, N}(mat)) # The order below is subtle: you might think that `initial_tfm` should be applied # to `x` first, so that the order of composition should be `lm ∘ initial_tfm`. # But what we're going to be doing is calculating an `lm` that aligns @@ -40,16 +40,16 @@ and the last `N^2` are for the linear portion (see [`RegisterQD.linmap`](@ref)). If `initial_tfm` is supplied, it will be composed with the one from `params`. """ -function aff(params, img::AbstractArray{T,N}, initial_tfm=IdentityTransformation()) where {T,N} +function aff(params, img::AbstractArray{T, N}, initial_tfm = IdentityTransformation()) where {T, N} # img isn't used *except* to get the dimensionality in an inferrable way, so that # the (static) sizes of the arrays in `tfm` are known. params = [params...] - length(params) == (N+N^2) || throw(DimensionMismatch("expected $(N+N^2) parameters, got $(length(params))")) + length(params) == (N + N^2) || throw(DimensionMismatch("expected $(N + N^2) parameters, got $(length(params))")) offs = Float64.(params[1:N]) - mat = Float64.(params[(N+1):end]) - mat = reshape(mat,N,N) + mat = Float64.(params[(N + 1):end]) + mat = reshape(mat, N, N) # The order below is explained in `linmap` - return initial_tfm ∘ AffineMap(SMatrix{N,N}(mat), SVector{N}(offs)) + return initial_tfm ∘ AffineMap(SMatrix{N, N}(mat), SVector{N}(offs)) end # This acts as the objective function during optimization, see qd_affine_coarse @@ -57,20 +57,20 @@ end # here params contains parameters of a linear map # The "fast" part means that it handles translation using `best_shift`, so this is # concerned only with the linear-map portion of the affine transformation. -function affine_mm_fast(params, mxshift, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation()) +function affine_mm_fast(params, mxshift, fixed, moving, thresh, SD; initial_tfm = IdentityTransformation()) tfm = arrayscale(linmap(params, moving, initial_tfm), SD) moving, fixed = warp_and_intersect(moving, fixed, tfm) - bshft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity) + bshft, mm = best_shift(fixed, moving, mxshift, thresh; normalization = :intensity) return mm end # Similar to the above but includes the translation in `params` # This is used for optimization of the complete affine transformation. -function affine_mm_slow(params, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation()) +function affine_mm_slow(params, fixed, moving, thresh, SD; initial_tfm = IdentityTransformation()) tfm = arrayscale(aff(params, moving, initial_tfm), SD) moving, fixed = warp_and_intersect(moving, fixed, tfm) - mm = mismatch_zeroshift(fixed, moving; normalization=:intensity) - return ratio(mm, thresh; fillval=Inf) + mm = mismatch_zeroshift(fixed, moving; normalization = :intensity) + return ratio(mm, thresh; fillval = Inf) end """ @@ -81,26 +81,30 @@ Compute an affine transformation `tfm` that coarsely minimizes the mismatch betw The translational component is handled at the level of integer-pixel shifts, which is why this is "coarse" optimization. """ -function qd_affine_coarse(fixed, moving, mxshift, linmins, linmaxs; - SD=I, - initial_tfm=IdentityTransformation(), - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - minwidth=default_lin_minwidths(moving), - maxevals=5e4, - kwargs...) +function qd_affine_coarse( + fixed, moving, mxshift, linmins, linmaxs; + SD = I, + initial_tfm = IdentityTransformation(), + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + minwidth = default_lin_minwidths(moving), + maxevals = 5.0e4, + kwargs... + ) # Define the objective function for the coarse stage - f(x) = affine_mm_fast(x, mxshift, fixed, moving, thresh, SD; initial_tfm=initial_tfm) + f(x) = affine_mm_fast(x, mxshift, fixed, moving, thresh, SD; initial_tfm = initial_tfm) # QuadDIRECT benefits from bounds on the search region, build those bounds upper = linmaxs lower = linmins # Try to find the global minimum (within the needed precision) - root, x0 = _analyze(f, lower, upper; - minwidth=minwidth, print_interval=100, maxevals=maxevals, kwargs..., atol=0, rtol=1e-3) + root, x0 = _analyze( + f, lower, upper; + minwidth = minwidth, print_interval = 100, maxevals = maxevals, kwargs..., atol = 0, rtol = 1.0e-3 + ) box = minimum(root) # Convert the coarse-minimum into a complete affine transformation params = position(box, x0) tfmcoarse0 = linmap(params, moving, initial_tfm) - best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=arrayscale(tfmcoarse0, SD)) + best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization = :intensity, initial_tfm = arrayscale(tfmcoarse0, SD)) tfmcoarse = tfmcoarse0 ∘ pscale(Translation(best_shft), SD) return tfmcoarse, mm end @@ -120,15 +124,17 @@ distinguishes this from the full-search-space `minwidth` accepted by [`RegisterQD.qd_affine_coarse`](@ref): here it names just the linear-map subspace, which is combined internally with the (fixed) translation resolution. """ -function qd_affine_fine(fixed, moving, linmins, linmaxs; - SD=I, - initial_tfm=IdentityTransformation(), - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - minwidth_mat=default_lin_minwidths(fixed)./10, - maxevals=5e4, - kwargs...) +function qd_affine_fine( + fixed, moving, linmins, linmaxs; + SD = I, + initial_tfm = IdentityTransformation(), + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + minwidth_mat = default_lin_minwidths(fixed) ./ 10, + maxevals = 5.0e4, + kwargs... + ) # Define the objective function - f(x) = affine_mm_slow(x, fixed, moving, thresh, SD; initial_tfm=initial_tfm) + f(x) = affine_mm_slow(x, fixed, moving, thresh, SD; initial_tfm = initial_tfm) # Limit translations to only a couple of pixels upper_shft = fill(2.0, ndims(fixed)) # Combine translation limits with the affine-transform limits @@ -138,8 +144,10 @@ function qd_affine_fine(fixed, moving, linmins, linmaxs; minwidth_shfts = fill(0.01, ndims(fixed)) minwidth = vcat(minwidth_shfts, minwidth_mat) # Try to find the global minimum - root, x0 = _analyze(f, lower, upper; - minwidth=minwidth, print_interval=100, maxevals=maxevals, kwargs...) + root, x0 = _analyze( + f, lower, upper; + minwidth = minwidth, print_interval = 100, maxevals = maxevals, kwargs... + ) box = minimum(root) # Convert the result to an affine transformation params = position(box, x0) @@ -190,16 +198,18 @@ The default value for `thresh` is 50% of the sum-of-squared-intensity of `fixed` used by `qd_translate` and `qd_rigid` because affine transformations have additional degrees of freedom (scaling and shear) that make degenerate low-overlap solutions more likely. """ -function qd_affine(fixed, moving, mxshift, linmins, linmaxs; - presmoothed=false, - SD=I, - # Affine's extra degrees of freedom (scaling/shear) make degenerate - # low-overlap solutions more likely, so the default requires more overlap - # (50%) than the 10% used by qd_translate/qd_rigid. - thresh=0.5*sum(abs2.(fixed[.!(isnan.(fixed))])), - initial_tfm=IdentityTransformation(), - print_interval=100, - kwargs...) +function qd_affine( + fixed, moving, mxshift, linmins, linmaxs; + presmoothed = false, + SD = I, + # Affine's extra degrees of freedom (scaling/shear) make degenerate + # low-overlap solutions more likely, so the default requires more overlap + # (50%) than the 10% used by qd_translate/qd_rigid. + thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))])), + initial_tfm = IdentityTransformation(), + print_interval = 100, + kwargs... + ) fixed, moving = float(fixed), float(moving) if presmoothed moving = qinterp(eltype(fixed), moving) @@ -208,19 +218,25 @@ function qd_affine(fixed, moving, mxshift, linmins, linmaxs; linmaxs = [linmaxs...] print_interval < typemax(Int) && print("Running coarse step\n") mw = default_lin_minwidths(moving) - tfm_coarse, mm_coarse = qd_affine_coarse(fixed, moving, mxshift, linmins, linmaxs; - SD = SD, minwidth=mw, initial_tfm=initial_tfm, thresh=thresh, print_interval=print_interval, kwargs...) + tfm_coarse, mm_coarse = qd_affine_coarse( + fixed, moving, mxshift, linmins, linmaxs; + SD = SD, minwidth = mw, initial_tfm = initial_tfm, thresh = thresh, print_interval = print_interval, kwargs... + ) print_interval < typemax(Int) && print("Running fine step\n") - mw = mw./100 + mw = mw ./ 100 linmins, linmaxs = scalebounds(linmins, linmaxs, 0.5) # should these be narrowed further? - final_tfm, final_mm = qd_affine_fine(fixed, moving, linmins, linmaxs; - SD = SD, minwidth_mat=mw, initial_tfm=tfm_coarse, thresh=thresh, print_interval=print_interval, kwargs...) + final_tfm, final_mm = qd_affine_fine( + fixed, moving, linmins, linmaxs; + SD = SD, minwidth_mat = mw, initial_tfm = tfm_coarse, thresh = thresh, print_interval = print_interval, kwargs... + ) return final_tfm, final_mm end -function qd_affine(fixed, moving, mxshift; - dmax = 0.05, ndmax = 0.05, - kwargs...) +function qd_affine( + fixed, moving, mxshift; + dmax = 0.05, ndmax = 0.05, + kwargs... + ) minb, maxb = default_linmap_bounds(fixed; dmax = dmax, ndmax = ndmax) return qd_affine(fixed, moving, mxshift, minb, maxb; kwargs...) end diff --git a/src/gridsearch.jl b/src/gridsearch.jl index 078be26..ff75d09 100644 --- a/src/gridsearch.jl +++ b/src/gridsearch.jl @@ -18,36 +18,36 @@ pixel spacing. function grid_rotations(maxradians, rgridsz, SD) rgridsz = [rgridsz...] maxradians = [maxradians...] - nd = size(SD,1) + nd = size(SD, 1) if !all(isodd, rgridsz) @warn("rgridsz should be odd; rounding up to the next odd integer") end - for i = 1:length(rgridsz) + for i in 1:length(rgridsz) if !isodd(rgridsz[i]) rgridsz[i] = max(round(Int, rgridsz[i]) + 1, 1) end end - grid_radius = map(x->div(x,2), rgridsz) + grid_radius = map(x -> div(x, 2), rgridsz) if nd > 2 - gridcoords = [range(-grid_radius[x]*maxradians[x], stop=grid_radius[x]*maxradians[x], length=rgridsz[x]) for x=1:nd] + gridcoords = [range(-grid_radius[x] * maxradians[x], stop = grid_radius[x] * maxradians[x], length = rgridsz[x]) for x in 1:nd] rotation_angles = Iterators.product(gridcoords...) else - rotation_angles = range(-grid_radius[1]*maxradians[1], stop=grid_radius[1]*maxradians[1], length=rgridsz[1]) + rotation_angles = range(-grid_radius[1] * maxradians[1], stop = grid_radius[1] * maxradians[1], length = rgridsz[1]) end - axmat = Matrix{Float64}(I,nd,nd) - axs = map(x->axmat[:,x], 1:nd) + axmat = Matrix{Float64}(I, nd, nd) + axs = map(x -> axmat[:, x], 1:nd) tfeye = tformeye(nd) output = typeof(tfeye)[] for ra in rotation_angles if nd > 2 - euler_rots = map(x->tformrotate(x...), zip(axs, ra)) - rot = foldr(∘, euler_rots, init=tfeye) + euler_rots = map(x -> tformrotate(x...), zip(axs, ra)) + rot = foldr(∘, euler_rots, init = tfeye) elseif nd == 2 rot = tformrotate(ra) else error("Unsupported dimensionality") end - push!(output, AffineMap(SD*rot.linear/SD , zeros(nd))) #account for sample spacing + push!(output, AffineMap(SD * rot.linear / SD, zeros(nd))) #account for sample spacing end return output end @@ -61,7 +61,7 @@ best rotation and shift out of the values searched, along with the mismatch valu For more on how the arguments `maxradians`, `rgridsz`, and `SD` influence the search, see the documentation for `grid_rotations`. """ -function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD = Matrix{Float64}(I,ndims(fixed),ndims(fixed))) +function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD = Matrix{Float64}(I, ndims(fixed), ndims(fixed))) rgridsz = [rgridsz...] nd = ndims(moving) @assert nd == ndims(fixed) @@ -74,9 +74,9 @@ function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz; SD = #calc mismatch #mm = mismatch(fixed, new_moving, maxshift; normalization=:pixels) mm = mismatch(fixed, new_moving, maxshift) - thresh = 0.1*maximum(x->x.denom, mm) + thresh = 0.1 * maximum(x -> x.denom, mm) best_i = argmin_mismatch(mm, thresh) - cur_best =ratio(mm[best_i], 0.0) + cur_best = ratio(mm[best_i], 0.0) if cur_best < best_mm best_mm = cur_best best_rot = rot diff --git a/src/rigid.jl b/src/rigid.jl index 9965761..32309dc 100644 --- a/src/rigid.jl +++ b/src/rigid.jl @@ -24,9 +24,13 @@ function rot(qx::Number, qy::Number, qz::Number) r2 = qx^2 + qy^2 + qz^2 r2 > 1 && return nothing qr = sqrt(1 - r2) - R = @SMatrix([1 - 2*(qy^2 + qz^2) 2*(qx*qy - qz*qr) 2*(qx*qz + qy*qr); - 2*(qx*qy + qz*qr) 1 - 2*(qx^2 + qz^2) 2*(qy*qz - qx*qr); - 2*(qx*qz - qy*qr) 2*(qy*qz + qx*qr) 1 - 2*(qx^2 + qy^2)]) + R = @SMatrix( + [ + 1 - 2 * (qy^2 + qz^2) 2 * (qx * qy - qz * qr) 2 * (qx * qz + qy * qr); + 2 * (qx * qy + qz * qr) 1 - 2 * (qx^2 + qz^2) 2 * (qy * qz - qx * qr); + 2 * (qx * qz - qy * qr) 2 * (qy * qz + qx * qr) 1 - 2 * (qx^2 + qy^2) + ] + ) return LinearMap(RotMatrix(R)) end @@ -40,24 +44,24 @@ tfmrigid(dx::Number, dy::Number, dz::Number, qx::Number, qy::Number, qz::Number) @noinline dimcheck(v, lv) = length(v) == lv || throw(DimensionMismatch("expected $lv parameters, got $(length(v))")) -function rot(theta, img::AbstractArray{T,2}) where {T} +function rot(theta, img::AbstractArray{T, 2}) where {T} dimcheck(theta, 1) return rot(theta[1]) end -function rot(qs, img::AbstractArray{T,3}) where {T} +function rot(qs, img::AbstractArray{T, 3}) where {T} dimcheck(qs, 3) qx, qy, qz = qs return rot(qx, qy, qz) end -function tfmrigid(params, img::AbstractArray{T,2}) where {T} +function tfmrigid(params, img::AbstractArray{T, 2}) where {T} dimcheck(params, 3) dx, dy, θ = params return tfmrigid(dx, dy, θ) end -function tfmrigid(params, img::AbstractArray{T,3}) where {T} +function tfmrigid(params, img::AbstractArray{T, 3}) where {T} dimcheck(params, 6) dx, dy, dz, qx, qy, qz = params return tfmrigid(dx, dy, dz, qx, qy, qz) @@ -66,54 +70,60 @@ end ## Computing mismatch #rotation + shift, fast because it uses fourier method for shift -function rigid_mm_fast(theta, mxshift, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation()) +function rigid_mm_fast(theta, mxshift, fixed, moving, thresh, SD; initial_tfm = IdentityTransformation()) # The reason this is `initial_tfm ∘ rot` rather than `rot ∘ initial_tfm` # is explained in `linmap` tfm = arrayscale(initial_tfm ∘ rot(theta, moving), SD) moving, fixed = warp_and_intersect(moving, fixed, tfm) - bshft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity) + bshft, mm = best_shift(fixed, moving, mxshift, thresh; normalization = :intensity) return mm end #rotation + shift, slow because it warps for every rotation and shift -function rigid_mm_slow(params, fixed, moving, thresh, SD; initial_tfm=IdentityTransformation()) +function rigid_mm_slow(params, fixed, moving, thresh, SD; initial_tfm = IdentityTransformation()) tfm = arrayscale(initial_tfm ∘ tfmrigid(params, moving), SD) moving, fixed = warp_and_intersect(moving, fixed, tfm) - mm = mismatch_zeroshift(fixed, moving; normalization=:intensity) - return ratio(mm, thresh; fillval=Inf) + mm = mismatch_zeroshift(fixed, moving; normalization = :intensity) + return ratio(mm, thresh; fillval = Inf) end -function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot; - SD=I, - initial_tfm=IdentityTransformation(), - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - kwargs...) +function qd_rigid_coarse( + fixed, moving, mxshift, mxrot, minwidth_rot; + SD = I, + initial_tfm = IdentityTransformation(), + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + kwargs... + ) #note: if a trial rotation results in image overlap < thresh for all possible shifts then QuadDIRECT throws an error - f(x) = rigid_mm_fast(x, mxshift, fixed, moving, thresh, SD; initial_tfm=initial_tfm) + f(x) = rigid_mm_fast(x, mxshift, fixed, moving, thresh, SD; initial_tfm = initial_tfm) upper = [mxrot...] lower = -upper - root_coarse, x0coarse = _analyze(f, lower, upper; - minwidth=minwidth_rot, print_interval=100, maxevals=5e4, kwargs..., atol=0, rtol=1e-3) + root_coarse, x0coarse = _analyze( + f, lower, upper; + minwidth = minwidth_rot, print_interval = 100, maxevals = 5.0e4, kwargs..., atol = 0, rtol = 1.0e-3 + ) box_coarse = minimum(root_coarse) tfmcoarse0 = initial_tfm ∘ rot(position(box_coarse, x0coarse), moving) - best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=arrayscale(tfmcoarse0, SD)) + best_shft, mm = best_shift(fixed, moving, mxshift, thresh; normalization = :intensity, initial_tfm = arrayscale(tfmcoarse0, SD)) tfmcoarse = tfmcoarse0 ∘ pscale(Translation(best_shft), SD) return tfmcoarse, mm end -function qd_rigid_fine(fixed, moving, mxrot, minwidth_rot; - SD=I, - initial_tfm=IdentityTransformation(), - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - kwargs...) - f(x) = rigid_mm_slow(x, fixed, moving, thresh, SD; initial_tfm=initial_tfm) +function qd_rigid_fine( + fixed, moving, mxrot, minwidth_rot; + SD = I, + initial_tfm = IdentityTransformation(), + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + kwargs... + ) + f(x) = rigid_mm_slow(x, fixed, moving, thresh, SD; initial_tfm = initial_tfm) upper_shft = fill(2.0, ndims(fixed)) upper_rot = mxrot upper = vcat(upper_shft, upper_rot) lower = -upper minwidth_shfts = fill(0.005, ndims(fixed)) minwidth = vcat(minwidth_shfts, minwidth_rot) - root, x0 = _analyze(f, lower, upper; minwidth=minwidth, print_interval=100, maxevals=5e4, kwargs...) + root, x0 = _analyze(f, lower, upper; minwidth = minwidth, print_interval = 100, maxevals = 5.0e4, kwargs...) box = minimum(root) tfmfine = initial_tfm ∘ tfmrigid(position(box, x0), moving) return tfmfine, value(box) @@ -179,14 +189,16 @@ Both the output `tfm` and any `initial_tfm` are represented in *physical* coordi as long as `initial_tfm` is a rigid transformation, `tfm` will be a pure rotation+translation. If `SD` is not the identity, use `arrayscale` before applying the result to `moving`. """ -function qd_rigid(fixed, moving, mxshift::VecLike, mxrot::Union{Number,VecLike}; - presmoothed=false, - SD=I, - minwidth_rot=default_minwidth_rot(CartesianIndices(fixed), SD), - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - initial_tfm=IdentityTransformation(), - print_interval=100, - kwargs...) +function qd_rigid( + fixed, moving, mxshift::VecLike, mxrot::Union{Number, VecLike}; + presmoothed = false, + SD = I, + minwidth_rot = default_minwidth_rot(CartesianIndices(fixed), SD), + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + initial_tfm = IdentityTransformation(), + print_interval = 100, + kwargs... + ) #TODO make sure isrotation(initial_tfm.linear) returns true, else throw a warning and an option to terminate? do I still allow the main block to run? if initial_tfm == IdentityTransformation() || isrotation(initial_tfm.linear) else @@ -198,8 +210,8 @@ function qd_rigid(fixed, moving, mxshift::VecLike, mxrot::Union{Number,VecLike}; end mxrot = [mxrot...] print_interval < typemax(Int) && print("Running coarse step\n") - tfm_coarse, mm_coarse = qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot; SD=SD, initial_tfm=initial_tfm, thresh=thresh, print_interval=print_interval, kwargs...) + tfm_coarse, mm_coarse = qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot; SD = SD, initial_tfm = initial_tfm, thresh = thresh, print_interval = print_interval, kwargs...) print_interval < typemax(Int) && print("Obtained mismatch $mm_coarse from coarse step, running fine step\n") - final_tfm, mm_fine = qd_rigid_fine(fixed, moving, mxrot./2, minwidth_rot; SD=SD, initial_tfm=tfm_coarse, thresh=thresh, print_interval=print_interval, kwargs...) + final_tfm, mm_fine = qd_rigid_fine(fixed, moving, mxrot ./ 2, minwidth_rot; SD = SD, initial_tfm = tfm_coarse, thresh = thresh, print_interval = print_interval, kwargs...) return final_tfm, mm_fine end diff --git a/src/translations.jl b/src/translations.jl index 219d3f7..cd5616c 100644 --- a/src/translations.jl +++ b/src/translations.jl @@ -1,27 +1,29 @@ #################### Translation Search ########################## -function tfmshift(params, img::AbstractArray{T,N}) where {T,N} +function tfmshift(params, img::AbstractArray{T, N}) where {T, N} length(params) == N || throw(DimensionMismatch("expected $N parameters, got $(length(params))")) return Translation(params...) end #slow because it warps for every shift instead of using fourier method -function translate_mm_slow(params, fixed, moving, thresh; initial_tfm=IdentityTransformation()) +function translate_mm_slow(params, fixed, moving, thresh; initial_tfm = IdentityTransformation()) tfm = initial_tfm ∘ tfmshift(params, moving) moving, fixed = warp_and_intersect(moving, fixed, tfm) - mm = mismatch_zeroshift(fixed, moving; normalization=:intensity) - return ratio(mm, thresh; fillval=Inf) + mm = mismatch_zeroshift(fixed, moving; normalization = :intensity) + return ratio(mm, thresh; fillval = Inf) end -function qd_translate_fine(fixed, moving; - initial_tfm=IdentityTransformation(), - minwidth=fill(0.01, ndims(fixed)), - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - kwargs...) - f(x) = translate_mm_slow(x, fixed, moving, thresh; initial_tfm=initial_tfm) +function qd_translate_fine( + fixed, moving; + initial_tfm = IdentityTransformation(), + minwidth = fill(0.01, ndims(fixed)), + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + kwargs... + ) + f(x) = translate_mm_slow(x, fixed, moving, thresh; initial_tfm = initial_tfm) upper = fill(1.0, ndims(fixed)) lower = -upper - root, x0 = _analyze(f, lower, upper; minwidth=minwidth, print_interval=100, kwargs...) + root, x0 = _analyze(f, lower, upper; minwidth = minwidth, print_interval = 100, kwargs...) box = minimum(root) tfmfine = initial_tfm ∘ tfmshift(position(box, x0), moving) return tfmfine, value(box) @@ -55,11 +57,13 @@ If the `crop` keyword arg is `true` then `fixed` is cropped by `mxshift` (after so that there will be complete overlap between `fixed` and `moving` for any evaluated shift. This avoids edge effects that can occur due to normalization when the transformed `moving` doesn't fully overlap with `fixed`. """ -function qd_translate(fixed, moving, mxshift; - presmoothed=false, - thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])), - initial_tfm=IdentityTransformation(), - minwidth=fill(0.01, ndims(fixed)), print_interval=100, crop=false, kwargs...) +function qd_translate( + fixed, moving, mxshift; + presmoothed = false, + thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])), + initial_tfm = IdentityTransformation(), + minwidth = fill(0.01, ndims(fixed)), print_interval = 100, crop = false, kwargs... + ) fixed, moving = float(fixed), float(moving) if presmoothed moving = qinterp(eltype(fixed), moving) @@ -68,7 +72,7 @@ function qd_translate(fixed, moving, mxshift; if crop #we enforce that moving is always bigger than fixed by amount 2*(maxshift+1) (the +1 is for the fine step) sz = size(fixed) .- (2 .* mxshift) - if any(size(moving) .< (2 .* (mxshift.+1))) + if any(size(moving) .< (2 .* (mxshift .+ 1))) error("Moving image size must be at least 2 * (mxshift+1) when crop_edges is set to true") end cropped_inds_f = crop_rng.(axes(fixed), mxshift) @@ -80,10 +84,10 @@ function qd_translate(fixed, moving, mxshift; fixed_inner = fixed moving_inner = moving_fine = moving end - best_shft, mm = best_shift(fixed_inner, moving_inner, mxshift, thresh; normalization=:intensity, initial_tfm=initial_tfm) + best_shft, mm = best_shift(fixed_inner, moving_inner, mxshift, thresh; normalization = :intensity, initial_tfm = initial_tfm) tfm_coarse = initial_tfm ∘ Translation(best_shft) print_interval < typemax(Int) && print("Running fine step\n") - return qd_translate_fine(fixed_inner, moving_fine; initial_tfm=tfm_coarse, thresh=thresh, minwidth=minwidth, print_interval=print_interval, kwargs...) + return qd_translate_fine(fixed_inner, moving_fine; initial_tfm = tfm_coarse, thresh = thresh, minwidth = minwidth, print_interval = print_interval, kwargs...) end -crop_rng(rng::AbstractUnitRange{Int}, amt::Int) = (first(rng)+amt):(last(rng)-amt) +crop_rng(rng::AbstractUnitRange{Int}, amt::Int) = (first(rng) + amt):(last(rng) - amt) diff --git a/src/util.jl b/src/util.jl index 5833615..58d4d30 100644 --- a/src/util.jl +++ b/src/util.jl @@ -24,12 +24,12 @@ expressed in pixels if `normalization=:pixels` and in sum-of-squared-intensity i One can compute an overall transformation by composing `initial_tfm` with the returned shift `I`. """ -function best_shift(fixed, moving, mxshift, thresh; normalization=:intensity, initial_tfm=IdentityTransformation()) +function best_shift(fixed, moving, mxshift, thresh; normalization = :intensity, initial_tfm = IdentityTransformation()) moving, fixed = warp_and_intersect(moving, fixed, initial_tfm) - mms = mismatch(fixed, moving, mxshift; normalization=normalization) + mms = mismatch(fixed, moving, mxshift; normalization = normalization) best_i = argmin_mismatch(mms, thresh) mm = mms[best_i] - return best_i.I, ratio(mm, thresh; fillval=typemax(eltype(mm))) + return best_i.I, ratio(mm, thresh; fillval = typemax(eltype(mm))) end """ @@ -66,15 +66,15 @@ See [`arrayscale`](@ref) for more information. """ pscale(t::Translation, SD::AbstractMatrix) = pscale(t, LinearMap(SD)) pscale(t::Translation, scale::LinearMap) = Translation(scale(t.translation)) -pscale(t::Translation, scale::UniformScaling) = Translation(scale.λ*t.translation) +pscale(t::Translation, scale::UniformScaling) = Translation(scale.λ * t.translation) #returns new minbounds and maxbounds with range sizes change by fac function scalebounds(minb, maxb, fac::Real) - orng = maxb.-minb - newradius = fac.*orng./2 - ctrs = minb.+(orng./2) - return ctrs.-newradius, ctrs.+newradius + orng = maxb .- minb + newradius = fac .* orng ./ 2 + ctrs = minb .+ (orng ./ 2) + return ctrs .- newradius, ctrs .+ newradius end """ @@ -90,12 +90,12 @@ along the diagonal and off the diagnonal, respectively. e.g. `dmax=0.05` implies that diagonal elements can range from 0.95 to 1.05. The space is centered on the identity matrix """ -function default_linmap_bounds(img::AbstractArray{T,N}; dmax=0.05, ndmax=0.05) where {T, N} - deltas = fill(abs(ndmax), N,N) - for i=1:N - deltas[i,i] = abs(dmax) +function default_linmap_bounds(img::AbstractArray{T, N}; dmax = 0.05, ndmax = 0.05) where {T, N} + deltas = fill(abs(ndmax), N, N) + for i in 1:N + deltas[i, i] = abs(dmax) end - return Matrix(1.0*I,N,N).-deltas, Matrix(1.0*I,N,N).+deltas + return Matrix(1.0 * I, N, N) .- deltas, Matrix(1.0 * I, N, N) .+ deltas end """ @@ -106,10 +106,10 @@ This can be useful for setting the `minwidth` parameter of QuadDIRECT when perfo full affine registration. `dmin` and `ndmin` set the tolerances for diagonal and off-diagonal elements of the linear transformation matrix, respectively. """ -function default_lin_minwidths(img::AbstractArray{T,N}; dmin=1e-5, ndmin=1e-5) where {T, N} - mat = fill(abs(ndmin), N,N) - for i=1:N - mat[i,i] = abs(dmin) +function default_lin_minwidths(img::AbstractArray{T, N}; dmin = 1.0e-5, ndmin = 1.0e-5) where {T, N} + mat = fill(abs(ndmin), N, N) + for i in 1:N + mat[i, i] = abs(dmin) end return mat[:] end @@ -122,28 +122,28 @@ Compute the rotation `θ` that results in largest change in coordinates of size `Δc` (in pixels) for any index in `ci`. When passed an array `img`, its `CartesianIndices` are used. """ -function default_minrot(ci::CartesianIndices, SD=I; Δc=0.1) +function default_minrot(ci::CartesianIndices, SD = I; Δc = 0.1) L = -Inf for x in CornerIterator(ci) - x′ = SD*SVector(Tuple(x)) # position of corner point in physical space + x′ = SD * SVector(Tuple(x)) # position of corner point in physical space L = max(L, norm(x′)) end - S2 = SD'*SD + S2 = SD' * SD if SD == I λ = 1 else F = eigen(S2) λ = minimum(F.values) end - ℓ = sqrt(λ)*Δc - return 2*asin(ℓ/(2*L)) + ℓ = sqrt(λ) * Δc + return 2 * asin(ℓ / (2 * L)) end -default_minrot(img::AbstractArray, SD=I; Δc=0.1) = default_minrot(CartesianIndices(img), SD; Δc) +default_minrot(img::AbstractArray, SD = I; Δc = 0.1) = default_minrot(CartesianIndices(img), SD; Δc) -default_minwidth_rot(ci::CartesianIndices{2}, SD=I; kwargs...) = +default_minwidth_rot(ci::CartesianIndices{2}, SD = I; kwargs...) = [default_minrot(ci, SD; kwargs...)] -function default_minwidth_rot(ci::CartesianIndices{3}, SD=I; kwargs...) +function default_minwidth_rot(ci::CartesianIndices{3}, SD = I; kwargs...) θ = default_minrot(ci, SD; kwargs...) return [θ, θ, θ] end @@ -173,15 +173,15 @@ julia> getSD(myimage) ``` """ getSD(A::AbstractArray) = getSD(spacedirections(A)) #takes advantage of how spacedirections automatically strips out the time component -function getSD(sd::NTuple{N,Tuple}) where N - SD = zeros(N,N) +function getSD(sd::NTuple{N, Tuple}) where {N} + SD = zeros(N, N) denom = oneunit(sd[1][1]) #dividing by a unit should remove unit inconsistancies - for i = 1:N - for j = 1:N - SD[i,j] = sd[i][j]/denom + for i in 1:N + for j in 1:N + SD[i, j] = sd[i][j] / denom end end - return SMatrix{N,N}(SD) #returns static matrix for faster type inferrence + return SMatrix{N, N}(SD) #returns static matrix for faster type inferrence end """ @@ -194,14 +194,14 @@ as an option. (Do not smooth `moving`.) `T` allows you to specify the output eltype (default `float(eltype(img))`). """ -function qsmooth(::Type{T}, img::AbstractArray{T2,N}) where {T,N,T2} - kern1 = centered(T[1/8, 3/4, 1/8]) # quadratic B-spline kernel - return imfilter(img, kernelfactors(ntuple(i->kern1, Val(N)))) +function qsmooth(::Type{T}, img::AbstractArray{T2, N}) where {T, N, T2} + kern1 = centered(T[1 / 8, 3 / 4, 1 / 8]) # quadratic B-spline kernel + return imfilter(img, kernelfactors(ntuple(i -> kern1, Val(N)))) end qsmooth(img::AbstractArray) = qsmooth(float(eltype(img)), img) -function qinterp(::Type{T}, moving) where T - widen1(ax) = first(ax)-1:last(ax)+1 +function qinterp(::Type{T}, moving) where {T} + widen1(ax) = (first(ax) - 1):(last(ax) + 1) # This sets up an interpolant object *without* using Interpolations.prefilter # Prefiltering is a form of "anti-smoothing," meaning that for quadratic interpolation @@ -241,7 +241,6 @@ end #sets splits based on lower and upper bounds function _analyze(f, lower, upper; kwargs...) - splits = ([[lower[i]; lower[i]+(upper[i]-lower[i])/2; upper[i]] for i=1:length(lower)]...,) - QuadDIRECT.analyze(f, splits, lower, upper; kwargs...) + splits = ([[lower[i]; lower[i] + (upper[i] - lower[i]) / 2; upper[i]] for i in 1:length(lower)]...,) + return QuadDIRECT.analyze(f, splits, lower, upper; kwargs...) end - diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 12c91fa..17e0f22 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -22,7 +22,7 @@ using Test #Helper to generate test image pairs function fixedmov(img, tfm) img = float(img) - img2 = warp(img,tfm) + img2 = warp(img, tfm) inds = OffsetArrays.IdentityUnitRange.(intersect.(axes(img), axes(img2))) fixed = img[inds...] moving = img2[inds...] @@ -40,9 +40,9 @@ function tfmtest(tfm, tfminv) diagtol = 0.005 offdiagtol = 0.005 vtol = 0.1 - @test all(x->(1-diagtol < x < 1+diagtol), diag(comp.linear)) - @test all(x->(-offdiagtol < x < offdiagtol), comp.linear.-Matrix(Diagonal(diag(comp.linear)))) - @test all(x-> x. (1 - diagtol < x < 1 + diagtol), diag(comp.linear)) + @test all(x -> (-offdiagtol < x < offdiagtol), comp.linear .- Matrix(Diagonal(diag(comp.linear)))) + return @test all(x -> x .< vtol, abs.(comp.translation)) end # helper for wrapping CuArrays in Offset Arrays @@ -52,78 +52,78 @@ end #rigid tests img = Float32.(testimage("cameraman")) SD = Matrix{Float64}(LinearAlgebra.I, 2, 2) -tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now +tfm = Translation(@SVector([14, 17])) ∘ LinearMap(RotMatrix(0.3)) #no distortion for now # fixed, moving = cu_wrap.(fixedmov(centered(img), tfm)); fixed, moving = fixedmov(centered(img), tfm) -mxshift = (100,100) #make sure this isn't too small +mxshift = (100, 100) #make sure this isn't too small mxrot = (0.5,) minwidth_rot = fill(0.002, 3) -tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002); +tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002); tfmtest(tfm, tform) -@test mm<0.001 +@test mm < 0.001 img3D = Float32.(testimage("mri-stack")) -SD3D = Matrix{Float64}(LinearAlgebra.I, 3,3 ) -tfm3D = Translation(@SVector([15, 10, 2]))∘LinearMap(RotXYZ(0.05, 0.02, 0)) #no distortion for now +SD3D = Matrix{Float64}(LinearAlgebra.I, 3, 3) +tfm3D = Translation(@SVector([15, 10, 2])) ∘ LinearMap(RotXYZ(0.05, 0.02, 0)) #no distortion for now fixed, moving = fixedmov(centered(img3D), tfm3D) mxshift = (30, 30, 5) #make sure this isn't too small mxrot = (0.1, 0.1, 0.1) minwidth_rot = fill(0.002, 3) -tform2, mm2 = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD3D, maxevals=1000, rtol=0, fvalue=0.0002); +tform2, mm2 = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD3D, maxevals = 1000, rtol = 0, fvalue = 0.0002); tfmtest(tfm3D, tform2) -@test mm2<0.01 +@test mm2 < 0.01 #Coppied over from qd_random @testset "QuadDIRECT tests with standard images" begin - img = testimage("cameraman"); + img = testimage("cameraman") #Translation (subpixel) tfm = Translation(@SVector([14.3, 17.6])) fixed, moving = fixedmov(img, tfm) - mxshift = (100,100) #make sure this isn't too small - tform, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, rtol=0, fvalue=0.0003) + mxshift = (100, 100) #make sure this isn't too small + tform, mm = qd_translate(fixed, moving, mxshift; maxevals = 1000, rtol = 0, fvalue = 0.0003) tfmtest(tfm, tform) #Rigid transform SD = Matrix{Float64}(LinearAlgebra.I, 2, 2) - tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now + tfm = Translation(@SVector([14, 17])) ∘ LinearMap(RotMatrix(0.3)) #no distortion for now fixed, moving = fixedmov(centered(img), tfm) - mxshift = (100,100) #make sure this isn't too small + mxshift = (100, 100) #make sure this isn't too small mxrot = (0.5,) minwidth_rot = fill(0.002, 3) - tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002) + tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002) tfmtest(tfm, tform) #with anisotropic sampling SD = Matrix(Diagonal([0.5; 1.0])) - tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.3)*SD) + tfm = Translation(@SVector([14.3, 17.8])) ∘ LinearMap(SD \ RotMatrix(0.3) * SD) fixed, moving = fixedmov(centered(img), tfm) - tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002) + tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002) tfmtest(tfm, arrayscale(tform, SD)) #Affine transform - tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.01)) + tfm = Translation(@SVector([14, 17])) ∘ LinearMap(RotMatrix(0.01)) #make it harder with nonuniform scaling scale = @SMatrix [1.005 0; 0 0.995] SD = Matrix{Float64}(LinearAlgebra.I, 2, 2) - tfm = AffineMap(tfm.linear*scale, tfm.translation) - mxshift = (100,100) #make sure this isn't too small + tfm = AffineMap(tfm.linear * scale, tfm.translation) + mxshift = (100, 100) #make sure this isn't too small fixed, moving = fixedmov(centered(img), tfm) - tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=1000, rtol=0, fvalue=0.0002) + tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002) tfmtest(tfm, tform) #with anisotropic sampling SD = Matrix(Diagonal([0.5; 1.0])) - tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(RotMatrix(0.1)) #Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.01)*SD) + tfm = Translation(@SVector([14.3, 17.8])) ∘ LinearMap(RotMatrix(0.1)) #Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.01)*SD) scale = @SMatrix [1.005 0; 0 0.995] - tfm = AffineMap(tfm.linear*scale, tfm.translation) + tfm = AffineMap(tfm.linear * scale, tfm.translation) tfm = arrayscale(tfm, SD) fixed, moving = fixedmov(centered(img), tfm) - tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=10000, rtol=0, fvalue=0.0002, ndmax = 0.25) + tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals = 10000, rtol = 0, fvalue = 0.0002, ndmax = 0.25) tform2 = arrayscale(tform, SD) tfmtest(tfm, tform2) end #tests with standard images @@ -131,33 +131,33 @@ end #tests with standard images @testset "Quadratic interpolation (issue #7)" begin samplefrom(n) = rand(Poisson(n)) - img = restrict(restrict(testimage("cameraman")))[2:end-1,2:end-1] + img = restrict(restrict(testimage("cameraman")))[2:(end - 1), 2:(end - 1)] # Convert to "photons" so we can mimic shot noise np = 1000 # maximum number of photons per pixel - img = round.(Int, np.*gray.(img)) - fixed = samplefrom.(img) + img = round.(Int, np .* gray.(img)) + fixed = samplefrom.(img) moving = samplefrom.(img) ff = qsmooth(fixed) - tform, mm = qd_translate(fixed, moving, (5, 5); print_interval=typemax(Int)) - tformq, mmq = qd_translate(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int)) + tform, mm = qd_translate(fixed, moving, (5, 5); print_interval = typemax(Int)) + tformq, mmq = qd_translate(ff, moving, (5, 5); presmoothed = true, print_interval = typemax(Int)) @test all(abs.(tformq.translation) .< abs.(tform.translation)) - tform, mm = qd_rigid(fixed, moving, (5, 5), 0.1; print_interval=typemax(Int)) - tformq, mmq = qd_rigid(ff, moving, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int)) - @test norm(tformq.linear-I) < norm(tform.linear-I) + tform, mm = qd_rigid(fixed, moving, (5, 5), 0.1; print_interval = typemax(Int)) + tformq, mmq = qd_rigid(ff, moving, (5, 5), 0.1; presmoothed = true, print_interval = typemax(Int)) + @test norm(tformq.linear - I) < norm(tform.linear - I) @test norm(tformq.translation) < norm(tform.translation) - tform, mm = qd_affine(fixed, moving, (5, 5); print_interval=typemax(Int)) - tformq, mmq = qd_affine(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int)) + tform, mm = qd_affine(fixed, moving, (5, 5); print_interval = typemax(Int)) + tformq, mmq = qd_affine(ff, moving, (5, 5); presmoothed = true, print_interval = typemax(Int)) @test mmq < mm # Test that we exactly reconstruct `qsmooth` with `presmoothed=true` - tformq, mmq = qd_translate(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int)) + tformq, mmq = qd_translate(ff, fixed, (5, 5); presmoothed = true, print_interval = typemax(Int)) @test all(iszero, tformq.translation) - @test mmq < 1e-10 - tformq, mmq = qd_rigid(ff, fixed, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int)) - @test mmq < 1e-8 - tformq, mmq = qd_affine(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int)) - @test mmq < 1e-6 # on 32-bit systems this can't be 1e-8, not quite sure why -end \ No newline at end of file + @test mmq < 1.0e-10 + tformq, mmq = qd_rigid(ff, fixed, (5, 5), 0.1; presmoothed = true, print_interval = typemax(Int)) + @test mmq < 1.0e-8 + tformq, mmq = qd_affine(ff, fixed, (5, 5); presmoothed = true, print_interval = typemax(Int)) + @test mmq < 1.0e-6 # on 32-bit systems this can't be 1e-8, not quite sure why +end diff --git a/test/gridsearch.jl b/test/gridsearch.jl index 83dcb29..8ac01f0 100644 --- a/test/gridsearch.jl +++ b/test/gridsearch.jl @@ -7,23 +7,23 @@ using Test @testset "Grid search rigid registration" begin ## 2D #note: if a is much smaller than this then it won't find the correct answer due to the mismatch normalization - a = rand(30,30) - b = transform(a, tformtranslate([2.0;0.0]) ∘ tformrotate(pi/6)) - tfm0 = tformtranslate([-2.0;0.0]) ∘ tformrotate(-pi/6) + a = rand(30, 30) + b = 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 # SD is now a keyword argument (was positional in v0.x) SD2 = Matrix{Float64}(I, 2, 2) - tfm, mm = RegisterQD.rotation_gridsearch(a, b, (11,11), [pi/6], [11]; SD=SD2) + tfm, mm = RegisterQD.rotation_gridsearch(a, b, (11, 11), [pi / 6], [11]; SD = SD2) @test tfm.translation == tfm0.translation @test tfm.linear == tfm0.linear ## 3D #note: if a is much smaller than this then it won't find the correct answer due to the mismatch normalization - a = rand(30,30,30) - b = transform(a, tformtranslate([2.0;0.0;0.0]) ∘ tformrotate([1.0;0;0], pi/4)) - tfm0 = tformtranslate([-2.0;0.0;0.0]) ∘ tformrotate([1.0;0;0], -pi/4) + a = rand(30, 30, 30) + b = transform(a, tformtranslate([2.0;0.0;0.0]) ∘ tformrotate([1.0;0;0], pi / 4)) + tfm0 = tformtranslate([-2.0;0.0;0.0]) ∘ tformrotate([1.0;0;0], -pi / 4) #note: maxshift must be GREATER than the true shift in order to find the true shift - tfm, mm = RegisterQD.rotation_gridsearch(a, b, (3,3,3), [pi/4, pi/4, pi/4], [5;5;5]) + tfm, mm = RegisterQD.rotation_gridsearch(a, b, (3, 3, 3), [pi / 4, pi / 4, pi / 4], [5;5;5]) @test tfm.translation == tfm0.translation @test tfm.linear == tfm0.linear end diff --git a/test/initial_tfm.jl b/test/initial_tfm.jl index 1e46160..54e9601 100644 --- a/test/initial_tfm.jl +++ b/test/initial_tfm.jl @@ -11,200 +11,200 @@ using Test g = 0.2:0.2:1.2 gradcube = g .* reshape(g, 1, 6) .* reshape(g, 1, 1, 6) -gradcube[1,4,2] = 0 -gradcube[4,5,2] = 0 #break up any potential rotational symmetry +gradcube[1, 4, 2] = 0 +gradcube[4, 5, 2] = 0 #break up any potential rotational symmetry -mxshift = (1,1,1); -mxrot = (0.01,0.01,0.01); +mxshift = (1, 1, 1); +mxrot = (0.01, 0.01, 0.01); -EYE =Matrix(1.0*I, 3,3) +EYE = Matrix(1.0 * I, 3, 3) @testset "initial_tfm improves rigid rotational alignment" begin - testimage1 = zeros(Float64,10,10,10) + testimage1 = zeros(Float64, 10, 10, 10) testimage1 .= NaN - testimage1[3:8,3:8,3:8] = gradcube + testimage1[3:8, 3:8, 3:8] = gradcube testimage1 = centered(testimage1) - tform = RotXYZ(0.1,0.1,0.1) - mytform = AffineMap(tform, [0,0,0]) + tform = RotXYZ(0.1, 0.1, 0.1) + mytform = AffineMap(tform, [0, 0, 0]) testimage2 = warp(testimage1, mytform, axes(testimage1)) # insufficient parameters + no initiat_tfm causes large misalignment - tformtest1, mm1 = qd_rigid(testimage2, testimage1, mxshift, mxrot; print_interval=typemax(Int)) # - @test !(mm1 < 1e-8) + tformtest1, mm1 = qd_rigid(testimage2, testimage1, mxshift, mxrot; print_interval = typemax(Int)) # + @test !(mm1 < 1.0e-8) @test !isapprox(tformtest1, mytform, atol = 0.1) #initial_tfm improves alignment - tformtest2, mm2 = qd_rigid(testimage2, testimage1, mxshift, mxrot; print_interval=typemax(Int), initial_tfm = mytform) # - #with the initial_tfm being the true tfm, this should give back the true rotation + tformtest2, mm2 = qd_rigid(testimage2, testimage1, mxshift, mxrot; print_interval = typemax(Int), initial_tfm = mytform) # + #with the initial_tfm being the true tfm, this should give back the true rotation - @test mm2 < 1e-8 + @test mm2 < 1.0e-8 @test isapprox(tformtest2, mytform, atol = 0.0001) @test isrotation(tformtest2.linear) end #Test initial_tfm improves rotational alignment for rigid @testset "initial_tfm improves rigid translational alignment" begin - testimage3 = zeros(10,20,10) + testimage3 = zeros(10, 20, 10) testimage3 .= NaN - testimage3[3:8,3:8,3:8] = gradcube + testimage3[3:8, 3:8, 3:8] = gradcube testimage3 = centered(testimage3) - testimage4 = zeros(10,20,10) + testimage4 = zeros(10, 20, 10) testimage4 .= NaN - testimage4[3:8,13:18,3:8] = gradcube + testimage4[3:8, 13:18, 3:8] = gradcube testimage4 = centered(testimage4) - init_tfm = AffineMap(diagm(0 => fill(1.0,3)), [0.0,10.0,0.0]) + init_tfm = AffineMap(diagm(0 => fill(1.0, 3)), [0.0, 10.0, 0.0]) tformtest3 = nothing mm3 = nothing #no initial_tfm should fail due to total lack of overlap @test try - tformtest3, mm3 = qd_rigid(testimage3, testimage4, mxshift, mxrot; print_interval=typemax(Int)) # + tformtest3, mm3 = qd_rigid(testimage3, testimage4, mxshift, mxrot; print_interval = typemax(Int)) # false catch err true # this should break end - tformtest4, mm4 = qd_rigid(testimage3, testimage4, mxshift, mxrot; print_interval=typemax(Int), initial_tfm = init_tfm) # + tformtest4, mm4 = qd_rigid(testimage3, testimage4, mxshift, mxrot; print_interval = typemax(Int), initial_tfm = init_tfm) # - @test mm4 < 1e-8 #this mismatch should be lower! + @test mm4 < 1.0e-8 #this mismatch should be lower! @test isapprox(tformtest4.translation, [0, 10, 0]) @test isapprox(tformtest4.linear, EYE) end @testset "Test initial_tfm creates real rotations." begin #some of this seems redundant - testimage1 = zeros(Float64,10,10,10) + testimage1 = zeros(Float64, 10, 10, 10) testimage1 .= NaN - testimage1[3:8,3:8,3:8] = gradcube + testimage1[3:8, 3:8, 3:8] = gradcube testimage1 = centered(testimage1) - tform = RotXYZ(0.1,0.1,0.1) - mytform = AffineMap(tform, [0,0,0]) + tform = RotXYZ(0.1, 0.1, 0.1) + mytform = AffineMap(tform, [0, 0, 0]) testimage2 = warp(testimage1, mytform, axes(testimage1)) - mxrot2 = (0.2,0.2,0.2); + mxrot2 = (0.2, 0.2, 0.2) minwidth_rot = RegisterQD.default_minwidth_rot(CartesianIndices(testimage2), EYE) # tests with equal spaces produces real rotations # TODO: the next test is machine-dependent, figure out why - tformtest0, mm0 = qd_rigid(testimage2, testimage1, mxshift, mxrot2; fvalue=1e-5, rtol=0, print_interval=typemax(Int)) + tformtest0, mm0 = qd_rigid(testimage2, testimage1, mxshift, mxrot2; fvalue = 1.0e-5, rtol = 0, print_interval = typemax(Int)) - @test mm0 < 1e-5 + @test mm0 < 1.0e-5 @test isapprox(tformtest0, mytform, atol = 0.1) @test isrotation(tformtest0.linear) - tformtest01, mm01 = RegisterQD.qd_rigid_coarse(testimage2, testimage1, mxshift, [mxrot2...], minwidth_rot; SD=EYE, print_interval=typemax(Int)) + tformtest01, mm01 = RegisterQD.qd_rigid_coarse(testimage2, testimage1, mxshift, [mxrot2...], minwidth_rot; SD = EYE, print_interval = typemax(Int)) @test isrotation(tformtest01.linear) tformtest02 = nothing mm02 = nothing - tformtest02, mm02 = RegisterQD.qd_rigid_fine(testimage2, testimage1, [mxrot2...]./2, minwidth_rot; SD=EYE, print_interval=typemax(Int)) + tformtest02, mm02 = RegisterQD.qd_rigid_fine(testimage2, testimage1, [mxrot2...] ./ 2, minwidth_rot; SD = EYE, print_interval = typemax(Int)) @test isrotation(tformtest02.linear) #with skewed spacing produces real rotations. - testimage5 = testimage1[:,:,-4:2:5] - testimage6 = testimage2[:,:,-4:2:5] + testimage5 = testimage1[:, :, -4:2:5] + testimage6 = testimage2[:, :, -4:2:5] ps = (1, 1, 2) - SD =SArray{Tuple{3,3}}(Diagonal(SVector(ps) ./ minimum(ps))) + SD = SArray{Tuple{3, 3}}(Diagonal(SVector(ps) ./ minimum(ps))) - tformtest5, mm5 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int)) - @test mm5 <1e-4 + tformtest5, mm5 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD = SD, print_interval = typemax(Int)) + @test mm5 < 1.0e-4 @test isapprox(tformtest5, mytform, atol = 1) @test isrotation(tformtest5.linear) - @test !isrotation(SD*tformtest5.linear*inv(SD)) + @test !isrotation(SD * tformtest5.linear * inv(SD)) - tformtest55, mm55 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int), initial_tfm = tformtest5) - @test mm55 <1e-4 + tformtest55, mm55 = qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD = SD, print_interval = typemax(Int), initial_tfm = tformtest5) + @test mm55 < 1.0e-4 @test isapprox(tformtest55, mytform, atol = 1) @test isrotation(tformtest55.linear) - @test !isrotation(SD*tformtest55.linear*inv(SD)) + @test !isrotation(SD * tformtest55.linear * inv(SD)) #coarse and fine produce real rotations - tformtest6, mm6 = RegisterQD.qd_rigid_coarse(testimage6, testimage5, mxshift, mxrot2, minwidth_rot; SD=SD, print_interval=typemax(Int)) - @test mm6 <1e-4 + tformtest6, mm6 = RegisterQD.qd_rigid_coarse(testimage6, testimage5, mxshift, mxrot2, minwidth_rot; SD = SD, print_interval = typemax(Int)) + @test mm6 < 1.0e-4 @test isrotation(tformtest6.linear) - @test !isrotation(SD*tformtest6.linear*inv(SD)) + @test !isrotation(SD * tformtest6.linear * inv(SD)) - tformtest66, mm66 = RegisterQD.qd_rigid_coarse(testimage6, testimage5, mxshift, mxrot2, minwidth_rot; SD=SD, print_interval=typemax(Int), initial_tfm = tformtest6) - @test mm66 <1e-4 + tformtest66, mm66 = RegisterQD.qd_rigid_coarse(testimage6, testimage5, mxshift, mxrot2, minwidth_rot; SD = SD, print_interval = typemax(Int), initial_tfm = tformtest6) + @test mm66 < 1.0e-4 @test isrotation(tformtest66.linear) - @test !isrotation(SD*tformtest66.linear*inv(SD)) + @test !isrotation(SD * tformtest66.linear * inv(SD)) - tformtest7, mm7 = RegisterQD.qd_rigid_fine(testimage6, testimage5, [mxrot2...]./2, minwidth_rot; SD=SD, print_interval=typemax(Int)) - @test mm7 <1e-4 + tformtest7, mm7 = RegisterQD.qd_rigid_fine(testimage6, testimage5, [mxrot2...] ./ 2, minwidth_rot; SD = SD, print_interval = typemax(Int)) + @test mm7 < 1.0e-4 @test isrotation(tformtest7.linear) - @test !isrotation(SD*tformtest7.linear*inv(SD)) + @test !isrotation(SD * tformtest7.linear * inv(SD)) - tformtest77, mm77 = RegisterQD.qd_rigid_fine(testimage6, testimage5, [mxrot2...]./2, minwidth_rot; SD=SD, print_interval=typemax(Int), initial_tfm = tformtest7) - @test mm77 <1e-4 + tformtest77, mm77 = RegisterQD.qd_rigid_fine(testimage6, testimage5, [mxrot2...] ./ 2, minwidth_rot; SD = SD, print_interval = typemax(Int), initial_tfm = tformtest7) + @test mm77 < 1.0e-4 @test isrotation(tformtest77.linear) - @test !isrotation(SD*tformtest77.linear*inv(SD)) + @test !isrotation(SD * tformtest77.linear * inv(SD)) #fails due to specific error in qd_rigid_fine # a non rigid initial_tfm does not return a rigid transformation and prints a warning tformtest8, mm8 = @test_logs (:warn, "initial_tfm is not a rigid transformation") begin - qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD=SD, print_interval=typemax(Int), initial_tfm = AffineMap(SD\tformtest5.linear*SD, SD\tformtest5.translation)) + qd_rigid(testimage6, testimage5, mxshift, mxrot2; SD = SD, print_interval = typemax(Int), initial_tfm = AffineMap(SD \ tformtest5.linear * SD, SD \ tformtest5.translation)) end - @test mm8 <1e-4 + @test mm8 < 1.0e-4 @test isapprox(tformtest8, mytform, atol = 1) @test !isrotation(tformtest8.linear) end @testset "Test initial_tfm improves translational alignment for affine" begin - testimage3 = zeros(10,20,10) + testimage3 = zeros(10, 20, 10) testimage3 .= NaN - testimage3[3:8,3:8,3:8] = gradcube + testimage3[3:8, 3:8, 3:8] = gradcube testimage3 = centered(testimage3) - testimage4 = zeros(10,20,10) + testimage4 = zeros(10, 20, 10) testimage4 .= NaN - testimage4[3:8,13:18,3:8] = gradcube + testimage4[3:8, 13:18, 3:8] = gradcube testimage4 = centered(testimage4) - init_tfm = AffineMap(diagm(0 => fill(1.0,3)), [0.0,10.0,0.0]) + init_tfm = AffineMap(diagm(0 => fill(1.0, 3)), [0.0, 10.0, 0.0]) tformtest3 = nothing mm3 = nothing @test try - tformtest3, mm3 = qd_affine(testimage3, testimage4, mxshift; print_interval=typemax(Int)) # + tformtest3, mm3 = qd_affine(testimage3, testimage4, mxshift; print_interval = typemax(Int)) # false catch err true # this should break end - tformtest4, mm4 = qd_affine(testimage3, testimage4, mxshift; print_interval = typemax(Int), initial_tfm = init_tfm, fvalue = 1e-5) # - @test mm4 < 1e-4#this mismatch should be lower! + tformtest4, mm4 = qd_affine(testimage3, testimage4, mxshift; print_interval = typemax(Int), initial_tfm = init_tfm, fvalue = 1.0e-5) # + @test mm4 < 1.0e-4 #this mismatch should be lower! @test isapprox(tformtest4.translation, [0.0, 10.0, 0.0]) end @testset "Test initial_tfm improves general alignment for affine" begin - testimage1 = zeros(Float64,10,10,10) + testimage1 = zeros(Float64, 10, 10, 10) testimage1 .= NaN - testimage1[3:8,3:8,3:8] = gradcube + testimage1[3:8, 3:8, 3:8] = gradcube testimage1 = centered(testimage1) - tform = RotXYZ(0.1,0.1,0.1) - mytform = AffineMap(tform, [0.0,0.0,0.0]) + tform = RotXYZ(0.1, 0.1, 0.1) + mytform = AffineMap(tform, [0.0, 0.0, 0.0]) testimage2 = warp(testimage1, mytform, axes(testimage1)) - tformtest1, mm1 = qd_affine(testimage2, testimage1, mxshift; print_interval=typemax(Int)) # + tformtest1, mm1 = qd_affine(testimage2, testimage1, mxshift; print_interval = typemax(Int)) # #as the max-rotation is set too low, this should give a bad mismatch, but should still work - @test !(mm1 < 1e-8) + @test !(mm1 < 1.0e-8) @test !isapprox(tformtest1, mytform, atol = 0.1) - tformtest2, mm2 = qd_affine(testimage2, testimage1, mxshift; print_interval=typemax(Int), initial_tfm = mytform, fvalue = 1e-5) # + tformtest2, mm2 = qd_affine(testimage2, testimage1, mxshift; print_interval = typemax(Int), initial_tfm = mytform, fvalue = 1.0e-5) # #with the initial_tfm being the true tfm, this should give back the true rotation - @test mm2 < 1e-8 + @test mm2 < 1.0e-8 @test isapprox(tformtest2, mytform, atol = 0.0001) #this may not work so well. see qd_standard. end #Test initial_tfm improves rotational alignment for rigid diff --git a/test/qd_random.jl b/test/qd_random.jl index 150fbfa..afb951c 100644 --- a/test/qd_random.jl +++ b/test/qd_random.jl @@ -18,153 +18,153 @@ using Test, TestImages ##### Translations #2D - moving = rand(50,50) + moving = rand(50, 50) tfm0 = Translation(-4.7, 5.1) #ground truth newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) fixed = etp(Base.axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) - mxshift = (10,10) + mxshift = (10, 10) - tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, rtol=0) + tfm, mm = qd_translate(fixed, moving, mxshift; maxevals = 1000, thresh = thresh, rtol = 0) - @test sum(abs.(tfm0.translation - tfm.translation)) < 1e-3 + @test sum(abs.(tfm0.translation - tfm.translation)) < 1.0e-3 #again with cropping - tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, crop=true, rtol=0) - @test sum(abs.(tfm0.translation - tfm.translation)) < 1e-3 + tfm, mm = qd_translate(fixed, moving, mxshift; maxevals = 1000, thresh = thresh, crop = true, rtol = 0) + @test sum(abs.(tfm0.translation - tfm.translation)) < 1.0e-3 #3D - moving = rand(30,30,30) - tfm0 = Translation(-0.9, 2.1,1.2) #ground truth + moving = rand(30, 30, 30) + tfm0 = Translation(-0.9, 2.1, 1.2) #ground truth newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) - mxshift = (5,5,5) + mxshift = (5, 5, 5) - tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, rtol=0) + tfm, mm = qd_translate(fixed, moving, mxshift; maxevals = 1000, thresh = thresh, rtol = 0) @test sum(abs.(tfm0.translation - tfm.translation)) < 0.1 ######Rotations + Translations #2D - moving = centered(rand(50,50)) - tfm0 = Translation(-4.0, 5.0) ∘ LinearMap(RotMatrix(pi/360)) #ground truth + moving = centered(rand(50, 50)) + tfm0 = Translation(-4.0, 5.0) ∘ LinearMap(RotMatrix(pi / 360)) #ground truth newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) - mxshift = (10,10) - mxrot = pi/90 + mxshift = (10, 10) + mxrot = pi / 90 minwidth_rot = [0.0002] SD = SDiagonal(@SVector(ones(ndims(fixed)))) - tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, thresh=thresh, maxevals=1000, rtol=0, fvalue=1e-8) + tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, thresh = thresh, maxevals = 1000, rtol = 0, fvalue = 1.0e-8) tfm = arrayscale(tfm, SD) - @test sum(abs.(tfm0.linear - tfm.linear)) < 1e-3 + @test sum(abs.(tfm0.linear - tfm.linear)) < 1.0e-3 #3D - moving = centered(rand(30,30,30)) - tfm0 = Translation(-1.0, 2.1,1.2) ∘ LinearMap(RotXYZ(pi/360, pi/180, pi/220)) #ground truth + moving = centered(rand(30, 30, 30)) + tfm0 = Translation(-1.0, 2.1, 1.2) ∘ LinearMap(RotXYZ(pi / 360, pi / 180, pi / 220)) #ground truth newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.1 * sum(abs2.(fixed[.!(isnan.(fixed))])) - mxshift = (5,5,5) - mxrot = [pi/90; pi/90; pi/90] + mxshift = (5, 5, 5) + mxrot = [pi / 90; pi / 90; pi / 90] minwidth_rot = fill(0.0002, 3) SD = SDiagonal(@SVector(ones(ndims(fixed)))) - tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, thresh=thresh, maxevals=Sys.iswindows() ? 2000 : 1000, rtol=0) + tfm, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, thresh = thresh, maxevals = Sys.iswindows() ? 2000 : 1000, rtol = 0) tfm = arrayscale(tfm, SD) @test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(RotXYZ(tfm.linear)[:], tfm.translation))) < 0.1 -#NOTE: the 2D test below fails rarely and the 3D test fails often, apparently because full affine is too difficult with these images + #NOTE: the 2D test below fails rarely and the 3D test fails often, apparently because full affine is too difficult with these images #####General Affine Transformations #2D - moving = centered(rand(50,50)) - shft = SArray{Tuple{2}}(rand(2).+2.0) + moving = centered(rand(50, 50)) + shft = SArray{Tuple{2}}(rand(2) .+ 2.0) #random displacement from the identity matrix - mat = SArray{Tuple{2,2}}([1 0; 0 1] + rand(2,2)./40 + -rand(2,2)./40) + mat = SArray{Tuple{2, 2}}([1 0; 0 1] + rand(2, 2) ./ 40 + -rand(2, 2) ./ 40) tfm0 = AffineMap(mat, shft) #ground truth newfixed = warp(moving, tfm0) itp = interpolate(newfixed, BSpline(Linear())) etp = extrapolate(itp, NaN) fixed = etp(axes(moving)...) #often the warped array has one-too-many pixels in one or more dimensions due to extrapolation thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))])) - mxshift = (5,5) + mxshift = (5, 5) SD = SDiagonal(@SVector(ones(ndims(fixed)))) - tfm, mm = qd_affine(fixed, moving, mxshift; SD=SD, thresh=thresh, rtol=0, fvalue=1e-3, maxevals=10^4) + tfm, mm = qd_affine(fixed, moving, mxshift; SD = SD, thresh = thresh, rtol = 0, fvalue = 1.0e-3, maxevals = 10^4) tfm = arrayscale(tfm, SD) @test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1 - #The tests below fail. Probably two factors contributing to failure: - # 1) Full affine 3D is a lot of parameters so it's just difficuilt (12 parameters) - # 2) The current way of computing mismatch may be flawed for scaling transformations - # because the algorithm output shows it tends to compress/expand the image in order to remove - # pixels from the denominator in the mismatch calculation (especially bad for small images). - #3D - #moving = centered(rand(20,20,20)); - #shft = SArray{Tuple{3}}(2.6, 0.1, -3.3); - ##random displacement from the identity matrix - #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./30 + -rand(3,3)./30); - #tfm0 = AffineMap(mat, shft); #ground truth - #newfixed = warp(moving, tfm0); - #inds = intersect.(axes(moving), axes(newfixed)) - #fixed = newfixed[inds...] - #moving = moving[inds...] - #thresh = 0.1 * (sum(abs2.(fixed[.!(isnan.(fixed))]))+sum(abs2.(moving[.!(isnan.(moving))]))); - #mxshift = (10,10,10) - #SD = eye(ndims(fixed)); - - #tfm, mm = qd_affine(centered(fixed), centered(moving), mxshift, SD; thresh=thresh, rtol=0, fvalue=1e-4); - - #@test mm <= 1e-4 - #@test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1 - - #not random - #moving = zeros(10,10,10); - #moving[5:7, 5:7, 5:7] = 1.0 - #moving = centered(moving) - #shft = SArray{Tuple{3}}(0.6, 0.1, -0.3); - ##random displacement from the identity matrix - #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./50 + -rand(3,3)./50); - #tfm00 = AffineMap(mat, shft); - ##tfm0 = recenter(tfm00, center(moving)); #ground truth - #tfm0 = tfm00 #ground truth - #newfixed = warp(moving, tfm0); - #inds = intersect.(axes(moving), axes(newfixed)) - #fixed = newfixed[inds...] - #moving = moving[inds...] - #thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))])); - #mxshift = (5,5,5) - #SD = eye(ndims(fixed)); - #@test RegisterOptimize.aff(vcat(tfm00.translation[:], tfm00.linear[:]), fixed, SD) == tfm0 - - #tfm, mm = qd_affine(centered(fixed), centered(moving), mxshift, SD; thresh=thresh, rtol=0, fvalue=1e-4); - - #@test mm <= 1e-4 - #@test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1 + #The tests below fail. Probably two factors contributing to failure: + # 1) Full affine 3D is a lot of parameters so it's just difficuilt (12 parameters) + # 2) The current way of computing mismatch may be flawed for scaling transformations + # because the algorithm output shows it tends to compress/expand the image in order to remove + # pixels from the denominator in the mismatch calculation (especially bad for small images). + #3D + #moving = centered(rand(20,20,20)); + #shft = SArray{Tuple{3}}(2.6, 0.1, -3.3); + ##random displacement from the identity matrix + #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./30 + -rand(3,3)./30); + #tfm0 = AffineMap(mat, shft); #ground truth + #newfixed = warp(moving, tfm0); + #inds = intersect.(axes(moving), axes(newfixed)) + #fixed = newfixed[inds...] + #moving = moving[inds...] + #thresh = 0.1 * (sum(abs2.(fixed[.!(isnan.(fixed))]))+sum(abs2.(moving[.!(isnan.(moving))]))); + #mxshift = (10,10,10) + #SD = eye(ndims(fixed)); + + #tfm, mm = qd_affine(centered(fixed), centered(moving), mxshift, SD; thresh=thresh, rtol=0, fvalue=1e-4); + + #@test mm <= 1e-4 + #@test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1 + + #not random + #moving = zeros(10,10,10); + #moving[5:7, 5:7, 5:7] = 1.0 + #moving = centered(moving) + #shft = SArray{Tuple{3}}(0.6, 0.1, -0.3); + ##random displacement from the identity matrix + #mat = SArray{Tuple{3,3}}(eye(3) + rand(3,3)./50 + -rand(3,3)./50); + #tfm00 = AffineMap(mat, shft); + ##tfm0 = recenter(tfm00, center(moving)); #ground truth + #tfm0 = tfm00 #ground truth + #newfixed = warp(moving, tfm0); + #inds = intersect.(axes(moving), axes(newfixed)) + #fixed = newfixed[inds...] + #moving = moving[inds...] + #thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))])); + #mxshift = (5,5,5) + #SD = eye(ndims(fixed)); + #@test RegisterOptimize.aff(vcat(tfm00.translation[:], tfm00.linear[:]), fixed, SD) == tfm0 + + #tfm, mm = qd_affine(centered(fixed), centered(moving), mxshift, SD; thresh=thresh, rtol=0, fvalue=1e-4); + + #@test mm <= 1e-4 + #@test sum(abs.(vcat(tfm0.linear[:], tfm0.translation) - vcat(tfm.linear[:], tfm.translation))) < 0.1 #TODO this test passes if run individually, but breaks when included in the testset. end #tests with random images @testset "Suppression of printing" begin a, b = rand(5, 5), rand(5, 5) ca, cb = centered(a), centered(b) - mxshift = (2,2) + mxshift = (2, 2) mxrot = 0.5 # minwidth_rot = 1e-3 mktemp() do path, io redirect_stdout(io) do - qd_translate(a, b, (2,2); print_interval=typemax(Int)) - qd_rigid(ca, cb, mxshift, mxrot; print_interval=typemax(Int)) - qd_affine(ca, cb, mxshift; print_interval=typemax(Int)) + qd_translate(a, b, (2, 2); print_interval = typemax(Int)) + qd_rigid(ca, cb, mxshift, mxrot; print_interval = typemax(Int)) + qd_affine(ca, cb, mxshift; print_interval = typemax(Int)) end flush(io) str = read(path, String) diff --git a/test/qd_standard.jl b/test/qd_standard.jl index 6f87a58..1e2800b 100644 --- a/test/qd_standard.jl +++ b/test/qd_standard.jl @@ -17,7 +17,7 @@ using Random #Helper to generate test image pairs function fixedmov(img, tfm) img = float(img) - img2 = warp(img,tfm) + img2 = warp(img, tfm) inds = OffsetArrays.IdentityUnitRange.(intersect.(axes(img), axes(img2))) fixed = img[inds...] moving = img2[inds...] @@ -35,9 +35,9 @@ function tfmtest(tfm, tfminv) diagtol = 0.005 offdiagtol = 0.005 vtol = 0.1 - @test all(x->(1-diagtol < x < 1+diagtol), diag(comp.linear)) - @test all(x->(-offdiagtol < x < offdiagtol), comp.linear.-Matrix(Diagonal(diag(comp.linear)))) - @test all(abs.(comp.translation) .< vtol) + @test all(x -> (1 - diagtol < x < 1 + diagtol), diag(comp.linear)) + @test all(x -> (-offdiagtol < x < offdiagtol), comp.linear .- Matrix(Diagonal(diag(comp.linear)))) + return @test all(abs.(comp.translation) .< vtol) end # tests with standard images @@ -45,50 +45,50 @@ end # answer to is the inverse of an input transformation. This seems to catch # a different set of errors than the tests above) @testset "QuadDIRECT tests with standard images" begin - img = testimage("cameraman"); + img = testimage("cameraman") #Translation (subpixel) tfm = Translation(@SVector([14.3, 17.6])) fixed, moving = fixedmov(img, tfm) - mxshift = (100,100) #make sure this isn't too small - tform, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, rtol=0, fvalue=0.0003) + mxshift = (100, 100) #make sure this isn't too small + tform, mm = qd_translate(fixed, moving, mxshift; maxevals = 1000, rtol = 0, fvalue = 0.0003) tfmtest(tfm, tform) #Rigid transform SD = Matrix{Float64}(LinearAlgebra.I, 2, 2) - tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now + tfm = Translation(@SVector([14, 17])) ∘ LinearMap(RotMatrix(0.3)) #no distortion for now fixed, moving = fixedmov(centered(img), tfm) - mxshift = (100,100) #make sure this isn't too small + mxshift = (100, 100) #make sure this isn't too small mxrot = (0.5,) minwidth_rot = fill(0.002, 3) - tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002) + tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002) tfmtest(tfm, tform) #with anisotropic sampling SD = Matrix(Diagonal([0.5; 1.0])) - tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.3)*SD) + tfm = Translation(@SVector([14.3, 17.8])) ∘ LinearMap(SD \ RotMatrix(0.3) * SD) fixed, moving = fixedmov(centered(img), tfm) - tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002) + tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002) tfmtest(tfm, arrayscale(tform, SD)) #Affine transform - tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.01)) + tfm = Translation(@SVector([14, 17])) ∘ LinearMap(RotMatrix(0.01)) #make it harder with nonuniform scaling scale = @SMatrix [1.005 0; 0 0.995] SD = Matrix{Float64}(LinearAlgebra.I, 2, 2) - tfm = AffineMap(tfm.linear*scale, tfm.translation) - mxshift = (100,100) #make sure this isn't too small + tfm = AffineMap(tfm.linear * scale, tfm.translation) + mxshift = (100, 100) #make sure this isn't too small fixed, moving = fixedmov(centered(img), tfm) - tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=1000, rtol=0, fvalue=0.0002) + tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals = 1000, rtol = 0, fvalue = 0.0002) tfmtest(tfm, tform) #with anisotropic sampling SD = Matrix(Diagonal([0.5; 1.0])) - tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(RotMatrix(0.1)) #Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.01)*SD) + tfm = Translation(@SVector([14.3, 17.8])) ∘ LinearMap(RotMatrix(0.1)) #Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.01)*SD) scale = @SMatrix [1.005 0; 0 0.995] - tfm = AffineMap(tfm.linear*scale, tfm.translation) + tfm = AffineMap(tfm.linear * scale, tfm.translation) tfm = arrayscale(tfm, SD) fixed, moving = fixedmov(centered(img), tfm) - tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=10000, rtol=0, fvalue=0.0002, ndmax = 0.25) + tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals = 10000, rtol = 0, fvalue = 0.0002, ndmax = 0.25) tform2 = arrayscale(tform, SD) tfmtest(tfm, tform2) end #tests with standard images @@ -96,33 +96,33 @@ end #tests with standard images @testset "Quadratic interpolation (issue #7)" begin samplefrom(n) = rand(Poisson(n)) - img = restrict(restrict(testimage("cameraman")))[2:end-1,2:end-1] + img = restrict(restrict(testimage("cameraman")))[2:(end - 1), 2:(end - 1)] # Convert to "photons" so we can mimic shot noise np = 1000 # maximum number of photons per pixel - img = round.(Int, np.*gray.(img)) - fixed = samplefrom.(img) + img = round.(Int, np .* gray.(img)) + fixed = samplefrom.(img) moving = samplefrom.(img) ff = qsmooth(fixed) - tform, mm = qd_translate(fixed, moving, (5, 5); print_interval=typemax(Int)) - tformq, mmq = qd_translate(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int)) + tform, mm = qd_translate(fixed, moving, (5, 5); print_interval = typemax(Int)) + tformq, mmq = qd_translate(ff, moving, (5, 5); presmoothed = true, print_interval = typemax(Int)) @test all(abs.(tformq.translation) .< abs.(tform.translation)) - tform, mm = qd_rigid(fixed, moving, (5, 5), 0.1; print_interval=typemax(Int)) - tformq, mmq = qd_rigid(ff, moving, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int)) - @test norm(tformq.linear-I) < norm(tform.linear-I) + tform, mm = qd_rigid(fixed, moving, (5, 5), 0.1; print_interval = typemax(Int)) + tformq, mmq = qd_rigid(ff, moving, (5, 5), 0.1; presmoothed = true, print_interval = typemax(Int)) + @test norm(tformq.linear - I) < norm(tform.linear - I) @test norm(tformq.translation) < norm(tform.translation) - tform, mm = qd_affine(fixed, moving, (5, 5); print_interval=typemax(Int)) - tformq, mmq = qd_affine(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int)) + tform, mm = qd_affine(fixed, moving, (5, 5); print_interval = typemax(Int)) + tformq, mmq = qd_affine(ff, moving, (5, 5); presmoothed = true, print_interval = typemax(Int)) @test mmq < mm # Test that we exactly reconstruct `qsmooth` with `presmoothed=true` - tformq, mmq = qd_translate(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int)) + tformq, mmq = qd_translate(ff, fixed, (5, 5); presmoothed = true, print_interval = typemax(Int)) @test all(iszero, tformq.translation) - @test mmq < 1e-10 - tformq, mmq = qd_rigid(ff, fixed, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int)) - @test mmq < 1e-8 - tformq, mmq = qd_affine(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int)) - @test mmq < 1e-6 # on 32-bit systems this can't be 1e-8, not quite sure why + @test mmq < 1.0e-10 + tformq, mmq = qd_rigid(ff, fixed, (5, 5), 0.1; presmoothed = true, print_interval = typemax(Int)) + @test mmq < 1.0e-8 + tformq, mmq = qd_affine(ff, fixed, (5, 5); presmoothed = true, print_interval = typemax(Int)) + @test mmq < 1.0e-6 # on 32-bit systems this can't be 1e-8, not quite sure why end diff --git a/test/util.jl b/test/util.jl index 0fa11f6..503d742 100644 --- a/test/util.jl +++ b/test/util.jl @@ -11,56 +11,56 @@ using Unitful: μm, mm, cm, km, s img = rand(3, 10) ci = CartesianIndices(img) θ = RegisterQD.default_minrot(ci) - @test θ ≈ 0.01 rtol=0.1 + @test θ ≈ 0.01 rtol = 0.1 θ = RegisterQD.default_minrot(ci, [1 0; 0 2]) - @test θ ≈ 0.005 rtol=0.1 + @test θ ≈ 0.005 rtol = 0.1 θ = RegisterQD.default_minrot(ci, [3 0; 0 1]) - @test θ ≈ 0.007 rtol=0.1 + @test θ ≈ 0.007 rtol = 0.1 θ = RegisterQD.default_minrot(ci, [10 0; 0 1]) - @test θ ≈ 0.01/3 rtol=0.1 + @test θ ≈ 0.01 / 3 rtol = 0.1 img = rand(3, 10, 5) ci = CartesianIndices(img) θ = RegisterQD.default_minrot(ci) - @test θ ≈ 0.1/sqrt(3^2 + 10^2 + 5^2) rtol=1e-3 + @test θ ≈ 0.1 / sqrt(3^2 + 10^2 + 5^2) rtol = 1.0e-3 # array overload matches the CartesianIndices form SD = [1 0 0; 0 2 0; 0 0 3] @test RegisterQD.default_minrot(img) == RegisterQD.default_minrot(CartesianIndices(img)) @test RegisterQD.default_minrot(img, SD) == RegisterQD.default_minrot(CartesianIndices(img), SD) - @test RegisterQD.default_minrot(img, SD; Δc=0.2) == RegisterQD.default_minrot(CartesianIndices(img), SD; Δc=0.2) + @test RegisterQD.default_minrot(img, SD; Δc = 0.2) == RegisterQD.default_minrot(CartesianIndices(img), SD; Δc = 0.2) end @testset "getSD" begin #test that getSD deals with arbitrary dimensions - A2 = rand(10,10) + A2 = rand(10, 10) @test getSD(A2) == I - A3 = rand(10,10,10) + A3 = rand(10, 10, 10) @test getSD(A3) == I - A5 = rand(10,10,10,10,10) + A5 = rand(10, 10, 10, 10, 10) @test getSD(A5) == I Ax3 = AxisArray(A3, 1:1:10, 1:2:20, 1:3:30) - @test getSD(Ax3) == Diagonal([1.0,2.0,3.0]) + @test getSD(Ax3) == Diagonal([1.0, 2.0, 3.0]) #test that getSD works with test images mri = testimage("mri-stack.tif") @test getSD(mri) == Diagonal([1.0, 1.0, 5.0]) #test that getSD deals with images with arbitrary space directions - skewed = ImageMeta(rand(10,10,10)) - skewmatrix = rand(3,3) - skewed.spacedirections = (Tuple(skewmatrix[1,:]), Tuple(skewmatrix[2,:]), Tuple(skewmatrix[3,:])) + skewed = ImageMeta(rand(10, 10, 10)) + skewmatrix = rand(3, 3) + skewed.spacedirections = (Tuple(skewmatrix[1, :]), Tuple(skewmatrix[2, :]), Tuple(skewmatrix[3, :])) getSD(skewed) == skewmatrix #test that getSD can reconcile units of different magnitudes - badsampling = AxisArray(rand(10,10,10), (:x,:y,:z), (1mm, 2km, 3.4cm)) + badsampling = AxisArray(rand(10, 10, 10), (:x, :y, :z), (1mm, 2km, 3.4cm)) badsampling = ImageMeta(badsampling) - @test getSD(badsampling) == Diagonal([1, 2e6, 34]) + @test getSD(badsampling) == Diagonal([1, 2.0e6, 34]) #test that getSD ignores the time-axis - timedarray = AxisArray(rand(10,10,10), (:x, :y, :time), (1μm, 1μm, 1s)) - @test size(getSD(timedarray)) == (2,2) + timedarray = AxisArray(rand(10, 10, 10), (:x, :y, :time), (1μm, 1μm, 1s)) + @test size(getSD(timedarray)) == (2, 2) end @testset "qsmooth eltype" begin