From mwirth@cs.umanitoba.ca Mon Oct 10 14:03:30 1994
X-VM-Message-Order:
	(1 6 2 3 4 5 7 8)
X-VM-Summary-Format: "%n %*%a %-17.17F %-3.3m %2d %4l/%-5c %I\"%s\"\n"
X-VM-Labels: nil
X-VM-VHeader: ("Resent-" "From:" "Sender:" "To:" "Apparently-To:" "Cc:" "Subject:" "Date:") nil
X-VM-Bookmark: 8
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["2235" "" " 4" "October" "1994" "20:43:44" "GMT" "Michael Wirth" "mwirth@cs.umanitoba.ca" "<36seq0$suj@canopus.cc.umanitoba.ca>" "74" "Image Registration Transformation Methods (File 1)" "^From:" nil nil "10" "1994100420:43:44" "Image Registration Transformation Methods (File 1)" (number " " mark "     Michael Wirth     Oct  4   74/2235  " thread-indent "\"Image Registration Transformation Methods (File 1)\"\n") nil]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!concert!gatech!swiss.ans.net!solaris.cc.vt.edu!MathWorks.Com!europa.eng.gtefsd.com!howland.reston.ans.net!math.ohio-state.edu!scipio.cyberstore.ca!vanbc.wimsey.com!unixg.ubc.ca!quartz.ucs.ualberta.ca!tribune.usask.ca!canopus.cc.umanitoba.ca!cs.umanitoba.ca!mwirth
Newsgroups: sci.image.processing
Organization: Computer Science, University of Manitoba, Winnipeg, Canada
Lines: 74
Distribution: world
Message-ID: <36seq0$suj@canopus.cc.umanitoba.ca>
NNTP-Posting-Host: silver.cs.umanitoba.ca
From: mwirth@cs.umanitoba.ca (Michael Wirth)
Subject: Image Registration Transformation Methods (File 1)
Date: 4 Oct 1994 20:43:44 GMT

There are five files: The first file contains 5 header files, which must be
separated and saved as the following files:

	imagematrix.h
	inverse.h
	average_filters.h
	gaussian.h
	error.h

The next four files contain the transformation methods:

	affine.c 	- Affine Transformation
	polyquad.c	- Polynomial Transformation - Quadratic
	polycubic.c	- Polynomial Transformation - Cubic
	spline.c	- Spline Transformation

Each file can be compiled using the standard Sun C++ compiler. Input
to each of the programmes is in the following format:

	a.out reference_file distorted_file transformed_file control_points

where:

	reference_file = the reference image
	distorted_file = the image to be transformed
	transformed_file = the transformed image
	control_points = a file containing the control point pairs, in the
			 following form:

			 U1 V1 X1 Y1
			 U2 V2 X2 Y2
			 .  .  .  .
			 .  .  .  .
			 .  .  .  . 
			 Un Vn Xn Yn

			 U and V are the x and y coordinates of the control
			 point pairs in the reference image and X and Y are
			 the corresponding x and y coordinates of the control
			 point pairs in the distorted image. n is the number
			 of control-point pairs.

	The images are treated as discrete 2D arrays, and the files are of the
 	PGM form:

		columns rows
		- - - - - - - - - - - - - - - - - - - - - 
		- - - - - - - - - - - - - - - - - - - - - 
		- - - - - - - - - - - - - - - - - - - - - 
		.......
		.......
		.......
		- - - - - - - - - - - - - - - - - - - - - 


	columns = the number of columns
	rows	= the number of rows
	-	= a grayscale value (0-255)

	The programme will produce a transformed image of the form shown
	above and also the Mean-Square-Root-Error estimate for each
	control-point pair.


	(c) September 1994, Michael Wirth

-- 
--------------------------------------------------------------------------------
Michael Wirth				Email: mwirth@silver.cs.umanitoba.ca
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------

