Skip to content

Fix abs2 #36

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 27, 2021
Merged

Fix abs2 #36

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
331ed7b
fix abs2
IanButterworth Jul 27, 2021
cd04bd8
add import
IanButterworth Jul 27, 2021
d2e5d40
manually remove type promotion
IanButterworth Jul 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/RegisterQD.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ function qd_rigid(fixed, moving, mxshift::VecLike, mxrot::Union{Number,VecLike},
end

function qd_affine(fixed, moving, mxshift, linmins, linmaxs, SD;
thresh=0.5*sum(abs2.(fixed[.!(isnan.(fixed))])),
thresh=0.5*sum(_abs2.(fixed[.!(isnan.(fixed))])),
initial_tfm=IdentityTransformation(),
print_interval=100,
kwargs...)
Expand Down
6 changes: 3 additions & 3 deletions src/affine.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ 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))])),
thresh=0.1*sum(_abs2.(fixed[.!(isnan.(fixed))])),
minwidth=default_lin_minwidths(moving),
maxevals=5e4,
kwargs...)
Expand Down Expand Up @@ -117,7 +117,7 @@ As a consequence, any large translations must be supplied with reasonable accura
function qd_affine_fine(fixed, moving, linmins, linmaxs;
SD=I,
initial_tfm=IdentityTransformation(),
thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])),
thresh=0.1*sum(_abs2.(fixed[.!(isnan.(fixed))])),
minwidth_mat=default_lin_minwidths(fixed)./10,
maxevals=5e4,
kwargs...)
Expand Down Expand Up @@ -184,7 +184,7 @@ overlap between the two images; with non-zero `thresh`, it is not permissible to
function qd_affine(fixed, moving, mxshift, linmins, linmaxs;
presmoothed=false,
SD=I,
thresh=0.5*sum(abs2.(fixed[.!(isnan.(fixed))])),
thresh=0.5*sum(_abs2.(fixed[.!(isnan.(fixed))])),
initial_tfm=IdentityTransformation(),
print_interval=100,
kwargs...)
Expand Down
6 changes: 3 additions & 3 deletions src/rigid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ end
function qd_rigid_coarse(fixed, moving, mxshift, mxrot, minwidth_rot;
SD=I,
initial_tfm=IdentityTransformation(),
thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])),
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)
Expand All @@ -104,7 +104,7 @@ end
function qd_rigid_fine(fixed, moving, mxrot, minwidth_rot;
SD=I,
initial_tfm=IdentityTransformation(),
thresh=0.1*sum(abs2.(fixed[.!(isnan.(fixed))])),
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))
Expand Down Expand Up @@ -179,7 +179,7 @@ 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))])),
thresh=0.1*sum(_abs2.(fixed[.!(isnan.(fixed))])),
initial_tfm=IdentityTransformation(),
print_interval=100,
kwargs...)
Expand Down
6 changes: 3 additions & 3 deletions src/translations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ end
function qd_translate_fine(fixed, moving;
initial_tfm=IdentityTransformation(),
minwidth=fill(0.01, ndims(fixed)),
thresh=0.1*sum(abs2.(fixed[.!(isnan.(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))
Expand Down Expand Up @@ -50,13 +50,13 @@ If you have a good initial guess at the solution, pass it with the `initial_tfm`
`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.

If the `crop` keyword arg is `true` then `fixed` is cropped by `mxshift` (after the optional `initial_tfm`) on all sides
If the `crop` keyword arg is `true` then `fixed` is cropped by `mxshift` (after the optional `initial_tfm`) on all sides
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))])),
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)
Expand Down
5 changes: 5 additions & 0 deletions src/util.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -240,3 +240,8 @@ 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...)
end

# abs2 was removed in https://github.com/JuliaGraphics/ColorVectorSpace.jl/pull/131
# this is the current recommended replacement
# TODO: follow https://github.com/JuliaGraphics/ColorVectorSpace.jl/issues/157 and adjust for changes
_abs2(c) = mapreducec(v->float(v)^2, +, float(zero(eltype(c))), c)
15 changes: 8 additions & 7 deletions test/qd_random.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ using RegisterQD.Images
using RegisterQD.CoordinateTransformations
using RegisterQD.Rotations
using RegisterQD.RegisterMismatch
using RegisterQD: _abs2
using Random

#import BlockRegistration, RegisterOptimize
Expand All @@ -23,7 +24,7 @@ using Test, TestImages
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))]))
thresh = 0.1 * sum(_abs2.(fixed[.!(isnan.(fixed))]))
mxshift = (10,10)

tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, rtol=0)
Expand All @@ -41,7 +42,7 @@ using Test, TestImages
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))]))
thresh = 0.1 * sum(_abs2.(fixed[.!(isnan.(fixed))]))
mxshift = (5,5,5)

tfm, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, thresh=thresh, rtol=0)
Expand All @@ -56,7 +57,7 @@ using Test, TestImages
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))]))
thresh = 0.1 * sum(_abs2.(fixed[.!(isnan.(fixed))]))
mxshift = (10,10)
mxrot = pi/90
minwidth_rot = [0.0002]
Expand All @@ -73,7 +74,7 @@ using Test, TestImages
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))]))
thresh = 0.1 * sum(_abs2.(fixed[.!(isnan.(fixed))]))
mxshift = (5,5,5)
mxrot = [pi/90; pi/90; pi/90]
minwidth_rot = fill(0.0002, 3)
Expand All @@ -95,7 +96,7 @@ using Test, TestImages
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))]))
thresh = 0.5 * sum(_abs2.(fixed[.!(isnan.(fixed))]))
mxshift = (5,5)
SD = SDiagonal(@SVector(ones(ndims(fixed))))

Expand All @@ -118,7 +119,7 @@ using Test, TestImages
#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))])));
#thresh = 0.1 * (sum(_abs2.(fixed[.!(isnan.(fixed))]))+sum(_abs2.(moving[.!(isnan.(moving))])));
#mxshift = (10,10,10)
#SD = eye(ndims(fixed));

Expand All @@ -141,7 +142,7 @@ using Test, TestImages
#inds = intersect.(axes(moving), axes(newfixed))
#fixed = newfixed[inds...]
#moving = moving[inds...]
#thresh = 0.5 * sum(abs2.(fixed[.!(isnan.(fixed))]));
#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
Expand Down