Subject: Image Registration - Affine Transformation

From: [email-masked] (Michael Wirth)
Date: 4 Oct 1994 20:48:44 GMT

#include <stdio.h>
#include <stdlib.h>
#include <fstream.h>
#include <math.h>
#include "imagematrix.h"
#include "inverse.h"
#include "error.h"

#define PI 3.1415926535979

//  Image Registration by Affine Transformation
//  RST - Rotation, Scaling (Uniform) and Translation 

void Find_Properties(MPoints A, MPoints B, int n, float &a, float &b, float &h, 
				float &k)
{
    matrix4x4 X;
    float Y[4][1], W[4][1];
    float s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0,s6=0.0,s7=0.0,temp,temp1;
    float theta,alpha;
    int i,j;
 
    for (i=0; i<n; i++)
    {
	s1 = s1 + (pow(A[i][0],2)+pow(A[i][1],2));
	s2 = s2 + A[i][0];
	s3 = s3 + A[i][1];
    }
    X[0][1]=0.0; X[1][0]=0.0; X[3][2]=0.0; X[2][3]=0.0;
    X[0][0]=s1; X[1][1]=s1; 
    X[2][0]=s2; X[3][1]=s2; X[0][2]=s2; X[1][3]=s2;
    X[3][0]=s3; X[2][1]=-s3; X[1][2]=-s3; X[0][3]=s3;
    X[2][2]=n; X[3][3]=n;

    for (j=0; j<n; j++)
    {
	s4 = s4 + ((A[j][0]*B[j][0])+(A[j][1]*B[j][1]));
	s5 = s5 + ((B[j][1]*A[j][0])-(A[j][1]*B[j][0]));
	s6 = s6 + B[j][0];
	s7 = s7 + B[j][1];
    }
    Y[0][0]=s4; Y[1][0]=s5; Y[2][0]=s6; Y[3][0]=s7;

    InverseMatrix(X,4);
    for (i=0; i<4; i++)
    {
	temp = 0.0;
	for (j=0; j<4; j++)
	    temp = temp + (X[i][j]*Y[j][0]);
	W[i][0] = temp;
    }
    a = W[0][0];
    b = W[1][0];
    h = W[2][0];
    k = W[3][0];

    temp1 = atan(b/a);
    theta = (180.0*temp1)/PI;

    alpha = sqrt(pow(a,2)+pow(b,2));

    cout << "Scale   = " << alpha << endl;
    cout << "Angle   = " << theta << endl;
    cout << "h       = " << h << endl;
    cout << "v       = " << k << endl;
}

void Match_Images(AImage A, int nr1, int nc1, AImage C, int nr3, int nc3, 
		  AImage B, int &nr2, int &nc2, MPoints u, MPoints v, int cpr)
{
    AImage tempImage;
    MPoints w;
    matrix2x2 AB;
    float theta, h,k,a,b,tempx,tempy,tx,ty,DSum=0.0;
    float tempX, tempY,tempA,tempB,tempC,Q,rx,ry;
    int nrow,ncolumn, maxx=0,minx=400,maxy=0,miny=400;
    int i,j,X,Y,count=0;

    Find_Properties(u,v,cpr,a,b,h,k);

    for (i=0; i<nr1; i++)
	for (j=0; j<nc1; j++)
	{
	    tempx = (i*a)-(j*b)+h;
	    tempy = (i*b)+(j*a)+k;
            if (tempx < minx)
		minx = (int)tempx;
	    if (tempx > maxx)
		maxx = (int)tempx;
            if (tempy < miny)
		miny = (int)tempy;
	    if (tempy > maxy)
		maxy = (int)tempy;
	} 

    cout << "MinX = " << minx << " MaxX = " << maxx << endl;
    cout << "MinY = " << miny << " MaxY = " << maxy << endl;

    for (i=0; i<nr1; i++)
    {
	w[i][0] = (u[i][0]*a)-(u[i][1]*b)+h;
	w[i][1] = (u[i][0]*b)+(u[i][1]*a)+k;
    }

    Match_Error(v,w,cpr);

    nrow = maxx-minx;
    ncolumn = maxy-miny;
    nr2 = nrow;
    nc2 = ncolumn;
 
    AB[0][0] = a;
    AB[0][1] = -b;
    AB[1][0] = b;
    AB[1][1] = a;

    InverseMatrix2x2(AB,2);

    for (i=0; i<nrow; i++)
	for (j=0; j<ncolumn; j++)
	{
	    tx = (i-abs(minx)) - h;
	    ty = (j-abs(miny)) - k;
	    tempX = (AB[0][0]*tx)+(AB[0][1]*ty);
	    tempY = (AB[1][0]*tx)+(AB[1][1]*ty);
	    X = (int) floor(tempX);
	    Y = (int) floor(tempY);
	    rx = tempX - X;
	    ry = tempY - Y;
	    if ((X>=0 && X<=nr1) && (Y>=0 && Y<=nc1))
	    {
		tempA = (1-rx)*(1-ry)*A[X][Y];
		tempB = (1-rx)*ry*A[X][Y+1];
		tempC = (1-ry)*rx*A[X+1][Y];
		Q = tempA + tempB + tempC + (rx*ry*A[X+1][Y+1]);
		tempImage[i][j] = Q;
		DSum = DSum + (1.0/(1.0+abs(Q-C[i][j])));
		count++;
	    }
	    else
		tempImage[i][j] = 0;
	}

    cout << "Number of Rows = " << nrow << endl;
    cout << "Number of Columns = " << ncolumn << endl;
    CopyImage(tempImage,B,nrow,ncolumn);

}

main(int argc, char **argv)
{
    AImage ref,dis,tra;
    int c1, r1, c2, r2, r3, c3, N;
    MPoints m1,m2;

    if (argc != 5)
    {
        cout << "Usage: affine reference distorted ";
        cout << "transformed control_points" << endl;
  	exit(1);
    }

    ReadImage(argv[1],ref,r1,c1);
    ReadImage(argv[2],dis,r2,c2);

    ReadControlPoints(argv[4],m1,m2,N);
    Match_Error(m1,m2,N);

    Match_Images(dis,r2,c2,ref,r1,c1,tra,r3,c3,m2,m1,N);

    WriteImage(argv[3],tra,r3,c3);
}

-- 
--------------------------------------------------------------------------------
Michael Wirth				Email: [email-masked]
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------