From mwirth@cs.umanitoba.ca Mon Oct 10 14:04:00 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["13972" "" " 5" "October" "1994" "18:22:40" "GMT" "Michael Wirth" "mwirth@cs.umanitoba.ca" "<36uqtg$n2q@canopus.cc.umanitoba.ca>" "681" "Image Registration Header Files" "^From:" nil nil "10" "1994100518:22:40" "Image Registration Header Files" (number " " mark "     Michael Wirth     Oct  5  681/13972 " thread-indent "\"Image Registration Header Files\"\n") nil]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!concert!inxs.ncren.net!taco.cc.ncsu.edu!gatech!howland.reston.ans.net!europa.eng.gtefsd.com!library.ucla.edu!news.mic.ucla.edu!unixg.ubc.ca!quartz.ucs.ualberta.ca!tribune.usask.ca!canopus.cc.umanitoba.ca!cs.umanitoba.ca!mwirth
Newsgroups: sci.image.processing
Organization: Computer Science, University of Manitoba, Winnipeg, Canada
Lines: 681
Distribution: world
Message-ID: <36uqtg$n2q@canopus.cc.umanitoba.ca>
NNTP-Posting-Host: silver.cs.umanitoba.ca
From: mwirth@cs.umanitoba.ca (Michael Wirth)
Subject: Image Registration Header Files
Date: 5 Oct 1994 18:22:40 GMT

Oops... seems the header files only posted locally, so here it is again.

======================================================================
Imagematrix.h
======================================================================

#include <stdio.h>
#include <fstream.h>
#define MAX(a, b) (a > b ? (a) : (b))

typedef int AImage[500][500];
typedef double Afilter[150][150];
typedef int Contour[600][2];
typedef int MPoints[50][2];

void ViewImageMatrix(AImage A, int nrow, int ncolumn)
{
    for (int i=0; i<nrow; i++)
    {
	for (int j=0; j<ncolumn; j++)
	    cout << A[i][j] << " ";
	cout << endl;
    }
}

void ReadImage(char *filename, AImage A, int &row, int &column)
{
    ifstream fp;

    fp.open(filename,ios::in);
    fp >> column >> row;

    for (int i=0; i<row; i++)
	for (int j=0; j<column; j++)
	    fp >> A[i][j];

    fp.close();
}

void CopyImage(AImage A, AImage B, int row, int column)
{
    for (int i=0; i<row; i++)
	for (int j=0; j<column; j++)
	    B[i][j] = A[i][j];
}

void WriteImage(char *filename, AImage A, int row, int column)
{
    ofstream fq;

    fq.open(filename,ios::out);
    fq << column << " " << row << endl;

    for (int i=0; i<row; i++)
    {
	for (int j=0; j<column; j++)
	    fq << A[i][j] << " ";
	fq << endl;
    }
    fq.close(); 
}

void WriteGrid(char *filename, AImage A, int row, int column)
{
    ofstream fq;

    fq.open(filename,ios::out);

    for (int i=0; i<row; i++)
    {
	for (int j=0; j<column; j++)
	    fq << A[i][j] << " ";
	fq << endl;
    }
    fq.close(); 
}

void ReadContour(char *filename, Contour C, int &N)
{
    ifstream fp;
    int temp;

    fp.open(filename,ios::in);

    fp >> N;
    for (int i=0; i<N; i++)
	fp >> temp >> C[i][0] >> C[i][1];
    fp.close(); 
}

void WriteContour(char *filename, Contour C, int N)
{
    ofstream fq;

    fq.open(filename,ios::out);

    fq << N << endl;
    for (int i=0; i<N; i++)
	fq << i+1 << " " << C[i][0] << " " << C[i][1] << endl;
    fq.close(); 
}

void ReadControlPoints(char *filename, MPoints M1, MPoints M2, int &N)
{
    ifstream fp;
    int i;

    fp.open(filename,ios::in);
    
    i = 0;
    do
    {
	fp >> M1[i][0] >> M1[i][1] >> M2[i][0] >> M2[i][1];
        i++;
    } while (!fp.eof());
    fp.close(); 
    N = i - 1;
}

