[RFC] Fix piracy (#51)

* Update to CoordinateTransformations.kabsch

`kabsch` has now been implemented in CoordinateTransformations. While
new features are ordinarily not breaking changes, for this package it is
due to name conflict. This PR resolves the name conflict, and bumps
`[compat]`.

Closes #49

* Eliminate piracy of CoordinateTransformations

This package extends several functions from CoordinateTransformations,
but it extends them in ways that do not involve types defined in this
package. Hence these methods are piracy.

`trans(A::AbstractMatrix)` could conceivably be defined in
CoordinateTransformations, so rather than merge this as-is it might be
better to consider submitting it there.

* Add Aqua tests

I've set these up to run only on CI, but if you prefer it's easier to
always run them. They just take a while.

* Eliminate piracy of CoordinateTransformations

This package extends several functions from CoordinateTransformations,
but it extends them in ways that do not involve types defined in this
package. Hence these methods are piracy.

`trans(A::AbstractMatrix)` could conceivably be defined in
CoordinateTransformations, so rather than merge this as-is it might be
better to consider submitting it there.

* Add Aqua tests

I've set these up to run only on CI, but if you prefer it's easier to
always run them. They just take a while.

* Pass Aqua tests

* Fix premature version bump

* Fix LinearAlgebra compat entry

* Switch to build_tform for local ICP methods

---------

Co-authored-by: Tom McGrath <[email protected]>

Author: timholy

Project.toml

@@ -1,7 +1,7 @@
 name = "GaussianMixtureAlignment"
 uuid = "f2431ed1-b9c2-4fdb-af1b-a74d6c93b3b3"
 authors = ["Tom McGrath <[email protected]> and contributors"]
-version = "0.2.3"
+version = "0.3.0"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -27,13 +27,16 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 GaussianMixtureAlignmentMakieExt = "Makie"
 
 [compat]
+Aqua = "0.8"
 Colors = "0.12, 0.13"
 CoordinateTransformations = "0.6.4"
 Distances = "0.10"
 ForwardDiff = "0.10"
 GenericLinearAlgebra = "0.3"
 GeometryBasics = "0.4, 0.5"
 Hungarian = "0.7"
+IntervalSets = "0.7"
+LinearAlgebra = "1.10"
 Makie = "0.21, 0.22"
 MakieCore = "0.6, 0.7, 0.8, 0.9"
 MutableConvexHulls = "0.2"
@@ -43,12 +46,14 @@ PairedLinkedLists = "0.2"
 Requires = "1.3"
 Rotations = "1.4, 1.5, 1.6, 1.7"
 StaticArrays = "1.5"
+Test = "1.10"
 julia = "1.10"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "ForwardDiff", "IntervalSets"]
+test = ["Aqua", "Test", "ForwardDiff", "IntervalSets"]

src/branchbound.jl

