How to stitch images with very little overlap?

Question

I am trying to create a panorama using images with very little overlap, but I know the angle of the camera so I know exactly how much overlap there is and I know the order of the images so I know where each one belongs in the panorama. As a first pass I simply cocantenated the images together but the result is not good enough. Is there a way to crop the Bitmaps to trapezoids to eliminate (most of the) the overlap, then stretch the Bitmaps back to Rectangles before the concantenation? I know this will produce distortion during the stretching, and that a trapezoid is just a close approximation to how the Bitmap actually needs to be cropped, but I am hoping this will be good enough.

Answer 1

See Automatic Image Stitching with Accord.NET, with demo and source code.

Answer 2

What you want is referred to as a matrix transformation.

Here are some simple examples in C#/GDI+.

MSDN has some more in-depth descriptions.

I believe in the end you will be looking for a "Perspective Transformation", here is an SO question about that that might lead you down the right path.

I'm not worried about the bounty, this is a complex (and fun) topic and I don't have the time to work out a solution, I just hope this information is helpful. :)

Answer 3

The technique you're looking for is called Image Registration using an Affine Transform. This can be achieved in software by computing the matrix transform for image B that maps it on to image A. I assume you are trying to do this in software using Windows Forms & GDI+? Then the matrices you have available are 3x3 matrices, capable of Scale, translate, Rotate and Skew. This is often enough to create simple image registration and I used this technique successfully in a commercial software package (however was WPF).

To achieve Image registration using an Affine Transform, firstly you need a collection of control points in a pair of images to be registered. From this we can compute the 2D transformation to register the images? I have done this in WPF, which has a 3x3 Matrix can be defined using the System.Windows.Media.Matrix class, which has the following constructor:

Matrix(double m11, double m12, double m21, double m22, 
         double offsetX, double offsetY)

Note: GDI+ has a Matrix class whose constructor may differ but the principle is the same

The constructor arguments form the matrix as follows:

M11     M12     0
M21     M22     0
OffsetX OffsetY 1

Now if the input points are called X,Y and output U,V, the affine matrix transform, T, that maps X,Y onto U,V can be computed as follows:

U        =  X         *  T

[U1 V1 1] = [X1 Y1 1]   [A  B  0]
[U2 V2 1] = [X2 Y2 1] * [C  D  0]
[U3 V3 1] = [X3 Y3 1]   [Tx Ty 1]

This can also be simplified as follows:

U      =  X         *  T

[U1 V1] = [X1 Y1 1]   [A  B ]
[U2 V2] = [X2 Y2 1] * [C  D ]
[U3 V3] = [X3 Y3 1]   [Tx Ty]

or

X^-1         * U       = T
[X1 Y1 1]^-1   [U1 V1]   [A  B ]
[X2 Y2 1]    * [U2 V2] = [C  D ]
[X3 Y3 1]      [U3 V3]   [Tx Ty]

In english what this means is, given a list of points X,Y in image 1 that correspond to image 2, the inverse of the matrix X containing XY points multiplied by the matrix of corresponding points in image 2 gives you your matrix transform from Image 1 to 2

The output transform T contains A,B,C,D and Tx,Ty which correspond to M11,M12,M21,M22,OffsetX,OffsetY in the 3x3 affine matrix class (constructor above). However, if the X matrix and U matrix have more than 3 points, the solution is overdetermined and a least squares fit must be found. This is acheived using the Moores-Penrose Psuedo-inverse to find X^-1.

What does this mean in code? Well I coded my own Matrix3x3, Matrix3x2 classes and Control Point (x,y point) to handle the transformation then applied this to a WPF MatrixTransform on an element. In GDI+ you can do the same by applying the Matrix to the graphics pipeline before calling Graphics.DrawImage. Let's see how we can compute the transformation matrix.

The first class we need is the Matrix3x3 class:

public class Matrix3x3 : ICloneable
{
    #region Local Variables

    private double [] coeffs;

    private const int _M11 = 0;
    private const int _M12 = 1;
    private const int _M13 = 2;
    private const int _M21 = 3;
    private const int _M22 = 4;
    private const int _M23 = 5;
    private const int _M31 = 6;
    private const int _M32 = 7;
    private const int _M33 = 8;

    #endregion