void WritePGM(char *filename, AImage A, int r, int c)
{
    int i,j,k,max=0,row,column;
    AImage B;
    ifstream fp;
    ofstream fq;
 
    WriteImage("image.tmp",A,r,c);

    fp.open("image.tmp",ios::in);

    fp >> row >> column;

    double shade[400][400];
  
    for (i=0; i<row; i++)
    	for (j=0; j<column; j++)
    	{
     	    fp >> shade[i][j];
      	    B[i][j] = (int) shade[i][j];
      	    max = MAX(max, B[i][j]);
    	}
    fp.close();

    fq.open(filename,ios::out);

    fq << "P2\n" << row << "  " << column << "  " << max << "\n";
    fq << "# gridtopgm  M.W. 1994\n";
  
    for (i=0; i<row; i++)
    {
    	for (j=0; j<column / 10; j++)
    	{
      	    for ( k = 0; k < 10; k++)
	    	fq << B[i][j * 10 + k] << " ";
      	    fq << "\n";
    	}
    	for (k = 0; k < column % 10; k++)
      	    fq << B[i][j * 10 + k] << " ";
    	fq << "\n\n";
    }

    fq.close();
}

==================================================================
inverse.h
==================================================================

/* shply.h
 * 
 * Matrix inversion include file.
 * Inverts a matrix in-place.
 *
 * Using the Shipley-Coleman technique with Kron reduction.
 *
 * Russ Patterson (tva!tprwp.chattd@mhs.attmail.com)
 * 
 * A is square matrix.
 * N is dimension of A.
 */


typedef double matrix8x8[8][8];
typedef double matrix6x6[6][6];
typedef double matrix4x4[4][4];
typedef double matrix2x2[2][2];
typedef double matrix8x1[8][1];

void Gaussian_Elimination(matrix8x8 A, matrix8x1 B, matrix8x1 C)
{
    int i,j,k,p,FOUND=0,n=8;
    double temp;
    matrix8x8 M;

    for (i=0; i<n; i++)
    {
	p = i;
	while (p>=i && p<n && !FOUND)
	{
	    if (A[p][i] != 0)
		FOUND = 1;
	    else
		p++;
	}

	if (p != i)
	{
	    for (j=0; j<n; j++)
	    {
		temp = A[i][j];
		A[i][j] = A[p][j];
		A[p][j] = temp;
	    }
	    temp = B[i][0];
	    B[i][0] = B[p][0];
	    B[p][0] = temp;
	}

	for (k=i+1; k<n; k++)
	{
	    M[k][i] = A[k][i]/A[i][i];
	    for (j=0; j<n; j++)
		A[k][j] = A[k][j] - M[k][i]*A[i][j];
	    B[k][0] = B[k][0] - M[k][i]*B[i][0];
	}
    }

    if (A[n-1][n-1] == 0.0)
	cout << "no unique solution exists" << endl;
	
    C[n-1][0] = B[n-1][0]/A[n-1][n-1];

    for (i=n-2; i>=0; i--)
    {
	temp = 0.0;
	for (j=i+1; j<n; j++)
	    temp = temp + (A[i][j]*C[j][0]);
	C[i][0] = (B[i][0] - temp)/A[i][i];
    }

    for (i=0; i<n; i++)
	cout << C[i][0] << endl;  
}


void InverseMatrix8x8(matrix8x8 A, int N)
{
	int i,j,p,n;     // N is dimension
	n=N;

	// Kron reduction, begin with pivot axis #1
    	for (p=0; p<N; p++)
	{
	    for (i=0; i<n; i++)
 		if(i != p)
		    for (j=0; j<n; j++)
		    {
			if(j != p)
			    A[i][j] = (A[i][j] - (A[i][p] * A[p][j] /
						A[p][p]) );
		    }
	    A[p][p] = (-1.0 / A[p][p] );
	    i=p;
  	    for (j=0; j<n; j++)
	    	if(i != j)
	    	{
		    A[i][j] = A[p][p] * A[i][j] ;
		    A[j][i] = A[p][p] * A[j][i] ;
	    	}
	}
// result of above code is -inv(A)
// following code negates every element in the -inv(A) matrix.
		
	for (i=0; i<N; i++)
	      for (j=0; j<n; j++)
		      A[i][j] = -1.0 * A[i][j] ;
}