@@ -73,7 +73,7 @@ end
 #     maxevals      - the maximum number of objective function evaluations allowed before search termination
 #     maxstagnant   - the maximum number of `Block` splits allowed without improvement before search termination
 """
-    result = branchbound(x, y; nsplits=2, searchspace=nothing, 
+    result = branchbound(x, y; nsplits=2, searchspace=nothing,
                          rot=nothing, trl=nothing, blockfun=fullBlock, objfun=alignment_objective,
                          rtol=0.01, maxblocks=5e8, maxeva ls=Inf, maxstagnant=Inf, threads=false)
 
@@ -82,12 +82,12 @@ and `y`, using the [GOGMA algorithm](https://arxiv.org/abs/1603.00150).
 
 Returns a `GlobalAlignmentResult` that contains the maximized overlap of the two GMMs (the upperbound on the objective function),
 a lower bound on the alignment objective function, an `AffineMap` which aligns `x` with `y`, and information about the
-number of evaluations during the alignment procedure. 
-""" 
+number of evaluations during the alignment procedure.
+"""
 function branchbound(xinput::AbstractModel, yinput::AbstractModel;
                      nsplits=2, searchspace=nothing, blockfun=UncertaintyRegion, R=RotationVec(0.,0.,0.), T=SVector{3}(0.,0.,0.),
-                     nextblockfun=lowestlbblock, centerinputs=false, boundsfun=tight_distance_bounds, localfun=local_align, tformfun=AffineMap,
-                     atol=0.1, rtol=0, maxblocks=5e8, maxsplits=Inf, maxevals=Inf, maxstagnant=Inf, separatesplit=false)
+                     nextblockfun=lowestlbblock, centerinputs=false, boundsfun=tight_distance_bounds, localfun=local_align, tformfun::TF=AffineMap,
+                     atol=0.1, rtol=0, maxblocks=5e8, maxsplits=Inf, maxevals=Inf, maxstagnant=Inf, separatesplit=false) where TF
     x = xinput
     y = yinput
     if isodd(nsplits)
@@ -127,12 +127,12 @@ function branchbound(xinput::AbstractModel, yinput::AbstractModel;
     ub, bestloc = localfun(x, y, searchspace)
 
     progress = [(0, ub, bestloc)]
-    
+
     # split cubes until convergence
     ndivisions = 0
     sinceimprove = 0
     evalsperdiv = rot_trl_split ? length(x)*length(y)*2*nsplits^3 : length(x)*length(y)*nsplits^ndims
-    
+
     while !isempty(hull)
         if (length(hull) > maxblocks) || (ndivisions*evalsperdiv > maxevals) || (sinceimprove > maxstagnant) || (ndivisions > maxsplits)
             break
@@ -150,9 +150,9 @@ function branchbound(xinput::AbstractModel, yinput::AbstractModel;
 
         # if the best solution so far is close enough to the best possible solution, end
         if abs((ub - lb)/lb) < rtol || abs(ub-lb) < atol
-            tform = tformfun(bestloc)
+            tform = build_tform(tformfun, bestloc)
             if centerinputs
-                tform = centerx_tform ∘ tform ∘  inv(centery_tform) 
+                tform = centerx_tform ∘ tform ∘  inv(centery_tform)
             end
             return GlobalAlignmentResult(x, y, ub, lb, tform, bestloc, ndivisions*evalsperdiv, ndivisions, length(hull), sinceimprove, progress, "optimum within tolerance")
         end
@@ -222,17 +222,17 @@ function branchbound(xinput::AbstractModel, yinput::AbstractModel;
         end
     end
     if isempty(hull)
-        tform = tformfun(bestloc)
+        tform = build_tform(tformfun, bestloc)
         if centerinputs
             tform = centerx_tform ∘ tform ∘  inv(centery_tform)
         end
-        return GlobalAlignmentResult(x, y, ub, lb, tformfun(bestloc), bestloc, ndivisions*evalsperdiv, ndivisions, length(hull), sinceimprove, progress, "priority queue empty")
+        return GlobalAlignmentResult(x, y, ub, lb, build_tform(tformfun, bestloc), bestloc, ndivisions*evalsperdiv, ndivisions, length(hull), sinceimprove, progress, "priority queue empty")
     else
-        tform = tformfun(bestloc)
+        tform = build_tform(tformfun, bestloc)
         if centerinputs
             tform = centerx_tform ∘ tform ∘  inv(centery_tform)
         end
-        return GlobalAlignmentResult(x, y, ub, lowestlbnode(hull).data[1], tformfun(bestloc), bestloc, ndivisions*evalsperdiv, ndivisions, length(hull), sinceimprove, progress, "terminated early")
+        return GlobalAlignmentResult(x, y, ub, lowestlbnode(hull).data[1], build_tform(tformfun, bestloc), bestloc, ndivisions*evalsperdiv, ndivisions, length(hull), sinceimprove, progress, "terminated early")
     end
 end
 
@@ -296,22 +296,22 @@ function planefit(mgmm::AbstractIsotropicMultiGMM, R)
     return planefit(R * ptsmat)
 end
 
-function tiv_branchbound(   x::AbstractModel, 
-                            y::AbstractModel, 
-                            tivx::AbstractModel, 
-                            tivy::AbstractModel; 
-                            boundsfun=tight_distance_bounds, 
-                            rot_boundsfun=boundsfun, 
-                            trl_boundsfun=boundsfun, 
-                            localfun=local_align, 
+function tiv_branchbound(   x::AbstractModel,
+                            y::AbstractModel,
+                            tivx::AbstractModel,
+                            tivy::AbstractModel;
+                            boundsfun=tight_distance_bounds,
+                            rot_boundsfun=boundsfun,
+                            trl_boundsfun=boundsfun,
+                            localfun=local_align,
                             rot_localfun=localfun,
                             trl_localfun=localfun,
                             kwargs...)
     t = promote_type(numbertype(x),numbertype(y))
     p = t(π)
     z = zero(t)
     zeroTranslation = SVector{3}(z,z,z)
-    
+
     rot_res = rot_branchbound(tivx, tivy; localfun=rot_localfun, boundsfun=rot_boundsfun, kwargs...)
     rotblock = RotationRegion(RotationVec(rot_res.tform_params...), zeroTranslation, p)
     rotscore, rotpos = rot_localfun(tivx, tivy, rotblock)
@@ -339,8 +339,8 @@ function tiv_branchbound(   x::AbstractModel,
         min = trl_res.upperbound
         bestpos = (rot_res.tform_params...,  trl_res.tform_params...)
     end
- 
-    return TIVAlignmentResult(x, y, min, trl_res.lowerbound, AffineMap(bestpos), bestpos, 
+
+    return TIVAlignmentResult(x, y, min, trl_res.lowerbound, build_tform(AffineMap, bestpos), bestpos,
                               rot_res.obj_calls+trl_res.obj_calls, rot_res.num_splits+trl_res.num_splits,
                               rot_res.num_blocks+trl_res.num_blocks,
                               rot_res, trl_res)

src/draw.jl

@@ -96,7 +96,16 @@ function plot!(gd::GMMDisplay{<:NTuple{<:Any,<:AbstractIsotropicGMM}})
     return gd
 end
 
-function plot!(gd::GMMDisplay{<:NTuple{<:Any,<:AbstractIsotropicMultiGMM{N,T,K}}}) where {N,T,K}
+@recipe(MultiGMMDisplay, g) do scene
+    Theme(
+        display = :wire,
+        palette = DEFAULT_COLORS, 
+        color = nothing,
+        label = "",
+    )
+end
+
+function plot!(gd::MultiGMMDisplay{<:NTuple{<:Any,<:AbstractIsotropicMultiGMM{N,T,K}}}) where {N,T,K}
     mgmms = [gd[i][] for i=1:length(gd)]
     disp = gd[:display][]
     color = gd[:color][]

src/gogma/bounds.jl

@@ -25,8 +25,17 @@ function pairwise_consts(gmmx::AbstractIsotropicGMM, gmmy::AbstractIsotropicGMM,
     return pσ, pϕ
 end
 
-function pairwise_consts(mgmmx::AbstractMultiGMM{N,T,K}, mgmmy::AbstractMultiGMM{N,S,K}, interactions::Union{Nothing,Dict{Tuple{K,K},V}}=nothing) where {N,T,S,K,V <: Number}
-    t = promote_type(numbertype(mgmmx),numbertype(mgmmy), isnothing(interactions) ? numbertype(mgmmx) : V)
+function pairwise_consts(mgmmx::AbstractMultiGMM{N,T,K}, mgmmy::AbstractMultiGMM{N,S,K}, interactions::Nothing=nothing) where {N,T,S,K}
+    t = promote_type(T, S)
+    self_interactions = Dict{Tuple{K,K},t}()
+    for key in keys(mgmmx.gmms) ∩ keys(mgmmy.gmms)
+        self_interactions[(key,key)] = one(t)
+    end 
+    pairwise_consts(mgmmx, mgmmy, self_interactions)
+end
+
+function pairwise_consts(mgmmx::AbstractMultiGMM{N,T,K}, mgmmy::AbstractMultiGMM{N,S,K}, interactions::Dict{Tuple{K,K},V}) where {N,T,S,K,V <: Number}
+    t = promote_type(T, S, isnothing(interactions) ? T : V)
     xkeys = keys(mgmmx.gmms)
     ykeys = keys(mgmmy.gmms)
     if isnothing(interactions)

src/goicp/icp.jl

@@ -15,14 +15,15 @@ function iterate_kabsch(P, Q, wp=ones(size(P,2)), wq=ones(size(Q,2)); iterations
 
         prevscore = score
         tform = kabsch_matches(P, Q, matches, wp, wq)
-        score = squared_deviation(tform(P),Q,matches)
+        Pt = transform_columns(tform, P)
+        score = squared_deviation(Pt,Q,matches)
         if prevscore < score
             matches = prevmatches
             break
         end
 
         prevmatches = matches
-        matches = correspondence(tform(P),Q)
+        matches = correspondence(Pt,Q)
         if matches == prevmatches
             break
         end
@@ -63,7 +64,7 @@ end
 function local_matching_alignment(x::AbstractPointSet, y::AbstractPointSet, block::TranslationRegion; matching_fun = iterative_hungarian, kwargs...)
     tformedx = block.R*x + block.T
     matches = matching_fun(tformedx, y; kwargs...)
-    tform = kabsch(x, y, matches)
+    tform = kabsch_matches(x, y, matches)
     score = squared_deviation(tform(x), y, matches)
     params = (tform.translation...,)
     return (score, params)

src/goicp/kabsch.jl

@@ -14,9 +14,15 @@ end
 kabsch_centered_matches(P::PointSet, Q::PointSet, matches::AbstractVector{<:Tuple{Int,Int}}) = kabsch_centered_matches(P.coords, Q.coords, matches, P.weights, Q.weights);
 
 # transform DxN matrices
-function (tform::Translation)(A::AbstractMatrix)
-    return hcat([tform(A[:,i]) for i=1:size(A,2)]...)
+function transform_columns(tform::Translation, A::AbstractMatrix)
+    return reduce(hcat, [tform(A[:,i]) for i=1:size(A,2)])
 end
+function transform_columns(tform::AffineMap, A::AbstractMatrix)
+    l = LinearMap(tform.linear)
+    t = Translation(tform.translation)
+    transform_columns(t, l(A))
+end
+transform_columns(tform::Union{Translation,AffineMap}, P::PointSet) = PointSet(transform_columns(tform, P.coords), P.weights)
 
 
 function kabsch_matches(P,Q,matches::AbstractVector{<:Tuple{Int,Int}},wp=ones(size(P,2)),wq=ones(size(Q,2)))
@@ -27,7 +33,7 @@ end
 
 kabsch_matches(P::PointSet, Q::PointSet, matches::AbstractVector{<:Tuple{Int,Int}}) = kabsch_matches(P.coords, Q.coords, matches, P.weights, Q.weights);
 
-function kabsch_matches(P::AbstractMultiPointSet{N,T,K}, Q::AbstractMultiPointSet{N,T,K}, matchesdict, wp = weights(P), wq = weights(Q)) where {N,T,K}
+function kabsch_matches(P::AbstractMultiPointSet{N,T,K}, Q::AbstractMultiPointSet{N,T,K}, matchesdict::Dict, wp = weights(P), wq = weights(Q)) where {N,T,K}
     matchedP, matchedQ = matched_points(P,Q,matchesdict)
     w = Vector{T}()
 

src/goicp/local.jl

@@ -11,7 +11,7 @@ function align_local_points(P, Q; maxevals=1000, tformfun=AffineMap)
 
     # local optimization within the block
     function f(X)
-        tform = tformfun((X...,))
+        tform = build_tform(tformfun, X)
         score = squared_deviation(tform(P), Q)
         return score
     end
@@ -31,7 +31,7 @@ function iterate_local_alignment(P, Q; correspondence = hungarian_assignment, it
         it += 1
         matchedP, matchedQ = matched_points(P,Q,matches)
         score, tformparams = align_local_points(matchedP, matchedQ; tformfun=tformfun, kwargs...)
-        tform = tformfun(tformparams)
+        tform = build_tform(tformfun, tformparams)
         prevmatches = matches
         matches = correspondence(tform(P), Q)
         if matches == prevmatches
@@ -45,7 +45,7 @@ function iterate_local_alignment(P, Q, block; tformfun=AffineMap, kwargs...)
     block_tform = AffineMap(block.R, block.T)
     tformedP = block_tform(P)
     score, opt_tformparams = iterate_local_alignment(tformedP, Q; tformfun=tformfun, kwargs...)
-    opt_tform = tformfun(opt_tformparams)
+    opt_tform = build_tform(tformfun, opt_tformparams)
     tform = opt_tform ∘ block_tform
     tform = AffineMap(RotationVec(tform.linear), tform.translation)
     tformparams = tformfun === LinearMap ? (tform.linear.sx, tform.linear.sy, tform.linear.sz) :

src/goicp/rmsd.jl

@@ -3,7 +3,7 @@ squared_deviation(P::AbstractMatrix, Q::AbstractMatrix, matches::AbstractVector{
 squared_deviation(P::AbstractPointSet, Q::AbstractPointSet, matches::AbstractVector{<:Tuple{Int,Int}}) = squared_deviation(P.coords, Q.coords, matches, P.weights, Q.weights)
 squared_deviation(P::AbstractPointSet, Q::AbstractPointSet) = squared_deviation(P.coords, Q.coords, hungarian_assignment(P.coords,Q.coords), P.weights, Q.weights)
 
-function squared_deviation(P::AbstractMultiPointSet{N,T,K}, Q::AbstractMultiPointSet{N,T,K}, matchesdict) where {N,T,K}
+function squared_deviation(P::AbstractMultiPointSet{N,T,K}, Q::AbstractMultiPointSet{N,T,K}, matchesdict::Dict) where {N,T,K}
     sqdev = zero(T)
     for (key, matches) in matchesdict
         sqdev += squared_deviation(P.pointsets[key], Q.pointsets[key], matches)

src/rocs/rocsalign.jl

@@ -5,14 +5,14 @@ struct ROCSAlignmentResult{D,S,T,F<:AbstractAffineMap,X<:AbstractGMM{D,S},Y<:Abs
     tform::F
 end
 
-""" 
+"""
     m = second_moment(gmm, center, dim1, dim2)
 
 Returns the second order moment of `gmm`
 """
