Prokrustes: optimal affine transformation between two aligned 2D-point sets

Prokrustes.java

package de.unikn.algo.nexus.exsite.impl.rendering.pointpatternmatching;

import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

/**
 * This computes registration parameters for an optimal affine transformation
 * between 2 matched point sets towards the least squared error.
 * <p>
 * The implementation is based on <br />
 * <i> Chang, Shih-Hsu, et al.
 * "Fast algorithm for point pattern matching: invariant to translations, rotations and scale changes."
 * Pattern recognition 30.2 (1997): 311-320.</i>
 * </p>
 * 
 * @author Daniel Weidele (daniel.weidele@uni-konstanz.de)
 * @version 14.11.2014
 */
public class Prokrustes
{
	public static void main(String[] args)
	{
		List<Point2D.Double> A = Arrays.asList(new Point2D.Double(1, 1), new Point2D.Double(5, 7));
		List<Point2D.Double> B = Arrays.asList(new Point2D.Double(8, 11), new Point2D.Double(8, 9));

		System.out.println(e(t(A, r(A, B)), B) < 0.000001d ? "SUCCESS" : "FAIL");
	}

	/**
	 * Computes the {@link Registration} parameters of an affine transformation
	 * that leads to least squared error between A and B. The lists of points
	 * must be equal in length. Points are matched according to positions in
	 * their respective lists.
	 * 
	 * @param A
	 *            list of points to match to B
	 * @param B
	 *            list of points to matched by A
	 * @return {@link Registration} parameters of the transformation that leads
	 *         to least squared error between A and B
	 */
	public static Registration r(List<Point2D.Double> A, List<Point2D.Double> B)
	{
		int k = A.size();
		double m_xa = 0, m_ya = 0, m_xb = 0, m_yb = 0, lAPlusB = 0, lAMinusB = 0, lA = 0;
		for (int i = 0; i < k; i++)
		{
			m_xa += A.get(i).x;
			m_ya += A.get(i).y;
			m_xb += B.get(i).x;
			m_yb += B.get(i).y;
			lAPlusB += A.get(i).x * B.get(i).x + A.get(i).y * B.get(i).y;
			lAMinusB += A.get(i).x * B.get(i).y - A.get(i).y * B.get(i).x;
			lA += Math.pow(A.get(i).x, 2) + Math.pow(A.get(i).y, 2);
		}
		double det = k * lA - Math.pow(m_xa, 2) - Math.pow(m_ya, 2);
		Registration r = new Registration();
		r.tx = (lA * m_xb - m_xa * lAPlusB + m_ya * lAMinusB) / det;
		r.ty = (lA * m_yb - m_ya * lAPlusB - m_xa * lAMinusB) / det;
		r.sct = (-m_xa * m_xb - m_ya * m_yb + k * lAPlusB) / det;
		r.sst = (m_ya * m_xb - m_xa * m_yb + k * lAMinusB) / det;
		return r;
	}

	/**
	 * Applies the affine transformation based on the {@link Registration}
	 * parameters.
	 * 
	 * @param A
	 *            the list of points to be transformed
	 * @param r
	 *            the {@link Registration} parameters
	 * @return a copy of transformed points A
	 */
	public static List<Point2D.Double> t(List<Point2D.Double> A, Registration r)
	{
		List<Point2D.Double> T = new ArrayList<Point2D.Double>();
		for (Point2D.Double a : A)
		{
			Point2D.Double t = new Point2D.Double();
			t.x = r.tx + a.x * r.sct - a.y * r.sst;
			t.y = r.ty + a.y * r.sct + a.x * r.sst;
			T.add(t);
		}
		return T;
	}

	/**
	 * Computes the least squared error between point sets A and B. The lists
	 * must have equal length.
	 * 
	 * @param A
	 *            list of points
	 * @param B
	 *            list of points
	 * @return the least squared error between A and B
	 */
	public static double e(List<Point2D.Double> A, List<Point2D.Double> B)
	{
		int k = A.size();
		double e = 0;
		for (int i = 0; i < k; i++)
		{
			e += Math.sqrt(Math.pow(A.get(i).x - B.get(i).x, 2) + Math.pow(A.get(i).y - B.get(i).y, 2));
		}
		return e;
	}

	/**
	 * Registration parameters for an affine transformation.
	 * 
	 * @author Daniel Weidele
	 * @version 14.11.2014
	 */
	public static class Registration
	{
		private double tx;
		private double ty;
		private double sct;
		private double sst;

		/**
		 * Getter for x translation
		 * 
		 * @return x translation for the affine transformation
		 */
		public double getTx()
		{
			return this.tx;
		}

		/**
		 * Getter for y translation
		 * 
		 * @return y translation for the affine transformation
		 */
		public double getTy()
		{
			return this.ty;
		}

		/**
		 * Getter for scaling
		 *
		 * @return scaling
		 */
		public double getScaling()
		{
			return Math.sqrt(Math.pow(this.sct, 2) + Math.pow(this.sst, 2));
		}

		/**
		 * Getter for rotation (in degrees)
		 *
		 * @return rotation (in degrees)
		 */
		public double getRotation()
		{
			return Math.toDegrees(Math.acos(this.sct / getScaling()));
		}

		@Override
		public String toString()
		{
			return this.tx + "," + this.ty + "," + getScaling() + "," + getRotation();
		}
	}
}

references.bib

@article{chang1997fast,
  title={Fast algorithm for point pattern matching: invariant to translations, rotations and scale changes},
  author={Chang, Shih-Hsu and Cheng, Fang-Hsuan and Hsu, Wen-Hsing and Wu, Guo-Zua},
  journal={Pattern recognition},
  volume={30},
  number={2},
  pages={311--320},
  year={1997},
  publisher={Elsevier}
}

@misc{weidele2014prokrustes,
    author={Daniel Weidele},
    title={{Optimal Affine Transformation of Aligned 2D-Point Sets in Java}},
    year={2014},
    howpublished={\url{https://gist.github.com/dweidele/83e4203141da2f2c8733}},
    note={[{Online at GitHub Gist}; accessed dd.MM.yyyy]}
}