void InverseMatrix6x6(matrix6x6 A, int N)
{
	int i,j,p,n;     // N is dimension
	n=N;

	// Kron reduction, begin with pivot axis #1
    	for (p=0; p<N; p++)
	{
	    for (i=0; i<n; i++)
 		if(i != p)
		    for (j=0; j<n; j++)
		    {
			if(j != p)
			    A[i][j] = (A[i][j] - (A[i][p] * A[p][j] /
						A[p][p]) );
		    }
	    A[p][p] = (-1.0 / A[p][p] );
	    i=p;
  	    for (j=0; j<n; j++)
	    	if(i != j)
	    	{
		    A[i][j] = A[p][p] * A[i][j] ;
		    A[j][i] = A[p][p] * A[j][i] ;
	    	}
	}
// result of above code is -inv(A)
// following code negates every element in the -inv(A) matrix.
		
	for (i=0; i<N; i++)
	      for (j=0; j<n; j++)
		      A[i][j] = -1.0 * A[i][j] ;
}

void InverseMatrix(matrix4x4 A, int N)
{
	int i,j,p,n;     // N is dimension
	n=N;

	// Kron reduction, begin with pivot axis #1
    	for (p=0; p<N; p++)
	{
	    for (i=0; i<n; i++)
 		if(i != p)
		    for (j=0; j<n; j++)
		    {
			if(j != p)
			    A[i][j] = (A[i][j] - (A[i][p] * A[p][j] /
						A[p][p]) );
		    }
	    A[p][p] = (-1.0 / A[p][p] );
	    i=p;
  	    for (j=0; j<n; j++)
	    	if(i != j)
	    	{
		    A[i][j] = A[p][p] * A[i][j] ;
		    A[j][i] = A[p][p] * A[j][i] ;
	    	}
	}
// result of above code is -inv(A)
// following code negates every element in the -inv(A) matrix.
		
	for (i=0; i<N; i++)
	      for (j=0; j<n; j++)
		      A[i][j] = -1.0 * A[i][j] ;
}

void InverseMatrix2x2(matrix2x2 A, int N)
{
	int i,j,p,n;     // N is dimension
	n=N;

	// Kron reduction, begin with pivot axis #1
    	for (p=0; p<N; p++)
	{
	    for (i=0; i<n; i++)
 		if(i != p)
		    for (j=0; j<n; j++)
		    {
			if(j != p)
			    A[i][j] = (A[i][j] - (A[i][p] * A[p][j] /
						A[p][p]) );
		    }
	    A[p][p] = (-1.0 / A[p][p] );
	    i=p;
  	    for (j=0; j<n; j++)
	    	if(i != j)
	    	{
		    A[i][j] = A[p][p] * A[i][j] ;
		    A[j][i] = A[p][p] * A[j][i] ;
	    	}
	}
// result of above code is -inv(A)
// following code negates every element in the -inv(A) matrix.
		
	for (i=0; i<N; i++)
	      for (j=0; j<n; j++)
		      A[i][j] = -1.0 * A[i][j] ;
}

==================================================================
Average_filters.h
==================================================================

#include <stdio.h>
#include <stdlib.h>
#include <fstream.h>

#define PI 3.1415926535979
	
void Bartlett_Filter(AImage data, AImage newImage, int nrow, int ncolumn)
{
    static int filter[3][3] = { {1,2,1},{2,4,2},{1,2,1} };
    float avg, temp;

    for (int x=0; x<nrow; x++)
	for (int y=0; y<ncolumn; y++)
	{
	    temp = 0.0;
	    avg = 0.0;
	    for (int i=-1; i<=1; i++) {
		for (int j=-1; j<=1; j++)
		{
		    if (i+x>=0 && j+y<ncolumn && i+x<nrow && j+y>=0)
		    {
			temp += data[i+x][j+y]*filter[i+1][j+1];
			avg += filter[i+1][j+1];
		    }
		}
	    }
	    newImage[x][y] = (int)(temp/avg);
	}
}

void Box_Filter(AImage data, AImage newImage, int nrow, int ncolumn)
{
    static int filter[3][3] = { {1,1,1},{1,1,1},{1,1,1} };
    float temp;

    for (int x=0; x<nrow; x++)
	for (int y=0; y<ncolumn; y++)
	{
	    temp = 0.0;
	    for (int i=-1; i<=1; i++) {
		for (int j=-1; j<=1; j++)
		{
		    if (i+x>=0 && j+y<ncolumn && i+x<nrow && j+y>=0)
			temp += data[i+x][j+y]*filter[i+1][j+1];
		}
	    }
	    newImage[x][y] = (int)temp/9;
	}
}

