Subject: Image Registration Method - Cubic Polynomial

From: [email-masked] (Michael Wirth)
Date: 4 Oct 1994 20:52:02 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"
#include "error.h"

#define PI 3.1415926535979

//  Image Registration by Polynomial Warping
//  Cubic Polynomial - Geometric Distortions

void Find_Properties(MPoints A, MPoints B, int n, int C, float &a0, float &a1, 
			float &a2, float &a3, float &a4, float &a5, float &a6,
			float &a7, float &a8, float &a9)
{
    matrix20x20 X;
    matrix20x1 Y, Z;
    float s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0,s6=0.0,s7=0.0,s8=0.0,s9=0.0,s0=0.0;
    float x=0.0,x2=0.0,x3=0.0,x4=0.0,x5=0.0,x6=0.0;
    float y=0.0,y2=0.0,y3=0.0,y4=0.0,y5=0.0,y6=0.0;
    float sumx=0.0,sumx2=0.0,sumx3=0.0,sumx4=0.0,sumx5=0.0,sumx6=0.0;
    float sumy=0.0,sumy2=0.0,sumy3=0.0,sumy4=0.0,sumy5=0.0,sumy6=0.0;
    float sumxy=0.0,sumx2y2=0.0,sumx3y3=0.0;
    float sumxy2=0.0,sumx2y=0.0,sumxy3=0.0,sumx3y=0.0,sumxy4=0.0,sumx4y=0.0;
    float sumxy5=0.0,sumx5y=0.0,sumx2y3=0.0,sumx3y2=0.0,sumx2y4=0.0,sumx4y2=0.0;
    float theta,alpha;
    int i,j;
 
    for (i=0; i<n; i++)
    {
	x = A[i][0];
	y = A[i][1];
	x2 = pow(A[i][0],2);
	y2 = pow(A[i][1],2);
	x3 = pow(A[i][0],3);
	y3 = pow(A[i][1],3);
	x4 = pow(A[i][0],4);
	y4 = pow(A[i][1],4);
	x5 = pow(A[i][0],5);
	y5 = pow(A[i][1],5);
	x6 = pow(A[i][0],6);
	y6 = pow(A[i][1],6);

	sumx = sumx + x; sumy = sumy + y;
	sumx2 = sumx2 + x2; sumy2 = sumy2 + y2;
	sumx3 = sumx3 + x3; sumy3 = sumy3 + y3;
	sumx4 = sumx4 + x4; sumy4 = sumy4 + y4;
	sumx5 = sumx5 + x5; sumy5 = sumy5 + y5;
	sumx6 = sumx6 + x6; sumy6 = sumy6 + y6;

	sumxy = sumxy + (x*y);
	sumx2y2 = sumx2y2 + (x2*y2);
	sumx3y3 = sumx3y3 + (x3*y3);
	
	sumxy2 = sumxy2 + (x*y2);
	sumx2y = sumx2y + (x2*y);
	sumxy3 = sumxy3 + (x*y3);
	sumx3y = sumx3y + (x3*y);
	sumxy4 = sumxy4 + (x*y4);
	sumx4y = sumx4y + (x4*y);
	sumxy5 = sumxy5 + (x*y5);
	sumx5y = sumx5y + (x5*y);

	sumx2y3 = sumx2y3 + (x2*y3);
	sumx3y2 = sumx3y2 + (x3*y2);
	sumx2y4 = sumx2y4 + (x2*y4);
	sumx4y2 = sumx4y2 + (x4*y2);

    }

    X[0][0]=n;
    X[0][1]=X[1][0]=sumx;
    X[0][2]=X[2][0]=sumy;
    X[0][3]=X[3][0]=X[1][2]=X[2][1]=sumxy;
    X[1][1]=X[0][4]=X[4][0]=sumx2;
    X[2][2]=X[0][5]=X[5][0]=sumy2;
    X[1][3]=X[3][1]=X[2][4]=X[4][2]=X[6][0]=X[0][6]=sumx2y;
    X[7][0]=X[0][7]=X[5][1]=X[1][5]=X[3][2]=X[2][3]=sumxy2;
    X[4][1]=X[1][4]=X[8][0]=X[0][8]=sumx3;
    X[9][0]=X[0][9]=X[2][5]=X[5][2]=sumy3;
    X[8][1]=X[1][8]=X[4][4]=sumx4;
    X[9][2]=X[2][9]=X[5][5]=sumy4;
    X[4][8]=X[8][4]=sumx5;
    X[9][5]=X[5][9]=sumy5;
    X[8][8]=sumx6;
    X[9][9]=sumy6;
    X[3][3]=X[2][6]=X[6][2]=X[1][7]=X[7][1]=X[4][5]=X[5][4]=sumx2y2;
    X[6][1]=X[1][6]=X[8][2]=X[2][8]=X[4][3]=X[3][4]=sumx3y;
    X[9][1]=X[1][9]=X[7][2]=X[2][7]=X[3][5]=X[5][3]=sumxy3;
    X[6][3]=X[3][6]=X[4][7]=X[7][4]=X[8][5]=X[5][8]=sumx3y2;
    X[7][3]=X[3][7]=X[6][5]=X[5][6]=X[9][4]=X[4][9]=sumx2y3;
    X[8][3]=X[3][8]=X[6][4]=X[4][6]=sumx4y;
    X[9][3]=X[3][9]=X[7][5]=X[5][7]=sumxy4;
    X[6][6]=X[8][7]=X[7][8]=sumx4y2;
    X[7][7]=X[9][6]=X[6][9]=sumx2y4;
    X[6][7]=X[7][6]=X[9][8]=X[8][9]=sumx3y3;
    X[8][6]=X[6][8]=sumx5y;
    X[7][9]=X[9][7]=sumxy5;

    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 + B[j][0]*pow(A[j][0],2);
	    s6 = s6 + B[j][0]*pow(A[j][1],2);
	    s7 = s7 + B[j][0]*pow(A[j][0],2)*A[j][1];
	    s8 = s8 + B[j][0]*A[j][0]*pow(A[j][1],2);
	    s9 = s9 + B[j][0]*pow(A[j][0],3);
	    s0 = s0 + B[j][0]*pow(A[j][1],3);
    	}
    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 + B[j][1]*pow(A[j][0],2);
	    s6 = s6 + B[j][1]*pow(A[j][1],2);
	    s7 = s7 + B[j][1]*pow(A[j][0],2)*A[j][1];
	    s8 = s8 + B[j][1]*A[j][0]*pow(A[j][1],2);
	    s9 = s9 + B[j][1]*pow(A[j][0],3);
	    s0 = s0 + B[j][1]*pow(A[j][1],3);
    	}

    Y[0][0]=s1; Y[1][0]=s2; Y[2][0]=s3; Y[3][0]=s4; Y[4][0]=s5; 
    Y[5][0]=s6; Y[6][0]=s7; Y[7][0]=s8; Y[8][0]=s9; Y[9][0]=s0;

    Gaussian_Elimination2(X,Y,Z,10);

    a0 = Z[0][0];
    a1 = Z[1][0];
    a2 = Z[2][0];
    a3 = Z[3][0];
    a4 = Z[4][0];
    a5 = Z[5][0];
    a6 = Z[6][0];
    a7 = Z[7][0];
    a8 = Z[8][0];
    a9 = Z[9][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,newImage;
    MPoints w;
    float theta,tempx,tempy,tx,ty,DSum=0.0;
    float a0,a1,a2,a3,a4,a5,a6,a7,a8,a9;
    float b0,b1,b2,b3,b4,b5,b6,b7,b8,b9;
    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,a6,a7,a8,a9);
    Find_Properties(u,v,cpr,2,b0,b1,b2,b3,b4,b5,b6,b7,b8,b9);

