Subject: Image Registration Method - Quadratic Polynomial

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

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

#define PI 3.1415926535979

//  Image Registration by Polynomial Warping
//  Quadratic Polynomial - General Affine Distortions

void Find_Properties(MPoints A, MPoints B, int n, int C, float &a0, float &a1, 
			float &a2, float &a3, float &a4, float &a5)
{
    matrix6x6 X;
    float Y[6][1], W[6][1];
    float s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0,s6=0.0,temp;
    float x=0.0,y=0.0,xsq=0.0,ysq=0.0,xy=0.0,xsqy=0.0,xysq=0.0;
    float xsqysq=0.0,xcb=0.0,ycb=0.0,xq=0.0,yq=0.0,xycb=0.0,xcby=0.0;
    float theta,alpha;
    int i,j;
 
    for (i=0; i<n; i++)
    {
	x = x + A[i][0];
	y = y + A[i][1];
	xsq = xsq + pow(A[i][0],2);
	ysq = ysq + pow(A[i][1],2);
	xy = xy + (A[i][0]*A[i][1]);
	xsqy = xsqy + (pow(A[i][0],2)*A[i][1]);
	xysq = xysq + (pow(A[i][1],2)*A[i][0]);
	xsqysq = xsqysq + (pow(A[i][0],2)*pow(A[i][1],2));
	xcb = xcb + pow(A[i][0],3);
	ycb = ycb + pow(A[i][1],3);
	xq = xq + pow(A[i][0],4);
	yq = yq + pow(A[i][1],4);
	xycb = xycb + (pow(A[i][1],3)*A[i][0]);
	xcby = xcby + (pow(A[i][0],3)*A[i][1]);
    }

    X[0][0]=n;
    X[0][1]=X[1][0]=x;
    X[0][2]=X[2][0]=y;
    X[0][3]=X[3][0]=X[1][2]=X[2][1]=xy;
    X[0][4]=X[4][0]=X[1][1]=xsq;
    X[0][5]=X[5][0]=X[2][2]=ysq;
    X[1][3]=X[2][4]=X[3][1]=X[4][2]=xsqy;
    X[1][5]=X[2][3]=X[5][1]=X[3][2]=xysq;
    X[1][4]=X[4][1]=xcb;
    X[2][5]=X[5][2]=ycb;
    X[4][4]=xq;
    X[5][5]=yq;
    X[3][3]=X[4][5]=X[5][4]=xsqysq;
    X[3][4]=X[4][3]=xcby;
    X[3][5]=X[5][3]=xycb;

    if (C == 1)
    	for (j=0; j<n; j++)
    	{
	    s1 = s1 + B[j][0];
	    s2 = s2 + (A[j][0]*B[j][0]);
	    s3 = s3 + (A[j][1]*B[j][0]);
	    s4 = s4 + (A[j][0]*A[j][1]*B[j][0]);
	    s5 = s5 + pow(A[j][0],2)*B[j][0];
	    s6 = s6 + pow(A[j][1],2)*B[j][0];
    	}
    else if (C == 2)
    	for (j=0; j<n; j++)
    	{
	    s1 = s1 + B[j][1];
	    s2 = s2 + (A[j][0]*B[j][1]);
	    s3 = s3 + (A[j][1]*B[j][1]);
	    s4 = s4 + (A[j][0]*A[j][1]*B[j][1]);
	    s5 = s5 + pow(A[j][0],2)*B[j][1];
	    s6 = s6 + pow(A[j][1],2)*B[j][1];
    	}

    Y[0][0]=s1; Y[1][0]=s2; Y[2][0]=s3; Y[3][0]=s4; Y[4][0]=s5; Y[5][0]=s6;

    InverseMatrix6x6(X,6);

    for (i=0; i<6; i++)
    {
	temp = 0.0;
	for (j=0; j<6; j++)
	    temp = temp + (X[i][j]*Y[j][0]);
	W[i][0] = temp;
    }
    a0 = W[0][0];
    a1 = W[1][0];
    a2 = W[2][0];
    a3 = W[3][0];
    a4 = W[4][0];
    a5 = W[5][0];

}

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;
    float cat[300][300][2];
    float theta,tempx,tempy,tx,ty,DSum=0.0;
    float a0,a1,a2,a3,a4,a5,b0,b1,b2,b3,b4,b5;
    float nsum;
    int nrow,ncolumn, maxx=0,minx=400,maxy=0,miny=400;
    int i,j,k,l,count=0,neighbors;

    Find_Properties(u,v,cpr,1,a0,a1,a2,a3,a4,a5);
    Find_Properties(u,v,cpr,2,b0,b1,b2,b3,b4,b5);