    #region Construction

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x3"/> class.
    /// </summary>
    public Matrix3x3()
    {
        coeffs = new double[9];
    }

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x3"/> class.
    /// </summary>
    /// <param name="coefficients">The coefficients to initialise. The number of elements of the array should
    /// be equal to 9, else an exception will be thrown</param>
    public Matrix3x3(double[] coefficients)
    {
        if (coefficients.GetLength(0) != 9)
            throw new Exception("Matrix3x3.Matrix3x3()", "The number of coefficients passed in to the constructor must be 9");

        coeffs = coefficients;
    }

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x3"/> class. 
    /// </summary>
    /// <param name="m11">The M11 coefficient</param>
    /// <param name="m12">The M12 coefficien</param>
    /// <param name="m13">The M13 coefficien</param>
    /// <param name="m21">The M21 coefficien</param>
    /// <param name="m22">The M22 coefficien</param>
    /// <param name="m23">The M23 coefficien</param>
    /// <param name="m31">The M31 coefficien</param>
    /// <param name="m32">The M32 coefficien</param>
    /// <param name="m33">The M33 coefficien</param>
    public Matrix3x3(double m11, double m12, double m13, double m21, double m22, double m23, double m31, double m32, double m33)
    {
        // The 3x3 matrix is constructed as follows
        //
        // | M11 M12 M13 | 
        // | M21 M22 M23 | 
        // | M31 M32 M33 | 

        coeffs = new double[] { m11, m12, m13, m21, m22, m23, m31, m32, m33 };
    }

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x3"/> class. The IAffineTransformCoefficients
    /// passed in is used to populate coefficients M11, M12, M21, M22, M31, M32. The remaining column (M13, M23, M33)
    /// is populated with homogenous values 0 0 1.
    /// </summary>
    /// <param name="affineMatrix">The IAffineTransformCoefficients used to populate M11, M12, M21, M22, M31, M32</param>
    public Matrix3x3(IAffineTransformCoefficients affineTransform)
    {
        coeffs = new double[] { affineTransform.M11, affineTransform.M12, 0, 
                                affineTransform.M21, affineTransform.M22, 0, 
                                affineTransform.OffsetX, affineTransform.OffsetY, 1};
    }

    #endregion

    #region Public Properties

