Subject: Image Registration Technique - Surface Spline Transformation

From: [email-masked] (Michael Wirth)
Date: 4 Oct 1994 20:52:43 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

typedef double Param[30];

//  Image Registration by Surface Spline Fitting

void Find_Properties(MPoints A, MPoints B, int n, int C, float &a0, float &a1, 
			float &a2, Param D)
{
    matrix25x25 X;
    matrix25x1 Y, Z;
    float s1=0.0,s2=0.0,s3=0.0,s4=0.0,s5=0.0,s6=0.0;
    float r;
    float theta,alpha;
    int i,j;
 
    for (i=0; i<n; i++)
    {
 	X[i][0] = 1;
	X[i][1] = A[i][0];
	X[i][2] = A[i][1];
	for (j=0; j<n; j++)
	{
	    r = pow((A[i][0]-A[j][0]),2) + pow((A[i][1]-A[j][1]),2);
	    if (r == 0)
		X[i][j+3] = 0.0;
	    else
		X[i][j+3] = r * log(r);
	}
    }

    for (i=n; i<n+3; i++)
    {
	X[i][0]=X[i][1]=X[i][2] = 0.0;
	for (j=3; j<n+3; j++)
	{
	    if (i == n)
		X[i][j] = 1;
	    else if (i == n+1)
		X[i][j] = A[j-3][0];
	    else if (i == n+2)
		X[i][j] = A[j-3][1];
	}
    }

    if (C == 1)
    {
    	for (j=0; j<n; j++)
	    Y[j][0] = (double)B[j][0];
	Y[n][0]=Y[n+1][0]=Y[n+2][0]=0.0;
    }
    else if (C == 2)
    {
    	for (j=0; j<n; j++)
	    Y[j][0] = (double)B[j][1];
	Y[n][0]=Y[n+1][0]=Y[n+2][0]=0.0;
    }

    Gaussian_Elimination3(X,Y,Z,n+3);

    a0 = Z[0][0];
    a1 = Z[1][0];
    a2 = Z[2][0];
    for (i=3; i<n+3; i++)
	D[i-3] = Z[i][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;
    Param F1,F2;
    float theta,tempx,tempy,tx,ty,DSum=0.0;
    float a0,a1,a2,b0,b1,b2;
    float nsum;
    int nrow,ncolumn, maxx=0,minx=400,maxy=0,miny=400;
    int i,j,k,l,n=5;
    int count=0,neighbors;
    double sumF1, sumF2,r;
    MPoints w;

    Find_Properties(u,v,cpr,1,a0,a1,a2,F1);
    Find_Properties(u,v,cpr,2,b0,b1,b2,F2);

    for (i=0; i<cpr; i++)
    	cout << F1[i] << " ";
    cout << endl;
    for (i=0; i<cpr; i++)
    	cout << F2[i] << " ";
    cout << endl << endl;

    cout << a0 << " " << a1 << " " << a2 << endl;
    cout << b0 << " " << b1 << " " << b2 << endl;

    for (i=0; i<nr1; i++)
	for (j=0; j<nc1; j++)
	{
	    sumF1 = 0.0;
	    sumF2 = 0.0;
	    for (k=0; k<cpr; k++)
	    {
	    	r = pow((i-u[k][0]),2) + pow((j-u[k][1]),2);
	    	if (r == 0)
		{
		    sumF1 = sumF1 + 0.0;
		    sumF2 = sumF2 + 0.0;
		}
	    	else
		{
		    sumF1 = sumF1 + (F1[k]*r*log(r));
		    sumF2 = sumF2 + (F2[k]*r*log(r));
		}
	    }

	    tempx = a0 + a1*i + a2*j + sumF1;
	    tempy = b0 + b1*i + b2*j + sumF2;

            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<cpr; i++)
    {
	sumF1 = 0.0;
	sumF2 = 0.0;
	for (k=0; k<cpr; k++)
	{
	    r = pow((u[i][0]-u[k][0]),2) + pow((u[i][1]-u[k][1]),2);
	    if (r == 0)
	    {
		sumF1 = sumF1 + 0.0;
		sumF2 = sumF2 + 0.0;
	    }
	    else
	    {
		sumF1 = sumF1 + (F1[k]*r*log(r));
		sumF2 = sumF2 + (F2[k]*r*log(r));
	    }
	}

	tempx = a0 + a1*u[i][0] + a2*u[i][1] + sumF1;
	tempy = b0 + b1*u[i][0] + b2*u[i][1] + sumF2;
	w[i][0] = tempx;
   	w[i][1] = tempy;
    } 

    Match_Error(v,w,cpr);

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

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

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

    for (i=0; i<nr1; i++)
	for (j=0; j<nc1; j++)
	{
	    sumF1 = 0.0;
	    sumF2 = 0.0;
	    for (k=0; k<cpr; k++)
	    {
	    	r = pow((i-u[k][0]),2) + pow((j-u[k][1]),2);
	    	if (r == 0)
		{
		    sumF1 = sumF1 + 0.0;
		    sumF2 = sumF2 + 0.0;
		}
	    	else
		{
		    sumF1 = sumF1 + (F1[k]*r*log(r));
		    sumF2 = sumF2 + (F2[k]*r*log(r));
		}
	    }

	    tempx = a0 + a1*i + a2*j + sumF1;
	    tempy = b0 + b1*i + b2*j + sumF2;

	    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);

    nr2 = nrow;
    nc2 = 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: spline 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                                
--------------------------------------------------------------------------------