Subject: Image Registration Header Files

From: [email-masked] (Michael Wirth)
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 ([email-masked])
 * 
 * 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: [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                                
--------------------------------------------------------------------------------