    /// <summary>
    /// Gets or sets the M11 coefficient
    /// </summary>
    /// <value>The M11</value>
    public double M11
    {
        get
        {
            return coeffs[_M11];
        }
        set
        {
            coeffs[_M11] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M12 coefficient
    /// </summary>
    /// <value>The M12</value>
    public double M12
    {
        get
        {
            return coeffs[_M12];
        }
        set
        {
            coeffs[_M12] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M13 coefficient
    /// </summary>
    /// <value>The M13</value>
    public double M13
    {
        get
        {
            return coeffs[_M13];
        }
        set
        {
            coeffs[_M13] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M21 coefficient
    /// </summary>
    /// <value>The M21</value>
    public double M21
    {
        get
        {
            return coeffs[_M21];
        }
        set
        {
            coeffs[_M21] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M22 coefficient
    /// </summary>
    /// <value>The M22</value>
    public double M22
    {
        get
        {
            return coeffs[_M22];
        }
        set
        {
            coeffs[_M22] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M23 coefficient
    /// </summary>
    /// <value>The M23</value>
    public double M23
    {
        get
        {
            return coeffs[_M23];
        }
        set
        {
            coeffs[_M23] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M31 coefficient
    /// </summary>
    /// <value>The M31</value>
    public double M31
    {
        get
        {
            return coeffs[_M31];
        }
        set
        {
            coeffs[_M31] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M32 coefficient
    /// </summary>
    /// <value>The M32</value>
    public double M32
    {
        get
        {
            return coeffs[_M32];
        }
        set
        {
            coeffs[_M32] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M33 coefficient
    /// </summary>
    /// <value>The M33</value>
    public double M33
    {
        get
        {
            return coeffs[_M33];
        }
        set
        {
            coeffs[_M33] = value;
        }
    }

    /// <summary>
    /// Gets the determinant of the matrix
    /// </summary>
    /// <value>The determinant</value>
    public double Determinant
    {
        get
        {
            //                                |a b c|
            // In general, for a 3X3 matrix   |d e f|
            //                                |g h i|
            //
            // The determinant can be found as follows:
            // a(ei-fh) - b(di-fg) + c(dh-eg)

            // Get coeffs
            double a = coeffs[_M11];
            double b = coeffs[_M12];
            double c = coeffs[_M13];
            double d = coeffs[_M21];
            double e = coeffs[_M22];
            double f = coeffs[_M23];
            double g = coeffs[_M31];
            double h = coeffs[_M32];
            double i = coeffs[_M33];
            double ei = e * i;
            double fh = f * h;
            double di = d * i;
            double fg = f * g;
            double dh = d * h;
            double eg = e * g;

            // Compute the determinant
            return (a * (ei - fh)) - (b * (di - fg)) + (c * (dh - eg));
        }
    }

    /// <summary>
    /// Gets a value indicating whether this matrix is singular. If it is singular, it cannot be inverted
    /// </summary>
    /// <value>
    ///     <c>true</c> if this instance is singular; otherwise, <c>false</c>.
    /// </value>
    public bool IsSingular
    {
        get
        {
            return Determinant == 0;
        }
    }

    /// <summary>
    /// Gets the inverse of this matrix. If the matrix is singular, this method will throw an exception
    /// </summary>
    /// <value>The inverse</value>
    public Matrix3x3 Inverse
    {
        get
        {
            // Taken from http://everything2.com/index.pl?node_id=1271704
            //                                                  a b c
            //In general, the inverse matrix of a 3X3 matrix    d e f
            //                                                  g h i

            //is 

            //        1                              (ei-fh)   (bi-ch)   (bf-ce)
            // -----------------------------   x     (fg-di)   (ai-cg)   (cd-af)
            // a(ei-fh) - b(di-fg) + c(dh-eg)        (dh-eg)   (bg-ah)   (ae-bd)

            // Get coeffs
            double a = coeffs[_M11];
            double b = coeffs[_M12];
            double c = coeffs[_M13];
            double d = coeffs[_M21];
            double e = coeffs[_M22];
            double f = coeffs[_M23];
            double g = coeffs[_M31];
            double h = coeffs[_M32];
            double i = coeffs[_M33];

            //// Compute often used components
            double ei = e * i;
            double fh = f * h;
            double di = d * i;
            double fg = f * g;
            double dh = d * h;
            double eg = e * g;
            double bi = b * i;
            double ch = c * h;
            double ai = a * i;
            double cg = c * g;
            double cd = c * d;
            double bg = b * g;
            double ah = a * h;
            double ae = a * e;
            double bd = b * d;
            double bf = b * f;
            double ce = c * e;
            double cf = c * d;
            double af = a * f;

            // Construct the matrix using these components
            Matrix3x3 tempMat = new Matrix3x3(ei - fh, ch - bi, bf - ce, fg - di, ai - cg, cd - af, dh - eg, bg - ah, ae - bd);

            // Compute the determinant
            double det = Determinant;

            if (det == 0.0)
            {
                throw new Exception("Matrix3x3.Inverse", "Unable to invert the matrix as it is singular");
            }

            // Scale the matrix by 1/determinant
            tempMat.Scale(1.0 / det);

            return tempMat;
        }
    }

    /// <summary>
    /// Gets a value indicating whether this matrix is affine. This will be true if the right column 
    /// (M13, M23, M33) is 0 0 1
    /// </summary>
    /// <value><c>true</c> if this instance is affine; otherwise, <c>false</c>.</value>
    public bool IsAffine
    {
        get
        {
            return (coeffs[_M13] == 0 && coeffs[_M23] == 0 && coeffs[_M33] == 1);
        }
    }

    #endregion

    #region Public Methods

    /// <summary>
    /// Multiplies the current matrix by the 3x3 matrix passed in
    /// </summary>
    /// <param name="rhs"></param>
    public void Multiply(Matrix3x3 rhs)
    {
        // Get coeffs
        double a = coeffs[_M11];
        double b = coeffs[_M12];
        double c = coeffs[_M13];
        double d = coeffs[_M21];
        double e = coeffs[_M22];
        double f = coeffs[_M23];
        double g = coeffs[_M31];
        double h = coeffs[_M32];
        double i = coeffs[_M33];

        double j = rhs.M11;
        double k = rhs.M12;
        double l = rhs.M13;
        double m = rhs.M21;
        double n = rhs.M22;
        double o = rhs.M23;
        double p = rhs.M31;
        double q = rhs.M32;
        double r = rhs.M33;

        // Perform multiplication. Formula taken from
        // http://www.maths.surrey.ac.uk/explore/emmaspages/option1.html

        coeffs[_M11] = a * j + b * m + c * p;
        coeffs[_M12] = a * k + b * n + c * q;
        coeffs[_M13] = a * l + b * o + c * r;
        coeffs[_M21] = d * j + e * m + f * p;
        coeffs[_M22] = d * k + e * n + f * q;
        coeffs[_M23] = d * l + e * o + f * r;
        coeffs[_M31] = g * j + h * m + i * p;
        coeffs[_M32] = g * k + h * n + i * q;
        coeffs[_M33] = g * l + h * o + i * r;
    }

    /// <summary>
    /// Scales the matrix by the specified scalar value
    /// </summary>
    /// <param name="scalar">The scalar.</param>
    public void Scale(double scalar)
    {
        coeffs[0] *= scalar;
        coeffs[1] *= scalar;
        coeffs[2] *= scalar;
        coeffs[3] *= scalar;
        coeffs[4] *= scalar;
        coeffs[5] *= scalar;
        coeffs[6] *= scalar;
        coeffs[7] *= scalar;
        coeffs[8] *= scalar;
    }

    /// <summary>
    /// Makes the matrix an affine matrix by setting the right column (M13, M23, M33) to 0 0 1
    /// </summary>
    public void MakeAffine()
    {
        coeffs[_M13] = 0;
        coeffs[_M23] = 0;
        coeffs[_M33] = 1;
    }

    #endregion

    #region ICloneable Members

    /// <summary>
    /// Creates a new object that is a copy of the current instance.
    /// </summary>
    /// <returns>
    /// A new object that is a copy of this instance.
    /// </returns>
    public object Clone()
    {
        double[] coeffCopy = (double[])coeffs.Clone();
        return new Matrix3x3(coeffCopy);
    }

    #endregion

    #region IAffineTransformCoefficients Members

    //
    // NB: M11, M12, M21, M22 members of IAffineTransformCoefficients are implemented within the
    // #region Public Properties directive
    //

    /// <summary>
    /// Gets or sets the Translation Offset in the X Direction
    /// </summary>
    /// <value>The M31</value>
    public double OffsetX
    {
        get
        {
            return coeffs[_M31];
        }
        set
        {
            coeffs[_M31] = value;
        }
    }

    /// <summary>
    /// Gets or sets the Translation Offset in the Y Direction
    /// </summary>
    /// <value>The M32</value>
    public double OffsetY
    {
        get
        {
            return coeffs[_M32];
        }
        set
        {
            coeffs[_M32] = value;
        }
    }

    #endregion
}

And a Matrix3x2 class

public class Matrix3x2 : ICloneable
{
    #region Local Variables

    private double[] coeffs;

    private const int _M11 = 0;
    private const int _M12 = 1;
    private const int _M21 = 2;
    private const int _M22 = 3;
    private const int _M31 = 4;
    private const int _M32 = 5;

    #endregion

    #region Construction

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x2"/> class.
    /// </summary>
    public Matrix3x2()
    {
        coeffs = new double[6];
    }

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x2"/> class.
    /// </summary>
    /// <param name="coefficients">The coefficients to initialise. The number of elements of the array should
    /// be equal to 6, else an exception will be thrown</param>
    public Matrix3x2(double[] coefficients)
    {
        if (coefficients.GetLength(0) != 6)
            throw new Exception("Matrix3x2.Matrix3x2()", 
                "The number of coefficients passed in to the constructor must be 6");

        coeffs = coefficients;
    }

    public Matrix3x2(double m11, double m12, double m21, double m22, double m31, double m32)
    {
        coeffs = new double[] { m11, m12, m21, m22, m31, m32 };
    }

    /// <summary>
    /// Initializes a new instance of the <see cref="Matrix3x2"/> class. The IAffineTransformCoefficients
    /// passed in is used to populate coefficients M11, M12, M21, M22, M31, M32.
    /// </summary>
    /// <param name="affineMatrix">The IAffineTransformCoefficients used to populate M11, M12, M21, M22, M31, M32</param>
    public Matrix3x2(IAffineTransformCoefficients affineTransform)
    {
        coeffs = new double[] { affineTransform.M11, affineTransform.M12, 
                                affineTransform.M21, affineTransform.M22, 
                                affineTransform.OffsetX, affineTransform.OffsetY};
    }

    #endregion

    #region Public Properties

    /// <summary>
    /// Gets or sets the M11 coefficient
    /// </summary>
    /// <value>The M11</value>
    public double M11
    {
        get
        {
            return coeffs[_M11];
        }
        set
        {
            coeffs[_M11] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M12 coefficient
    /// </summary>
    /// <value>The M12</value>
    public double M12
    {
        get
        {
            return coeffs[_M12];
        }
        set
        {
            coeffs[_M12] = value;
        }
    }


    /// <summary>
    /// Gets or sets the M21 coefficient
    /// </summary>
    /// <value>The M21</value>
    public double M21
    {
        get
        {
            return coeffs[_M21];
        }
        set
        {
            coeffs[_M21] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M22 coefficient
    /// </summary>
    /// <value>The M22</value>
    public double M22
    {
        get
        {
            return coeffs[_M22];
        }
        set
        {
            coeffs[_M22] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M31 coefficient
    /// </summary>
    /// <value>The M31</value>
    public double M31
    {
        get
        {
            return coeffs[_M31];
        }
        set
        {
            coeffs[_M31] = value;
        }
    }

    /// <summary>
    /// Gets or sets the M32 coefficient
    /// </summary>
    /// <value>The M32</value>
    public double M32
    {
        get
        {
            return coeffs[_M32];
        }
        set
        {
            coeffs[_M32] = value;
        }
    }

    #endregion

    #region Public Methods

    /// <summary>
    /// Transforms the the ILocation passed in and returns the result in a new ILocation
    /// </summary>
    /// <param name="location">The location to transform</param>
    /// <returns>The transformed location</returns>
    public ILocation Transform(ILocation location)
    {
        // Perform the following equation:
        //
        // | x y 1 |   | M11 M12 |   |(xM11 + yM21 + M31) (xM12 + yM22 + M32)|
        //           * | M21 M22 | = 
        //             | M31 M32 | 

        double x = location.X * coeffs[_M11] + location.Y * coeffs[_M21] + coeffs[_M31];
        double y = location.X * coeffs[_M12] + location.Y * coeffs[_M22] + coeffs[_M32];

        return new Location(x, y);
    }

    /// <summary>
    /// Multiplies the 3x3 matrix passed in with the current 3x2 matrix
    /// </summary>
    /// <param name="x">The 3x3 Matrix X</param>
    public void Multiply(Matrix3x3 lhs)
    {
        // Multiply the 3x3 matrix with the 3x2 matrix and store inside the current 2x3 matrix
        // 
        // [a b c]   [j k]   [(aj + bl + cn) (ak + bm + co)]
        // [d e f] * [l m] = [(dj + el + fn) (dk + em + fo)]
        // [g h i]   [n o]   [(gj + hl + in) (gk + hm + io)]

        // Get coeffs
        double a = lhs.M11;
        double b = lhs.M12;
        double c = lhs.M13;
        double d = lhs.M21;
        double e = lhs.M22;
        double f = lhs.M23;
        double g = lhs.M31;
        double h = lhs.M32;
        double i = lhs.M33;

        double j = coeffs[_M11];
        double k = coeffs[_M12];
        double l = coeffs[_M21];
        double m = coeffs[_M22];
        double n = coeffs[_M31];
        double o = coeffs[_M32];

        coeffs[_M11] = a * j + b * l + c * n;
        coeffs[_M12] = a * k + b * m + c * o;
        coeffs[_M21] = d * j + e * l + f * n;
        coeffs[_M22] = d * k + e * m + f * o;
        coeffs[_M31] = g * j + h * l + i * n;
        coeffs[_M32] = g * k + h * m + i * o;
    }

    #endregion

    #region ICloneable Members

    /// <summary>
    /// Creates a new object that is a copy of the current instance.
    /// </summary>
    /// <returns>
    /// A new object that is a copy of this instance.
    /// </returns>
    public object Clone()
    {
        double[] coeffCopy = (double[])coeffs.Clone();
        return new Matrix3x2(coeffCopy);
    }

    #endregion

    #region IAffineTransformCoefficients Members

    //
    // NB: M11, M12, M21, M22 members of IAffineTransformCoefficients are implemented within the
    // #region Public Properties directive
    //

    /// <summary>
    /// Gets or sets the Translation Offset in the X Direction
    /// </summary>
    /// <value>The M31</value>
    public double OffsetX
    {
        get
        {
            return coeffs[_M31];
        }
        set
        {
            coeffs[_M31] = value;
        }
    }

    /// <summary>
    /// Gets or sets the Translation Offset in the Y Direction
    /// </summary>
    /// <value>The M32</value>
    public double OffsetY
    {
        get
        {
            return coeffs[_M32];
        }
        set
        {
            coeffs[_M32] = value;
        }
    }

    #endregion
}

From these we can perform image registration with a list of points that correspond to the two images. To clarify what this means, lets say your panorama shots have certain features that are the same. Both have a cathedral spire, both have a tree. The points that register images A to B would be the X,Y locations in each image that correspond, ie: the XY location of the spire in both images would be one pair of points.

Now with this list of points we can compute our transform:

    public Matrix3x2 ComputeForwardTransform(IList<Point> baselineLocations, IList<Point> registerLocations)
{
    if (baselineLocations.Count < 3 || registerLocations.Count < 3)
    {
        throw new Exception("ComputeForwardTransform()",
            "Unable to compute the forward transform. A minimum of 3 control point pairs are required");
    }

    if (baselineLocations.Count != registerLocations.Count)
    {
        throw new Exception("ComputeForwardTransform()",
            "Unable to compute the forward transform. The number of control point pairs in baseline and registration results must be equal");
    }

    // To compute 
    //    Transform = ((X^T * X)^-1 * X^T)U = (X^T * X)^-1 (X^T * U)

    //    X^T * X =
    //    [ Sum(x_i^2)   Sum(x_i*y_i) Sum(x_i) ]
    //    [ Sum(x_i*y_i) Sum(y_i^2)   Sum(y_i) ]
    //    [ Sum(x_i)     Sum(y_i)     Sum(1)=n ]

    //    X^T * U =
    //    [ Sum(x_i*u_i) Sum(x_i*v_i) ]
    //    [ Sum(y_i*u_i) Sum(y_i*v_i) ]
    //    [ Sum(u_i)     Sum(v_i) ]

    IList<Point> xy = baselineLocations;
    IList<Point> uv = registerLocations;

    double a = 0, b = 0, c = 0, d = 0, e = 0, f = 0, g = 0, h = 0, n = xy.Count;
    double p = 0, q = 0, r = 0, s = 0, t = 0, u = 0;

    for (int i = 0; i < n; i++)
    {
        // Compute sum of squares for X^T * X
        a += xy[i].X * xy[i].X;
        b += xy[i].X * xy[i].Y;
        c += xy[i].X;
        d += xy[i].X * xy[i].Y;
        e += xy[i].Y * xy[i].Y;
        f += xy[i].Y;
        g += xy[i].X;
        h += xy[i].Y;

        // Compute sum of squares for X^T * U
        p += xy[i].X * uv[i].X;
        q += xy[i].X * uv[i].Y;
        r += xy[i].Y * uv[i].X;
        s += xy[i].Y * uv[i].Y;
        t += uv[i].X;
        u += uv[i].Y;
    }

    // Create matrices from the coefficients
    Matrix3x2 uMat = new Matrix3x2(p, q, r, s, t, u);
    Matrix3x3 xMat = new Matrix3x3(a, b, c, d, e, f, g, h, n);

    // Invert X
    Matrix3x3 xInv = xMat.Inverse;

    // Perform the multiplication to get the transform
    uMat.Multiply(xInv);

    // Matrix uMat now holds the image registration transform to go from the current result to baseline
    return uMat;
}

Finally, the above can be called as follows:

// where xy1, xy2, xy3 are control points in first image, and uv1, uv2, uv3 are // corresponding pairs in the second image Matrix3x2 result = ComputeForwardTransform(new [] {xy1, xy2, xy3}. new [] {uv1, uv2, uv3});

Anyway, I hope this is helpful to you. I realise its not GDI+ specific but does discuss how to register images using 3x3 transforms in detail, which can be used both in GDI+ and WPF. I actually have a code example deep down somewhere on my hard drive and would be happy to talk more if you need clarification on the above.

Below: Demo showing stiched images