Computes a Matrix3D homography from source and destination coplanar points

HomographyUtil.as

package math{

	import flash.geom.Matrix3D;
	import flash.geom.Point;
	/**
	 * 
	 * @author Mark Lundin
	 * 
	 * Based upon code provided by nicoptere - http://www.nicoptere.net/AS3/homographie/blog/Homography.as
	 * The findHomography now computes a Matrix3D that maps the transformation between two sets of complanar points.
	 * The transformation finds the a mapping from the source quad to a unit square, 
	 * and a second mapping from a unit square to the destination quad.
	 * An adjoint is then found for the second matrix and the two multiplied together. 
	 * 
	 * Standard vectors3d points can be transformed against this matrix. 
	 * Divide by the z component to map to a coordinate in the uv plane.
	 * 
	 */
	public class HomographyUtil
	extends Matrix3D
	{

		public static function findHomography( source:Vector.<Point>, destination:Vector.<Point> ):Matrix3D
		{
			/*
			 * This can probably be optimized. Its using a very simple and straightforward equation
			 * 
			 * Mapping is quad -> unit square -> quad. This should probably skip out the unit square and go straight from quad to quad
			 */
			
			var sourceHomography : Matrix3D = HomographyUtil.getSystem( destination );
			var destHomography : Matrix3D 	= HomographyUtil.adjoint( HomographyUtil.getSystem( source ) );
			destHomography.append( sourceHomography );
			
         	return destHomography;
			
		}

		public static function getSystem( points:Vector.<Point> ):Matrix3D
		{
			
			var sx:Number = (points[0].x - points[1].x) + (points[2].x - points[3].x);
			var sy:Number = (points[0].y - points[1].y) + (points[2].y - points[3].y);
			
			var dx1:Number = points[1].x - points[2].x;
			var dx2:Number = points[3].x - points[2].x;
			var dy1:Number = points[1].y - points[2].y;
			var dy2:Number = points[3].y - points[2].y;
		 
			var z:Number = (dx1 * dy2) - (dy1 * dx2);
			var g:Number = ((sx * dy2) - (sy * dx2)) / z;
			var h:Number = ((sy * dx1) - (sx * dy1)) / z;
		 
			var a:Number = points[1].x - points[0].x + g * points[1].x;
			var b:Number = points[3].x - points[0].x + h * points[3].x;
			var c:Number = points[0].x;
			var d:Number = points[1].y - points[0].y + g * points[1].y;
			var e:Number = points[3].y - points[0].y + h * points[3].y;	
			var f:Number = points[0].y;
			
			var homography:Vector.<Number> = new Vector.<Number>( 16, true );
			homography[0] = a;	homography[4] = b;	homography[8] = c;	homography[12] = 0;
			homography[1] = d;	homography[5] = e;	homography[9] = f;	homography[13] = 0;
			homography[2] = g;	homography[6] = h;	homography[10] = 1;	homography[14] = 0;
			homography[3] = 0;	homography[7] = 0;	homography[11] = 0;	homography[15] = 0;
         	
         	return new Matrix3D( homography );
         				
		}
		
		public static function adjoint ( matrix : Matrix3D ) : Matrix3D
		{
			
			var adj:Vector.<Number> = new Vector.<Number>( 16, true);
			var m:Vector.<Number> = matrix.rawData;
			
			var a:Number = m[0];
			var b:Number = m[1];
			var c:Number = m[2];
			
			var d:Number = m[4];
			var e:Number = m[5];
			var f:Number = m[6];
			
			var g:Number = m[8];
			var h:Number = m[9];
			
         	adj [ 0 ] 	= e-f*h; 	adj [ 4 ] 	= f*g-d; 	adj [ 8 ] 	= d*h-e*g;	adj [ 12 ]	= 0;
         	adj [ 1 ] 	= c*h-b; 	adj [ 5 ] 	= a-c*g; 	adj [ 9 ] 	= b*g-a*h;	adj [ 13 ]	= 0;
         	adj [ 2 ] 	= b*f-c*e; 	adj [ 6 ] 	= c*d-a*f;	adj [ 10 ] 	= a*e-b*d;	adj [ 14 ]	= 0;
         	adj [ 3 ]	= 0;		adj [ 7 ]	= 0;		adj	[ 11 ]	= 0;		adj [ 15 ]	= 1;
         	
         	return new Matrix3D( adj ); 
         	
		}
		
	}
}