void Add_Snow_Noise(AImage data, int nrow, int ncolumn)
{
    double randtemp;
    for (int i=0; i<nrow; i++)
	for (int j=0; j<ncolumn; j++)
	{
	    randtemp = drand48()*10;
	    if (randtemp < 0.5)
		data[i][j] = 128;
	}
}

void Add_Random_Noise(AImage data, int nrow, int ncolumn, int noise)
{
    double randtemp;
    for (int i=0; i<nrow; i++)
	for (int j=0; j<ncolumn; j++)
	{
	    randtemp = drand48()*10;
	    if (randtemp <= 5.0)
	    {
		if (data[i][j] <= 255-noise)
	            data[i][j] += noise;
	    }
 	    else {
		if (data[i][j] >= noise)
		    data[i][j] -= noise;
	    }
	}
}

void ReSample(AImage A, AImage B, int nrow, int ncolumn)
{
    int i,j;
    float maxvalue=0.0, minvalue=500.0;
    
    for (i=0; i<nrow; i++)
	for (j=0; j<ncolumn; j++)
	{
	    if (A[i][j] > maxvalue)
		maxvalue = A[i][j];
	    if (A[i][j] < minvalue)
		minvalue = A[i][j];
	}

    for (i=0; i<nrow; i++)
	for (j=0; j<ncolumn; j++)
	{
            B[i][j] = (int)(((A[i][j]-minvalue)/maxvalue)*255.0); 
	}
}

void HighPass_Filter(AImage data, AImage newImage, int nrow, int ncolumn)
{
    static int Hfilter[3][3] = { {-1,-1,-1},{-1,9,-1},{-1,-1,-1} };
    float avg, temp;
    AImage tempImage;

    for (int x=0; x<nrow; x++)
	for (int y=0; y<ncolumn; y++)
	{
	    temp = 0.0;
	    avg = 0.0;
	    for (int i=-1; i<=1; i++)
		for (int j=-1; j<=1; j++)
		{
		    if (i+x>=0 && j+y<ncolumn && i+x<nrow && j+y>=0)
		    {
			temp += data[i+x][j+y]*Hfilter[i+1][j+1];
			avg += Hfilter[i+1][j+1];
		    }
		}

	    if (temp > 0)
	        tempImage[x][y] = (int)temp;
	    else
	        tempImage[x][y] = 0;
	}
    ReSample(tempImage,newImage,nrow,ncolumn);
}

==================================================================
gaussian.h
==================================================================

typedef double matrix20x20[20][20];
typedef double matrix20x1[20][1];
typedef double matrix25x25[25][25];
typedef double matrix25x1[25][1];

void Gaussian_Elimination2(matrix20x20 A, matrix20x1 B, matrix20x1 C, int n)
{
    int i,j,k,p,FOUND=0;
    double temp;
    matrix20x20 M;

    for (i=0; i<n; i++)
    {
	p = i;
	while (p>=i && p<n && !FOUND)
	{
	    if (A[p][i] != 0)
		FOUND = 1;
	    else
		p++;
	}

	if (p != i)
	{
	    for (j=0; j<n; j++)
	    {
		temp = A[i][j];
		A[i][j] = A[p][j];
		A[p][j] = temp;
	    }
	    temp = B[i][0];
	    B[i][0] = B[p][0];
	    B[p][0] = temp;
	}

	for (k=i+1; k<n; k++)
	{
	    M[k][i] = A[k][i]/A[i][i];
	    for (j=0; j<n; j++)
		A[k][j] = A[k][j] - M[k][i]*A[i][j];
	    B[k][0] = B[k][0] - M[k][i]*B[i][0];
	}
    }

    if (A[n-1][n-1] == 0.0)
	cout << "no unique solution exists" << endl;
    
    C[n-1][0] = B[n-1][0]/A[n-1][n-1];

    for (i=n-2; i>=0; i--)
    {
	temp = 0.0;
	for (j=i+1; j<n; j++)
	    temp = temp + (A[i][j]*C[j][0]);
	C[i][0] = (B[i][0] - temp)/A[i][i];
    }
}

