Determining affine transforms from three points

Status: Done
Confidence: Certain

Math in this page not rendering? See the fix

One of my labmates needed to visit a grid of points on a microscopy stage, something the software provided with the microscope wouldn’t handle. I wrote a program to generate grids of points for her. The heart of it is this calculation that three distinct points in the plane determine an affine transformation. This is true in two dimensions. In general, you need a nondegenerate simplex, i.e., in three dimensions you need four non-coplanar points.

In what follows, the $i$thth component of vector $v$ is written $v.i$, where $i = 0$ or 1 (we have only two indices since we are working in two dimensions). All the source code is Scheme.

We have a grid of positions with $N+1$ columns labelled 0 to $N$ and $M+1$ rows labelled 0 to $M$. It is worth labelling from 0, as it simplifies the mathematics dramatically. The only reason not to is that MATLAB begins counting at 1 in the old FORTRAN manner instead of at 0 as in all other languages. I will show the version starting at 1 and going to $N$ or $M$ at the end. We transform the labels for the grid positions into a grid of coordinates on a chip we are to observe under the microscope. The grid may be stretched, translated, and rotated from our labelling coordinate system. We can always make a translation tau to make the origins of the two systems coincide. Then we have two vector spaces (since we have found an origin) and the transformation between them is a linear function $T$. The chip coordinates $p$ are related to the labels $q$ by a transformation of the form

$p = T \cdot q + \tau$

or in coordinate form

$p.i = \langle\sum j : j = 0,1 : T.i.j \cdot q.j \rangle + \tau.i$

This is known as an affine transform. There are six unknown numbers in this equation, four components for T and two for tau. Therefore we need three points to calculate them. We choose the bottom left, bottom right, and top right labels of our grid. Their coordinates in the two systems are

These points give us a set of linear equations

$\begin{pmatrix} a.i \\ b.i \\ c.i \end{pmatrix} = T.i.0 \cdot \begin{pmatrix}0 \\ N \\ N \end{pmatrix} + T.i.1 \cdot \begin{pmatrix} 0 \\ 0 \\ M \end{pmatrix} + \tau.i \cdot \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}$

Now we turn to actually programming the transform. I am representing points in Scheme as functions of an integer argument just as I do mathematically

(define (point x y)
  (lambda (i)
    (if (= i 0) x y)))
(define (show-point p) 
  (list (p 0) (p 1)))

Similarly, $T$ will be a function of two arguments and tau a function of one. Given such functions, we can produce a function which does the point transform with

(define (generate-point-transform T tau)
  (lambda (p)
    (apply point (map (lambda (i) (+ (* (T i 0) (p 0)) 
                     (* (T i 1) (p 1))
                     (tau i))) '(0 1)))))

To build $T$’s and $\tau$’s, we use

(define (generate-T a b c M N)
  (lambda (i j)
    (if (= j 0) (/ (- (b i) (a i)) N)
    (/ (- (c i) (b i)) M))))

(define (generate-tau a b c M N)
  (lambda (i)
    (a i)))

We need to be able to produce grids of points to work with. generate-grid does exactly that.

(define (rem a b)
  (- a (* (floor (/ a b)) b)))

(define (seq a b)
  (if (< b a) '() (cons b (seq a (- b 1)))))

(define (generate-grid max-row max-col)
  (map (lambda (n) (point (floor (/ n (+ max-col 1))) (rem n (+ max-col 1)))) 
       (seq 0 (- (* (+ max-row 1) (+ max-col 1)) 1))))

In the following we create a transform with generate-point-transform and bind it to F, then use it to transform a number a list of points. Note that we pass $M=1$ and $N=1$ to generate-T and generate-tau. Remember that $M$ and $N$ are the last label of the grid, not the number of rows and columns!

(let ((F (generate-point-transform
      (generate-T (point 0 0) 
              (point (/ 1 (sqrt 2)) (/ 1 (sqrt 2))) 
              (point 0 (sqrt 2)) 
              1 1)
      (generate-tau (point 0 0)
                 (point (/ 1 (sqrt 2)) (/ 1 (sqrt 2)))
                 (point 0 (sqrt 2))
                 1 1))))
  (map show-point (map F (generate-grid 3 3))))

This approach is much easier than trying to figure out the geometry by hand, as it requires not particular insight or labor, just a clean calculation.

Appendix: Beginning at 1 instead of 0

If you are stuck in an environment where you must begin counting at 1, the matrices above become rather more complicated. To begin with, our points change to

and the equation that defines T and $\tau$ becomes

$\begin{pmatrix}a.i \\ b.i \\ c.i\end{pmatrix} = \begin{pmatrix}1 & 1 & 1 \\ N & 1 & 1 \\ N & M & 1\end{pmatrix} \cdot \begin{pmatrix}T.i.0 \\ T.i.1 \\ \tau.i\end{pmatrix}$
= { take the inverse of the central matrix in Maxima, and left multiply both sides with it }
$\begin{pmatrix} T.i.0 \\ T.i.1 \\ \tau.i \end{pmatrix} = \frac{1}{ M\cdot N - N - M + 1} \cdot \begin{pmatrix} 1-M & -(1-M) & 0 \\ 0 & 1-N & -(1-N) \\ (M-1)\cdot N & N-M & 1-N \end{pmatrix} \cdot \begin{pmatrix}a.i \\ b.i \\ c.i\end{pmatrix}$
= { carry out the multiplication }
$\begin{pmatrix} T.i.0 \\ T.i.1 \\ \tau.i \end{pmatrix} = \frac{1}{M\cdot N - N - M + 1} \cdot \begin{pmatrix} (M-1) \cdot (b.i - a.i) \\ (N-1) \cdot (c.i - b.i) \\ (M\cdot N - N) \cdot a.i + (N-M)\cdot b.i + (1-N)\cdot c.i \end{pmatrix}$

And we can implement it by providing functions to generate $T$ and $\tau$.

(define (denominator M N) (+ (* M N) (- N) (- M) 1))

(define (generate-base-1-T a b c M N)
  (lambda (i j)
    (/ (if (= j 0) (* (- M 1) (- (b i) (a i)))
       (* (- N 1) (- (c i) (b i))))
       (denominator M N))))

(define (generate-base-1-tau a b c M N)
  (lambda (i)
    (/ (+ (* (a i) (- (* M N) N)) (* (b i) (- N M)) (* (c i) (- 1 N)))
       (denominator M N))))

At this point the reader should be wondering what person in their right mind would count from 1 instead of 0. Yet it does work. For example,

(let ((F (generate-point-transform
      (generate-base-1-T (point 0 0) 
              (point (/ 1 (sqrt 2)) (/ 1 (sqrt 2))) 
              (point 0 (sqrt 2)) 
              2 2)
      (generate-base-1-tau (point 0 0)
                 (point (/ 1 (sqrt 2)) (/ 1 (sqrt 2)))
                 (point 0 (sqrt 2))
                 2 2))))
  (map show-point (map F (list (point 1 1) (point 2 1) (point 2 2) (point 1 2)))))

« Back to Miscellany | Home