-function mass_matrix(positions::AbstractMatrix{<:Real}, 
-                     weights=ones(eltype(positions),size(positions,2)), 
-                     widths=zeros(eltype(positions),size(positions,2)), 
+function mass_matrix(positions::AbstractMatrix{<:Real},
+                     weights=ones(eltype(positions),size(positions,2)),
+                     widths=zeros(eltype(positions),size(positions,2)),
                      center=centroid(positions,weights))
     t = eltype(positions)
     npts = size(positions,2)
@@ -46,7 +46,7 @@ is made diagonal, and the GMM center of mass is made the origin.
 """
 function inertial_transforms(positions::AbstractMatrix{<:Real},
                              weights=ones(eltype(positions),size(positions,2)),
-                             widths=zeros(eltype(positions),size(positions,2)); 
+                             widths=zeros(eltype(positions),size(positions,2));
                              invert = false)
     com = centroid(positions, weights / sum(weights))
     massmat = mass_matrix(positions, weights, widths, com)
@@ -60,7 +60,7 @@ function inertial_transforms(positions::AbstractMatrix{<:Real},
 
     # first, align the the eigenvectors to the coordinate system axes
     # make sure that a reflection is not performed
-    if det(evecs) < 0 
+    if det(evecs) < 0
         evecs[:,end] = -evecs[:,end]
     end
 
@@ -100,11 +100,11 @@ function rocs_align(gmmmoving::AbstractGMM, gmmfixed::AbstractGMM; kwargs...)
 
     # combine the inertial transform with the subsequent alignment transform for the best result
     minoverlap, mindex = findmin([r[1] for r in results])
-    tformmoving = AffineMap(results[mindex][2]) ∘ tformsmoving[mindex]
+    tformmoving = build_tform(AffineMap, results[mindex][2]) ∘ tformsmoving[mindex]
 
-    # Apply the inverse of `tformfixed` to the optimized transformation 
+    # Apply the inverse of `tformfixed` to the optimized transformation
     alignment_tform = inv(tformfixed) ∘ tformmoving
-    
+
     # return the result
     return ROCSAlignmentResult(gmmmoving, gmmfixed, minoverlap, alignment_tform)
 end

src/tforms.jl

@@ -2,9 +2,9 @@ import CoordinateTransformations.AffineMap
 import CoordinateTransformations.LinearMap
 import CoordinateTransformations.Translation
 
-AffineMap(params::NTuple{6}) = AffineMap(RotationVec(params[1:3]...), SVector{3}(params[4:6]...))
-LinearMap(params::NTuple{3}) = LinearMap(RotationVec(params...))
-Translation(params::NTuple{3}) = Translation(SVector{3}(params))
+build_tform(::Type{AffineMap}, params::NTuple{6}) = AffineMap(RotationVec(params[1:3]...), SVector{3}(params[4:6]...))
+build_tform(::Type{LinearMap}, params::NTuple{3}) = LinearMap(RotationVec(params...))
+build_tform(::Type{Translation}, params::NTuple{3}) = Translation(SVector{3}(params))
 
 function affinemap_to_params(tform::AffineMap)
     R = RotationVec(tform.linear)