void Gaussian_Elimination3(matrix25x25 A, matrix25x1 B, matrix25x1 C, int n)
{
    int i,j,k,p,FOUND=0;
    double temp;
    matrix25x25 M;

    for (i=0; i<n; i++)
    {
	p = i;
	while (p>=i && p<n && !FOUND)
	{
	    if (A[p][i] != 0)
		FOUND = 1;
	    else
		p++;
	}

	if (p != i)
	{
	    for (j=0; j<n; j++)
	    {
		temp = A[i][j];
		A[i][j] = A[p][j];
		A[p][j] = temp;
	    }
	    temp = B[i][0];
	    B[i][0] = B[p][0];
	    B[p][0] = temp;
	}

	for (k=i+1; k<n; k++)
	{
	    M[k][i] = A[k][i]/A[i][i];
	    for (j=0; j<n; j++)
		A[k][j] = A[k][j] - M[k][i]*A[i][j];
	    B[k][0] = B[k][0] - M[k][i]*B[i][0];
	}
    }

    if (A[n-1][n-1] == 0.0)
	cout << "no unique solution exists" << endl;
	
    C[n-1][0] = B[n-1][0]/A[n-1][n-1];

    for (i=n-2; i>=0; i--)
    {
	temp = 0.0;
	for (j=i+1; j<n; j++)
	    temp = temp + (A[i][j]*C[j][0]);
	C[i][0] = (B[i][0] - temp)/A[i][i];
    }
}

==================================================================
error.h
==================================================================

#include <stdio.h>
#include <stdlib.h>
#include <fstream.h>

void Match_Error(MPoints A, MPoints B, int n)
{
    int i;
    float error,temp1,temp2,msre=0.0;
   
    cout << "Reference" << "\t\t" << "Distorted" << endl;
 
    for (i=0; i<n; i++)
    {
   	temp1 = A[i][0] - B[i][0];
	temp2 = A[i][1] - B[i][1];
	error = sqrt(pow(temp1,2)+pow(temp2,2));
        msre = msre + error;
	
	cout << A[i][0] << "," << A[i][1] << "\t=>\t" << B[i][0] << "," << B[i][1];
        cout << "\t" << error << endl; 
    }
    cout << "MSRE = " << msre/n << endl;
}

==================================================================
==================================================================

Please direct any comments to the following address, as I am returning to
Australia, and as yet do not have an email address:

	Michael Wirth
	7 Pipit Parade 
	Burleigh Waters
	Queensland 4220
	AUSTRALIA	

Copies of my Master's thesis can be obtained from Dafoe Library at the 
University of Manitoba.

regards

Michael		
-- 
--------------------------------------------------------------------------------
Michael Wirth				Email: mwirth@silver.cs.umanitoba.ca
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------

From mwirth@cs.umanitoba.ca Mon Oct 10 14:05:03 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["4477" "" " 4" "October" "1994" "20:48:44" "GMT" "Michael Wirth" "mwirth@cs.umanitoba.ca" "<36sf3c$suj@canopus.cc.umanitoba.ca>" "183" "Image Registration - Affine Transformation" "^From:" nil nil "10" "1994100420:48:44" "Image Registration - Affine Transformation" (number " " mark "     Michael Wirth     Oct  4  183/4477  " thread-indent "\"Image Registration - Affine Transformation\"\n") nil]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!concert!news.duke.edu!MathWorks.Com!europa.eng.gtefsd.com!library.ucla.edu!news.mic.ucla.edu!unixg.ubc.ca!quartz.ucs.ualberta.ca!tribune.usask.ca!canopus.cc.umanitoba.ca!cs.umanitoba.ca!mwirth
Newsgroups: sci.image.processing
Organization: Computer Science, University of Manitoba, Winnipeg, Canada
Lines: 183
Distribution: world
Message-ID: <36sf3c$suj@canopus.cc.umanitoba.ca>
NNTP-Posting-Host: silver.cs.umanitoba.ca
From: mwirth@cs.umanitoba.ca (Michael Wirth)
Subject: Image Registration - Affine Transformation
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: mwirth@silver.cs.umanitoba.ca
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------