cout << a0 << " " << a1 << " " << a2 << " " << a3 << " " << a4 << " " << a5 << endl;
cout << b0 << " " << b1 << " " << b2 << " " << b3 << " " << b4 << " " << b5 << endl;

    for (i=0; i<nr1; i++)
	for (j=0; j<nc1; j++)
	{
	    tempx = a0 + a1*i + a2*j + a3*i*j + a4*pow(i,2) + a5*pow(j,2);
	    tempy = b0 + b1*i + b2*j + b3*i*j + b4*pow(i,2) + b5*pow(j,2);

            if (tempx < minx)
		minx = (int)tempx;
	    if (tempx > maxx)
		maxx = (int)tempx;
            if (tempy < miny)
		miny = (int)tempy;
	    if (tempy > maxy)
		maxy = (int)tempy;
	} 
    
    for (i=0; i<nr1; i++)
	for (j=0; j<nc1; j++)
	{
	    tempx = a0 + a1*i + a2*j + a3*i*j + a4*pow(i,2) + a5*pow(j,2);
	    tempy = b0 + b1*i + b2*j + b3*i*j + b4*pow(i,2) + b5*pow(j,2);
	    tx = (tempx+abs(minx));
	    ty = (tempy+abs(miny));
	    cat[(int)tx][(int)ty][0] = i;
	    cat[(int)tx][(int)ty][1] = j;
	}

	cout << "Transformed Control Points: " << endl;

	for (j=0; j<5; j++)
	{
	    tempx = a0 + a1*u[j][0] + a2*u[j][1] + a3*u[j][0]*u[j][1] + 				a4*pow(u[j][0],2) + a5*pow(u[j][1],2);
	    tempy = b0 + b1*u[j][0] + b2*u[j][1] + b3*u[j][0]*u[j][1] + b4*pow(u[j][0],2) + b5*pow(u[j][1],2);
 	} 

cout << minx << " " << maxx << " " << miny << " " << maxy << endl;

    nrow = maxx-minx;
    ncolumn = maxy-miny;
    nr2 = nrow;
    nc2 = ncolumn;

    for (i=0; i<nrow; i++)
	for (j=0; j<ncolumn; j++)
	    tempImage[i][j] = -1;

    for (i=0; i<nr1; i++)
	for (j=0; j<nc1; j++)
	{
	    tempx = a0 + a1*i + a2*j + a3*i*j + a4*pow(i,2) + a5*pow(j,2);
	    tempy = b0 + b1*i + b2*j + b3*i*j + b4*pow(i,2) + b5*pow(j,2);
	    tx = (tempx+abs(minx));
	    ty = (tempy+abs(miny));
	    tempImage[(int)tx][(int)ty] = A[i][j];
	}
	
    for (i=0; i<nrow; i++)
	for (j=0; j<ncolumn; j++)
	{
	    neighbors = 0;
	    nsum = 0.0;
	    if (tempImage[i][j] == -1)
	    {
	     	for (k=-1; k<=1; k++)
		    for (l=-1; l<=1; l++)
			if (tempImage[i+k][j+l] != -1)
			{
			    neighbors++;
			    nsum = nsum + tempImage[i+k][j+l];
			}
	    	tempImage[i][j] = (int)(nsum/neighbors);
	    }
	}

    Box_Filter(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: poly reference distorted transformed ";
        cout << "control_points" << endl;
  	exit(1);
    }

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

    ReadControlPoints(argv[4],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                                
--------------------------------------------------------------------------------