cout << a0 << " " << a1 << " " << a2 << " " << a3 << " " << a4 << " " << a5;
cout << a6 << " " << a7 << " " << a8 << " " << a9 << endl;
cout << b0 << " " << b1 << " " << b2 << " " << b3 << " " << b4 << " " << b5;
cout << b6 << " " << b7 << " " << b8 << " " << b9 << 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) +
a6*pow(i,2)*j + a7*i*pow(j,2) + a8*pow(i,3) + a9*pow(j,3);
	    tempy = b0 + b1*i + b2*j + b3*i*j + b4*pow(i,2) + b5*pow(j,2) +
b6*pow(i,2)*j + b7*i*pow(j,2) + b8*pow(i,3) + b9*pow(j,3);

            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) +
a6*pow(i,2)*j + a7*i*pow(j,2) + a8*pow(i,3) + a9*pow(j,3);
	    tempy = b0 + b1*i + b2*j + b3*i*j + b4*pow(i,2) + b5*pow(j,2) +
b6*pow(i,2)*j + b7*i*pow(j,2) + b8*pow(i,3) + b9*pow(j,3);
	    tx = (tempx+abs(minx));
	    ty = (tempy+abs(miny));
	}

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

	for (j=0; j<cpr; 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) + a6*pow(u[j][0],2)*u[j][1] +
a7*u[j][0]*pow(u[j][1],2) + a8*pow(u[j][0],3) + a9*pow(u[j][1],3);
	    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) + b6*pow(u[j][0],2)*u[j][1] +
b7*u[j][0]*pow(u[j][1],2) + b8*pow(u[j][0],3) + b9*pow(u[j][1],3);
cout << v[j][0] << "\t" << v[j][1] << "\t" << tempx << "\t" << tempy << endl;   
	w[j][0] = tempx;
	w[j][1] = tempy;
 	} 

    Match_Error(v,w,cpr);
    
    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++)
	{
	    newImage[i][j] = 0;
	    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) +
a6*pow(i,2)*j + a7*i*pow(j,2) + a8*pow(i,3) + a9*pow(j,3);
	    tempy = b0 + b1*i + b2*j + b3*i*j + b4*pow(i,2) + b5*pow(j,2) +
b6*pow(i,2)*j + b7*i*pow(j,2) + b8*pow(i,3) + b9*pow(j,3);
	    tx = (tempx+abs(minx));
	    ty = (tempy+abs(miny));
  	    newImage[(int)tx][(int)ty] = A[i][j];
	    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];
			}
		}
		if (neighbors == 0)
		    newImage[i][j] = 0;
	    	else
		    newImage[i][j] = (int)(nsum/neighbors);
	    }
 	}
    }

CopyImage(newImage,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_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                                
--------------------------------------------------------------------------------