From mwirth@cs.umanitoba.ca Mon Oct 10 14:05:09 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["6029" "" " 4" "October" "1994" "20:51:13" "GMT" "Michael Wirth" "mwirth@cs.umanitoba.ca" "<36sf81$suj@canopus.cc.umanitoba.ca>" "229" "Image Registration Method - Quadratic Polynomial" "^From:" nil nil "10" "1994100420:51:13" "Image Registration Method - Quadratic Polynomial" (number " " mark "     Michael Wirth     Oct  4  229/6029  " thread-indent "\"Image Registration Method - Quadratic Polynomial\"\n") nil]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!concert!news.duke.edu!MathWorks.Com!europa.eng.gtefsd.com!library.ucla.edu!news.mic.ucla.edu!unixg.ubc.ca!quartz.ucs.ualberta.ca!tribune.usask.ca!canopus.cc.umanitoba.ca!cs.umanitoba.ca!mwirth
Newsgroups: sci.image.processing
Organization: Computer Science, University of Manitoba, Winnipeg, Canada
Lines: 229
Distribution: world
Message-ID: <36sf81$suj@canopus.cc.umanitoba.ca>
NNTP-Posting-Host: silver.cs.umanitoba.ca
From: mwirth@cs.umanitoba.ca (Michael Wirth)
Subject: Image Registration Method - Quadratic Polynomial
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: mwirth@silver.cs.umanitoba.ca
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------

From mwirth@cs.umanitoba.ca Mon Oct 10 14:05:17 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["8900" "" " 4" "October" "1994" "20:52:02" "GMT" "Michael Wirth" "mwirth@cs.umanitoba.ca" "<36sf9i$suj@canopus.cc.umanitoba.ca>" "311" "Image Registration Method - Cubic Polynomial" "^From:" nil nil "10" "1994100420:52:02" "Image Registration Method - Cubic Polynomial" (number " " mark "     Michael Wirth     Oct  4  311/8900  " thread-indent "\"Image Registration Method - Cubic Polynomial\"\n") nil]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!concert!news.duke.edu!MathWorks.Com!europa.eng.gtefsd.com!library.ucla.edu!news.mic.ucla.edu!unixg.ubc.ca!quartz.ucs.ualberta.ca!tribune.usask.ca!canopus.cc.umanitoba.ca!cs.umanitoba.ca!mwirth
Newsgroups: sci.image.processing
Organization: Computer Science, University of Manitoba, Winnipeg, Canada
Lines: 311
Distribution: world
Message-ID: <36sf9i$suj@canopus.cc.umanitoba.ca>
NNTP-Posting-Host: silver.cs.umanitoba.ca
From: mwirth@cs.umanitoba.ca (Michael Wirth)
Subject: Image Registration Method - Cubic Polynomial
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: mwirth@silver.cs.umanitoba.ca
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------

From mwirth@cs.umanitoba.ca Mon Oct 10 14:05:24 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["5603" "" " 4" "October" "1994" "20:52:43" "GMT" "Michael Wirth" "mwirth@cs.umanitoba.ca" "<36sfar$suj@canopus.cc.umanitoba.ca>" "270" "Image Registration Technique - Surface Spline Transformation" "^From:" nil nil "10" "1994100420:52:43" "Image Registration Technique - Surface Spline Transformation" (number " " mark "     Michael Wirth     Oct  4  270/5603  " thread-indent "\"Image Registration Technique - Surface Spline Transformation\"\n") nil]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!concert!news.duke.edu!zombie.ncsc.mil!MathWorks.Com!europa.eng.gtefsd.com!library.ucla.edu!news.mic.ucla.edu!unixg.ubc.ca!quartz.ucs.ualberta.ca!tribune.usask.ca!canopus.cc.umanitoba.ca!cs.umanitoba.ca!mwirth
Newsgroups: sci.image.processing
Organization: Computer Science, University of Manitoba, Winnipeg, Canada
Lines: 270
Distribution: world
Message-ID: <36sfar$suj@canopus.cc.umanitoba.ca>
NNTP-Posting-Host: silver.cs.umanitoba.ca
From: mwirth@cs.umanitoba.ca (Michael Wirth)
Subject: Image Registration Technique - Surface Spline Transformation
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: mwirth@silver.cs.umanitoba.ca
Department of Computer Science          Snail: Box 654 Tache Residence
University of Manitoba			       University of Manitoba 
Winnipeg, Manitoba                             Winnipeg, MB  R3T 2N1 
CANADA  R3T 2N1                                
--------------------------------------------------------------------------------

