Skip to content

Apply runic formatting #56

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
kdw503 merged 1 commit into from
May 26, 2026
Merged

Apply runic formatting #56

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
14 changes: 7 additions & 7 deletions src/RegisterQD.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
120 changes: 68 additions & 52 deletions src/affine.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -40,37 +40,37 @@ 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
# which defines `f(x)` in terms of this function.
# 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

"""
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
26 changes: 13 additions & 13 deletions src/gridsearch.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
Loading
Loading