Fix abs2 (#36)

Author: IanButterworth

src/RegisterQD.jl

@@ -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...)

src/affine.jl

@@ -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...)
@@ -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...)
@@ -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...)

src/rigid.jl

@@ -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)
@@ -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))
@@ -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...)

src/translations.jl

@@ -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))
@@ -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)

src/util.jl

@@ -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)

test/qd_random.jl

@@ -5,6 +5,7 @@ using RegisterQD.Images
 using RegisterQD.CoordinateTransformations
 using RegisterQD.Rotations
 using RegisterQD.RegisterMismatch
+using RegisterQD: _abs2
 using Random
 
 #import BlockRegistration, RegisterOptimize
@@ -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)
@@ -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)
@@ -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]
@@ -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)
@@ -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))))
 
@@ -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));
 
@@ -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