From gary@maple.mc.com Fri Oct 28 00:53:43 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["650" "" "21" "October" "1994" "21:11:59" "GMT" "Gary Von Colln" "gary@maple.mc.com" "<389aqv$8ve@jericho.mc.com>" "23" "Re: TIFF library at Uni. of California, Berkeley" "^From:" nil nil "10" "1994102121:11:59" "TIFF library at Uni. of California, Berkeley" (number " " mark "     Gary Von Colln    Oct 21   23/650   " thread-indent "\"Re: TIFF library at Uni. of California, Berkeley\"\n") "<3838kc$go8@sun4.bham.ac.uk>"]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!darwin.sura.net!spool.mu.edu!uwm.edu!news.alpha.net!news.mathworks.com!noc.near.net!jericho.mc.com!maple!gary
Newsgroups: sci.image.processing
Organization: Mercury Computer Systems, Chelmsford MA 01824
Lines: 23
Distribution: world
Message-ID: <389aqv$8ve@jericho.mc.com>
References: <3838kc$go8@sun4.bham.ac.uk>
NNTP-Posting-Host: maple.mc.com
From: gary@maple.mc.com (Gary Von Colln)
Sender: gary@maple (Gary Von Colln)
Subject: Re: TIFF library at Uni. of California, Berkeley
Date: 21 Oct 1994 21:11:59 GMT

In article <3838kc$go8@sun4.bham.ac.uk>, S.S.Koh@BHAM.AC.UK (S.S.Koh) writes:
|> Dear netter, 
|>                   Could anybody direct me to the ftp location of the University of 
|> California, Berkeley. I am looking for Sam Leffler who wrote the TIFF library.
|> I would like to contact him concerning the TIFF library. Thank you in advance.
|> 
|> 
|> Best regards,
|> Edmund Koh
|> (s.s.koh@bham.ac.uk)


I have a copy of the TIFF library code.  The README file gives the
address sam@sgi.com

The ftp sites for the code are

    sgi.com			graphics/tiff/v3.0.tar.Z
	(192.48.153.1)
    ucbvax.berkeley.edu		pub/tiff/v3.0.tar.Z
	(128.32.130.12)



From metherall@aol.com Fri Oct 28 00:55:49 1994
Status: RO
X-VM-v5-Data: ([nil nil nil nil nil nil nil nil nil]
	["320" "" "25" "October" "1994" "08:34:06" "-0400" "metherall@aol.com" "metherall@aol.com" "<38itvu$evt@newsbf01.news.aol.com>" "8" "Re: BMP File Format Searched" "^From:" nil nil "10" "1994102512:34:06" "BMP File Format Searched" (number " " mark "     metherall@aol.com Oct 25    8/320   " thread-indent "\"Re: BMP File Format Searched\"\n") "<1994Oct17.135735.1999@cm.cf.ac.uk>"]
	nil)
Path: saips.cv.nrao.edu!hearst.acc.Virginia.EDU!caen!math.ohio-state.edu!howland.reston.ans.net!news.sprintlink.net!EU.net!uunet!newstf01.cr1.aol.com!newsbf01.news.aol.com!not-for-mail
Newsgroups: sci.image.processing
Organization: America Online, Inc. (1-800-827-6364)
Lines: 8
Message-ID: <38itvu$evt@newsbf01.news.aol.com>
References: <1994Oct17.135735.1999@cm.cf.ac.uk>
NNTP-Posting-Host: newsbf01.news.aol.com
From: metherall@aol.com (METHERALL)
Sender: news@newsbf01.news.aol.com
Subject: Re: BMP File Format Searched
Date: 25 Oct 1994 08:34:06 -0400

In article <1994Oct17.135735.1999@cm.cf.ac.uk>, spemk1@thor.cf.ac.uk
(Mathias Koehrer) writes:

There is a book called "Encyclopedia of Graphic File Formats" by J.murray
and W. VanRyper and published by O'Reilly & Assoc.,Inc.(Sebastopol,CA),
that has information on the BMP file format (along with many others).

Jeanne