Construct affine transforms from point-pairs (#87)

Given sets of point pairs `from_points => to_points`, compute the affine
transformation that best reproduces the observed mapping.

Closes #30

Author: timholy

README.md

@@ -171,6 +171,19 @@ defined by a composition of a translation and a linear transformation. An
 and a translation, e.g. `Translation(v) ∘ LinearMap(v)` (or any combination of
 `LinearMap`, `Translation` and `AffineMap`).
 
+`AffineMap`s can be constructed to fit point pairs `from_points => to_points`:
+
+```julia
+julia> from_points = [[0, 0], [1, 0], [0, 1]];
+
+julia> to_points   = [[1, 1], [3, 1], [1.5, 3]];
+
+julia> AffineMap(from_points => to_points)
+AffineMap([1.9999999999999996 0.4999999999999999; -5.551115123125783e-16 2.0], [0.9999999999999999, 1.0000000000000002])
+```
+
+The points can be supplied as a collection of vectors or as a matrix with points as columns.
+
 #### Perspective transformations
 
 The perspective transformation maps real-space coordinates to those on a virtual

src/affine.jl

@@ -209,3 +209,22 @@ recenter(trans::AbstractMatrix, origin::Union{AbstractVector, Tuple}) = recenter
 
 transform_deriv(trans::AffineMap, x) = trans.linear
 # TODO transform_deriv_params
+
+"""
+    AffineMap(from_points => to_points) → trans
+
+Create an Affine transformation that approximately maps the `from` points to the `to` points.
+At least `n+1` non-degenerate points are required to map an `n`-dimensional space.
+If there are more points than this, the transformation will be over-determined and a least-squares
+solution will be computed.
+"""
+function AffineMap((from_points,to_points)::Pair)
+    M = column_matrix(to_points) * pinv(column_matrix(from_points, 1))
+    AffineMap(M[:, 1:end-1], M[:, end])
+end
+
+column_matrix(points::AbstractMatrix) = points
+column_matrix(points) = reduce(hcat, points)
+
+column_matrix(points::AbstractMatrix, lastval) = vcat(points, fill(lastval, 1, axes(points, 2)))
+column_matrix(points, lastval) = reduce(hcat, [vcat(point, lastval) for point in points])

test/affine.jl

@@ -129,4 +129,20 @@ end
         @test expected.origin ≈ result.origin
         @test expected.direction ≈ result.direction
     end
+
+    @testset "construction from points" begin
+        M = [1.0 2.0; 3.0 4.0]
+        v = [-1.0, 1.0]
+        A = AffineMap(M,v)
+        from_points = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
+        to_points = map(A, from_points)
+        A2 = AffineMap(from_points => to_points)
+        @test A2 ≈ A
+        A2 = AffineMap(reduce(hcat, from_points) => reduce(hcat, to_points))
+        @test A2 ≈ A
+        from_points = ([0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0])
+        to_points = map(A, from_points)
+        A2 = AffineMap(from_points => to_points)
+        @test A2 ≈ A
+    end
 end