//***************************************************************************
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>             /* ANSI C standard library routines */
#include <string.h>             /* ANSI standard string routines */
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include "image.h"
#include "matrix.h"
//***************************************************************************
Matrix::Matrix ( )      // Constructor
{
allocateFlag = 0;
Optimization2D = 0;
}
//***************************************************************************
Matrix::Matrix (int m, int n)   // Constructor
{
Init (m, n);
}
//***************************************************************************
Matrix::~Matrix ( )     // Destructor
{
deAllocateMem( );
}
//***************************************************************************
void Matrix::Init (int m, int n)        // Called by Constructor
{
allocateFlag = 0;
Optimization2D = 0;
nx = m;
ny = n;
allocateMem (nx, ny);
}
//***************************************************************************
void Matrix::allocateMem (int m, int n)
{
if (allocateFlag)
        {
         cerr << "Matrix::allocateMem: Memory already allocated!!" << endl;
         return;
        }
allocateFlag = 1;

p = new float *[m];
for (int i = 0; i < m; i++)
        {
         p[i] = new float [n];
         if (!p[i])
                {
                 cerr << "Matrix: Not enough memory for allocateMem!!"
                        << endl;
                 exit(-1);
                }
         for (int j = 0; j < n; j++)
                p[i][j] = 0.0;
        }
}
//***************************************************************************
void Matrix::deAllocateMem ( )
{
if (!allocateFlag)
        return;

for (int i = 0; i < nx; i++)
        delete p[i];
delete p;
allocateFlag = 0;
}
//***************************************************************************
void Matrix::Identity ()
{
if (nx != ny)
        {
         cerr << "Identity(): Nx not equal to Ny !!" << endl;
         exit(-1);
        }
for (int i = 0; i < nx; i++)
        {
         for (int j = 0; j < ny; j++)
                p[i][j] = 0.0;
         p[i][i] = 1.0;
        }
}
//***************************************************************************
void Matrix::translate (float x, float y, float z)
{
*this = *this * Translate(x,y,z);
}
//***************************************************************************
void Matrix::scale (float x, float y, float z)
{
*this = *this * Scale(x,y,z);
}
//***************************************************************************
void Matrix::rotate (float angle, float x, float y, float z)
{
*this = *this * Rotate(angle,x,y,z);
}
//***************************************************************************
Matrix &Matrix::Transpose ()
{
static Matrix *M = new Matrix(ny,nx);

for (int i = 0; i < ny; i++)
        for (int j = 0; j < nx; j++)
                M->p[i][j] = p[j][i];
return *M;
}
//***************************************************************************
int Matrix::operator== (Matrix &M)
{
if (M.nx != nx || M.ny != ny)
        {
         cerr << "operator ==: Both matrices are not of the same size!!" << endl;
         exit(-1);
        }

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                if (p[i][j] != M.p[i][j])
                        return 0;

return 1;
}
//***************************************************************************
void Matrix::operator+= (Matrix &M)
{
if (M.nx != nx || M.ny != ny)
        {
         cerr << "operator +=: Both matrices are not of the same size!!" << endl;
         exit(-1);
        }

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] += M.p[i][j];

}
//***************************************************************************
Matrix &Matrix::operator= (const Matrix &M)
{
if (this == &M)
        return *this;

nx = M.nx;
ny = M.ny;

if (!allocateFlag)
        allocateMem (nx, ny);

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] = M.p[i][j];

return *this;
}
//***************************************************************************
Matrix &Matrix::operator+ (Matrix &M1)
{
static Matrix *M = new Matrix;

if (M1.nx != nx || M1.ny != ny)
        {
         cerr << "operator +: Both matrices are not of the same size!!" << endl;
         exit(-1);
        }
*M = M1;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                M->p[i][j] = p[i][j] + M1.p[i][j];
return *M;
}
//***************************************************************************
Matrix &Matrix::operator- (Matrix &M1)
{
static Matrix *M = new Matrix;

if (M1.nx != nx || M1.ny != ny)
        {
         cerr << "operator -: Both matrices are not of the same size!!" << endl;
         exit(-1);
        }
*M = M1;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                M->p[i][j] = p[i][j] - M1.p[i][j];
return *M;
}
//***************************************************************************
// Pixel by pixel multiplication
Matrix &Matrix::operator^ (Matrix &M1)
{
static Matrix *M = new Matrix;

if (M1.nx != nx || M1.ny != ny)
        {
         cerr << "operator ^: Both matrices are not of the same size!!" << endl;
         exit(-1);
        }
*M = M1;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                M->p[i][j] = p[i][j] * M1.p[i][j];
return *M;
}
//***************************************************************************
// Pixel by pixel division
Matrix &Matrix::operator/ (Matrix &M1)
{
static Matrix *M = new Matrix;

if (M1.nx != nx || M1.ny != ny)
        {
         cerr << "operator /: Both matrices are not of the same size!!" << endl;
         exit(-1);
        }
*M = M1;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                M->p[i][j] = p[i][j] / M1.p[i][j];
return *M;
}
//***************************************************************************
Matrix &Matrix::operator* (float k)
{
static Matrix *M = new Matrix(nx,ny);

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                M->p[i][j] = p[i][j]*k;
return *M;
}
//***************************************************************************
Matrix &Matrix::operator/ (float k)
{
static Matrix *M = new Matrix(nx,ny);

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                M->p[i][j] = p[i][j]/k;
return *M;
}
//***************************************************************************
Matrix &Matrix::operator* (Matrix &M1)
{
if (ny != M1.nx)
        {
         cerr << "operator *: Cannot multiply " << nx << "x" << ny
                << " matrix by " << M1.nx << "x" << M1.ny << " matrix!!"
                << endl;
         exit(-1);
        }
static int mx = nx, my = M1.ny;
static Matrix *M = new Matrix(mx,my);
if (mx != nx || my != M1.ny)
        {
         mx = nx; my = M1.ny;
         delete M;
         M = new Matrix(mx,my);
        }
for (int i = 0; i < M->nx; i++)
        for (int j = 0; j < M->ny; j++)
                {
                 M->p[i][j] = 0.0;
                 for (int k = 0; k < M1.nx; k++)
                        M->p[i][j] += p[i][k]*M1.p[k][j];
                }
return *M;
}
//***************************************************************************
ostream &operator<< (ostream &Out, Matrix &M)
{
for (int i = 0; i < M.nx; i++)
        {
         for (int j = 0; j < M.ny; j++)
                {
                 Out.setf(ios::fixed,ios::floatfield);
                 Out.precision(4);
                 if (M.p[i][j] >= 0.0)
                        Out << " ";
                 Out << M.p[i][j] << " ";
                }
         Out << endl;
        }
Out << endl;
Out.precision(6);
Out.setf(0,ios::floatfield);

return Out;
}
//***************************************************************************
float Matrix::operator() (int i, int j)
{
return p[i][j];
}
//***************************************************************************
float *Matrix::arrayform()
{
float *a = new float [nx*ny];
float *b;

for (int j = 0; j < ny; j++)
        for (int i = 0; i < nx; i++)
                {
                 a[j*nx+i] = p[i][j];
                }
b = a;
return b;
}
//***************************************************************************
Matrix &Transpose (Matrix &M)
{
return M.Transpose();
}
//***************************************************************************
Matrix &Identity (int n)
{
static Matrix *M = new Matrix(n,n);

M->Identity();

return *M;
}
//***************************************************************************
Matrix &Rotate (float angle, float x, float y, float z)
{
float r = sqrt(x*x+y*y+z*z);
if (fabs(r-1.0) > 1.0e-6)
        {
         x = x/r;
         y = x/r;
         z = x/r;
        }
float c = cos(angle);
float s = sin(angle);

static Matrix *M = new Matrix(4,4);

M->p[0][0] = x*x*(1.0-c)+c;
M->p[0][1] = x*y*(1.0-c)-z*s;
M->p[0][2] = x*z*(1.0-c)+y*s;
M->p[0][3] = 0.0;
M->p[1][0] = y*x*(1.0-c)+z*s;
M->p[1][1] = y*y*(1.0-c)+c;
M->p[1][2] = y*z*(1.0-c)-x*s;
M->p[1][3] = 0.0;
M->p[2][0] = z*x*(1.0-c)-y*s;
M->p[2][1] = z*y*(1.0-c)+x*s;
M->p[2][2] = z*z*(1.0-c)+c;
M->p[2][3] = 0.0;
M->p[3][0] = 0.0;
M->p[3][1] = 0.0;
M->p[3][2] = 0.0;
M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &xRotate (float angle)
{
static Matrix *M = new Matrix(4,4);

M->p[1][1] = M->p[2][2] = cos(angle);
M->p[1][2] = sin(angle);
M->p[2][1] = -sin(angle);
M->p[0][0] = M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &yRotate (float angle)
{
static Matrix *M = new Matrix(4,4);

M->p[0][0] = M->p[2][2] = cos(angle);
M->p[0][2] = sin(angle);
M->p[2][0] = -sin(angle);
M->p[1][1] = M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &zRotate (float angle)
{
static Matrix *M = new Matrix(4,4);

M->p[0][0] = M->p[1][1] = cos(angle);
M->p[0][1] = sin(angle);
M->p[1][0] = -sin(angle);
M->p[2][2] = M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &Translate (float x, float y, float z)
{
static Matrix *M = new Matrix(4,4);

M->p[0][3] = x;
M->p[1][3] = y;
M->p[2][3] = z;
M->p[0][0] = M->p[1][1] = M->p[2][2] = M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &Scale (float x, float y, float z)
{
static Matrix *M = new Matrix(4,4);

M->p[0][0] = x;
M->p[1][1] = y;
M->p[2][2] = z;
M->p[0][0] = M->p[1][1] = M->p[2][2] = M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &Perspective (float lambda)
{
static Matrix *M = new Matrix(4,4);

M->p[3][2] = -1.0/lambda;
M->p[0][0] = M->p[1][1] = M->p[2][2] = M->p[3][3] = 1.0;

return *M;
}
//***************************************************************************
Matrix &TransformToImagePlane (float x0, float y0, float z0, float xC,
        float yC, float zC, float r1, float r2, float r3, float lambda)
{
static Matrix *M = new Matrix;

float x = xC - x0;
float y = yC - y0;
float z = zC - z0;
float r = sqrt(x*x+y*y+z*z);
float theta = atan2(x,-y);
float alpha = acos(z/r);

if (x == 0.0 && y == 0.0)       // bugfix for Linux egcs-c++ atan2(,) function
        theta = 0.0;

Matrix P = Perspective(lambda);
Matrix Tw0 = Translate(-x0,-y0,-z0);
Matrix R_theta = zRotate(theta);
Matrix R_alpha = xRotate(alpha);
Matrix Tr = Translate(-r1,-r2,-r3);

Matrix T = P * Tr * R_alpha * R_theta * Tw0;

//cout << "alpha = " << alpha/DEG2RAD << " deg" << endl;
//cout << "theta = " << theta/DEG2RAD << " deg" << endl;
 
*M = T;
return *M;
}
//***************************************************************************
void PerspectiveTransform (Matrix T, float x, float y, float z, 
        float *xp, float *yp)
{
Matrix w(4,1), c_h(4,1);

w.p[0][0] = x;
w.p[1][0] = y;
w.p[2][0] = z;
w.p[3][0] = 1.0;

c_h = T * w;

if (fabs(c_h(3,0)) > 1.0e-4)
        {
         *xp = c_h(0,0)/c_h(3,0); 
         *yp = c_h(1,0)/c_h(3,0); 
        }
else
        {
         *xp = 1.0e+10;
         *yp = 1.0e+10;
        }
}
//***************************************************************************
void Matrix::translate_and_rotate2D (float x, float y, float angle)
{
p[0][0] = p[1][1] = cos(angle);
p[1][0] = sin(angle);
p[0][1] = -p[1][0];
p[0][2] = x;
p[1][2] = y;
}
//***************************************************************************
void Matrix::transformXY (int x, int y, int *xf, int *yf)
{
if (Optimization2D)
        {
         *xf = (int) (x*p[0][0] + y*p[0][1] + p[0][3]);
         *yf = (int) (x*p[1][0] + y*p[1][1] + p[1][3]);
        }
else
        {
         Matrix w(4,1), v(4,1);

         w.p[0][0] = x;
         w.p[1][0] = y;
         w.p[2][0] = 0.0;
         w.p[3][0] = 1.0;

         v = *this * w;

         *xf = (int) (v(0,0)/v(3,0));
         *yf = (int) (v(1,0)/v(3,0));
        }
}
//***************************************************************************
// Or and Oi are changed
void MatrixFFT2D(Matrix *Or, Matrix *Oi, Matrix *Ir, Matrix *Ii)
{
MatrixFFT2D(Or, Oi, Ir, Ii, 1);
}
//***************************************************************************
// Or and Oi are changed
void MatrixIFFT2D(Matrix *Or, Matrix *Oi, Matrix *Ir, Matrix *Ii)
{
MatrixFFT2D(Or, Oi, Ir, Ii, -1);
}
//***************************************************************************
static void four1(float *data, int nn, int isign);
//***************************************************************************
// 2D Fourier Transform of data with real part stored in Ir 
// and imaginary part in Ii.
// isign =  1 => Forward FFT; isign = -1 => Inverse FFT
// Or and Oi are changed
void MatrixFFT2D(Matrix *Or, Matrix *Oi, Matrix *Ir, Matrix *Ii, int isign)
{
int nx = Ir->nx;
int ny = Ir->ny;

Matrix T(2*nx,ny);
float *tmp1 = new float [2*ny+1];
float *tmp2 = new float [2*nx+1];

for (int i = 1; i <= nx; i++)
        {
         for (int j = 1; j <= ny; j++)
                {
                 tmp1[j*2-1] = Ir->p[i-1][j-1];
                 tmp1[j*2] = Ii->p[i-1][j-1];
                }
         four1(tmp1, ny, isign);
         for (int j = 1; j <= ny; j++)
                {
                 T.p[i*2-2][j-1] = tmp1[j*2-1];
                 T.p[i*2-1][j-1] = tmp1[j*2];
                }
        }

for (int i = 1; i <= ny; i++)
        {
         for (int j = 1; j <= nx; j++)
                {
                 tmp2[j*2-1] = T.p[j*2-2][i-1];
                 tmp2[j*2] = T.p[j*2-1][i-1];
                }
         four1(tmp2,nx,isign);
         for (int j = 1; j <= nx; j++)
                {
                 Or->p[j-1][i-1] = tmp2[j*2-1];
                 Oi->p[j-1][i-1] = tmp2[j*2];
                }
        }

delete tmp1;
delete tmp2;
}
//***************************************************************************
// From NRC
static void four1(float *data, int nn, int isign)
{
float wtemp, wr, wpr, wpi, wi, theta;
float tempr;
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

int n = nn << 1;
int j = 1;

for (int i = 1; i < n; i += 2)
        {
         if (j > i)
                {
                 SWAP(data[j],data[i]);
                 SWAP(data[j+1],data[i+1]);
                }
         int m = n >> 1;
         while (m >= 2 && j > m) 
                {
                 j -= m;
                 m >>= 1;
                }
         j += m;
        }

int mmax = 2;
while (n > mmax)
        {
         int istep = 2*mmax;
         theta = 2*Pi/(isign*mmax);
         wtemp = sin(0.5*theta);
         wpr = -2.0*wtemp*wtemp;
         wpi = sin(theta);
         wr = 1.0;
         wi = 0.0;
         for (int m = 1; m < mmax; m += 2)
                {
                 for (int i = m; i <= n; i += istep)
                        {
                         j = i+mmax;
                         tempr = wr*data[j]-wi*data[j+1];
                         float tempi = wr*data[j+1]+wi*data[j];
                         data[j] = data[i]-tempr;
                         data[j+1] = data[i+1]-tempi;
                         data[i] += tempr;
                         data[i+1] += tempi;
                        }
                 wr = (wtemp=wr)*wpr-wi*wpi+wr;
                 wi = wi*wpr+wtemp*wpi+wi;
                }
         mmax = istep;
        }
#undef SWAP
}
//***************************************************************************
// Pixel by Pixel product A = B * C
// A is changed, B and C are not
void MatProduct(Matrix *A, Matrix *B, Matrix *C)
{
for (int h = 0; h < A->nx; h++)
         for (int w = 0; w < A->ny; w++)
                A->p[h][w] = B->p[h][w] * C->p[h][w];
}
//***************************************************************************
// Pixel by Pixel product A = B + C
// A is changed, B and C are not
void MatSum(Matrix *A, Matrix *B, Matrix *C)
{
for (int h = 0; h < A->nx; h++)
         for (int w = 0; w < A->ny; w++)
                A->p[h][w] = B->p[h][w] + C->p[h][w];
}
//***************************************************************************
// Pixel by Pixel product A = B - C
// A is changed, B and C are not
void MatDiff(Matrix *A, Matrix *B, Matrix *C)
{
for (int h = 0; h < A->nx; h++)
         for (int w = 0; w < A->ny; w++)
                A->p[h][w] = B->p[h][w] - C->p[h][w];
}
//***************************************************************************
// v is changed, m is not
void Matrix2Vector(Matrix *m, float *v)
{
int nh = m->nx;
int nw = m->ny;
for (int h = 0; h < nh; h++) 
        for (int w = 0; w < nw; w++) 
                v[h*nw+w] = (float) m->p[h][w];
}
//***************************************************************************
void MatCopy(Matrix *A, Matrix *B, int h_target, int w_target, 
                int h_begin, int w_begin, int h_end, int w_end)
{
int Anx = A->nx, Any = A->ny;
int Bnx = B->nx, Bny = B->ny;

if ((h_target >= 0) && (h_target < Anx) && (w_target >= 0) && (w_target < Any))
        {
         if ((h_begin >= 0) && (h_begin < Bnx) && (w_begin >= 0) && (w_begin < Bny))
                {
                 int h = h_end-h_begin+1;
                 int w = w_end-w_begin+1;
                 if ((h >= 1) && (w >= 1))
                        {
                         int h_done = h_target+h-1;
                         int w_done = w_target+w-1;
                         if ((h_done < Anx) && (w_done < Any))
                                {
                                 for (int i = 0; i < h; i++)
                                        {
                                         for (int j = 0; j < w; j++)
                                                {
                                                 A->p[i+h_target][j+w_target] = B->p[i+h_begin][j+w_begin];
                                                }
                                        }
                                }
                        }
                } 
        }
else
        {
         cerr << "MatCopy: Matrix dimension error (Source matrix is " 
                << A->nx << "x" << A->ny << ", Dest matrix is "
                << B->nx << "x" << B->ny << ")!" << endl;
         cerr << "\tParameters are (" << h_target << "," << w_target
                << "," << h_begin << "," << w_begin << "," << h_end
                << "," << w_end << ")" << endl;

         exit(-1);
        }
}           
//***************************************************************************
float log2(double x)
{
return log10(x)/log10(2.0);
}
//***************************************************************************
void Matrix::WriteImage(char *filename, char *fmt)
{
GrayImage img;
img.ReadFromMatrix(*this);
img.WriteToFile(filename, fmt);
}
//***************************************************************************
void Matrix::WriteImage(char *filename)
{
GrayImage img;
img.ReadFromMatrix(*this);
img.WriteToFile(filename);
}
//***************************************************************************
void Matrix::WriteFFTImage(char *filename, char *fmt)
{
GrayImage img;
Matrix tmp(nx, ny);

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                tmp.p[i][j] = log10(1.0+fabs(p[i][j]));

img.ReadFromMatrix(tmp);
img.WriteToFile(filename, fmt);
}
//***************************************************************************
void Matrix::WriteFFTImage(char *filename)
{
GrayImage img;
Matrix tmp(nx, ny);

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                tmp.p[i][j] = log10(1.0+fabs(p[i][j]));

img.ReadFromMatrix(tmp);
img.WriteToFile(filename);
}
//***************************************************************************
void Matrix::Abs()
{
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] = fabs(p[i][j]);
}
//***************************************************************************
// Real is replaced by Magnitude, Imag is replaced by Phase
void Complex2Polar(Matrix *Real, Matrix *Imag)
{
int nx = Real->nx;
int ny = Real->ny;

if (nx != Imag->nx || ny != Imag->ny)
        {
         cerr << "Magnitude: Dimensions of the input and output matrices"
                << " do not match!!" << endl;
        }

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                {
                 double magnitude = sqrt(Real->p[i][j]*Real->p[i][j]
                        + Imag->p[i][j]*Imag->p[i][j]);
                 double phase = atan2(Imag->p[i][j], Real->p[i][j]);
                 Real->p[i][j] = magnitude;
                 Imag->p[i][j] = phase;
                }
/*
// Do phase unwrapping to remove discontinuities at +Pi/-Pi.
#define TOL 1.0e-3
for (int i = 1; i < nx; i++)
        for (int j = 0; j < ny; j++)
                {
                 if (fabs((Imag->p[i][j]-Imag->p[i-1][j]) - TwoPi) < TOL)
                        Imag->p[i-1][j] += TwoPi;
                }
*/
}
//***************************************************************************
float Matrix::Norm()
{
double norm = 0.0;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                norm += p[i][j]*p[i][j];
norm = sqrt(norm);

return (float) norm;
}
//***************************************************************************
// Normalize between 0 to 1
void Matrix::Normalize()
{
double minp = 1.0e+300, maxp = -1.0e+300;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                {
                 if (p[i][j] > maxp)
                        maxp = p[i][j]; 
                 if (p[i][j] < minp) 
                        minp = p[i][j];
                }
double range = maxp-minp;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] = (p[i][j]-minp)/range;
}
//***************************************************************************
// Normalize based on enery (norm2)
void Matrix::NormalizeEnergy()
{
float energy = 0.0;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                energy += fabs(p[i][j]);
//cout << "Energy factor = " << energy << endl;

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] /= energy;
}
//***************************************************************************
// Normalize based on enery (norm)
void Matrix::Normalize2()
{
double norm = 0.0;
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                norm += fabs(p[i][j]);
norm = sqrt(norm);

//cout << "Normalization factor = " << norm << endl;

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] /= (float) norm;
}
//***************************************************************************
void Matrix::scale(float k)
{
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] *= k;
}
//***************************************************************************
void Matrix::Zero()
{
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                p[i][j] = 0.0;
}
//***************************************************************************
void Matrix::Threshold(float k)
{
for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                {
                 if (p[i][j] > k)
                        p[i][j] = 1.0;
                 else
                        p[i][j] = 0.0;
                }
}
//***************************************************************************
Matrix *allocateMatrixArray(int n, int height, int width)
{
Matrix *M;
M = new Matrix [n];
for (int i = 0; i < n; i++)
        {
         M[i].Init(height, width);
        }
return M;
}
//**************************************************************************
void deallocateMatrixArray(Matrix *M, int n)
{
for (int i = 0; i < n; i++)
        {
         M[i].deAllocateMem();
        }
delete M;
}
//**************************************************************************
void printMatrixArray(Matrix *M, char *filename, int n, int height, int width)
{
Matrix Big(n*height, width);

for (int i = 0; i < n; i++)
        MatCopy(&Big, &M[i], i*height, 0, 0, 0, height-1, width-1);

Big.WriteImage(filename, "gif");
}
//**************************************************************************
Matrix **allocateMatrixArray(int nx, int ny, int height, int width)
{
Matrix **M;
M = new Matrix *[nx];
for (int i = 0; i < nx; i++)
        {
         M[i] = new Matrix [ny];
         for (int j = 0; j < ny; j++)
                M[i][j].Init(height, width);
        }
return M;
}
//**************************************************************************
void deallocateMatrixArray(Matrix **M, int nx, int ny)
{
for (int i = 0; i < nx; i++)
        {
         for (int j = 0; j < ny; j++)
                {
                 M[i][j].deAllocateMem();
                }
         delete M[i];
        }
delete M;
}
//**************************************************************************
void printMatrixArray(Matrix **M, char *filename, int nx, int ny, int height, int width)
{
Matrix Big(nx*height, ny*width);

for (int i = 0; i < nx; i++)
        for (int j = 0; j < ny; j++)
                MatCopy(&Big, &M[i][j], i*height, j*width, 0, 0, height-1, width-1);

Big.WriteImage(filename, "gif");
}
//**************************************************************************
int Matrix::WriteToFile(char *FileName)
{
ofstream out(FileName);

if (!out)
        return -1;              // Invalid Filename

cout << "Writing " << nx << "x" << ny << " matrix to file \"" <<
        FileName << "\"...";

out << nx << " " << ny << endl;
for (int i = 0; i < nx; i++)
        out.write((char *) p[i], sizeof(float)*ny);

out.close();

cout << "done!" << endl;

return 0;               // No error encountered....successful completion
}
//***************************************************************************
int Matrix::ReadFromFile(char *FileName)
{
const int buf_size = 40;
char buf[buf_size];
ifstream inp(FileName);

if (!inp)       // Invalid FileName
        {
         cerr << "Matrix::ReadFromFile: Invalid filname: \"" << FileName << "\" does not exist!!"
                << endl;
         exit(-1);
        }

inp.getline(buf,buf_size);
while (buf[0] == '#')   // ignore comments
        inp.getline(buf,buf_size);

nx = atoi(buf);
ny = atoi(strpbrk(buf, " \t"));

cout << "Reading " << nx << "x" << ny << " matrix from file \"" 
        << FileName << "\"...";

if (!allocateFlag)
        allocateMem(nx, ny);

for (int i = 0; i < nx; i++)
        inp.read((char *) p[i], sizeof(float)*ny);

cout << "done!" << endl;

return 0;       // No error encountered....successful completion
}
//***************************************************************************
void Write2DTecplotFile(Matrix *R, Matrix *I, char *filename)
{
ofstream out(filename);

cout << "Writing two " << R->nx << "x" << R->ny << " matrices to a 2D TECPLOT file \"" <<
        filename << "\"...";

int height = R->nx;
int width = R->ny;

out << "Title = \"Anirudh's Image Library\"" << endl;
out << "Variables = \"X\", \"Y\", \"Ex\", \"Ey\"" << endl;
out << "zone i = " << height << " j = " << width << endl;

for (int j = 0; j < width; j++)
        {
         for (int i = 0; i < height; i++)
                {
                 out << i << " " << j << " " << R->p[i][j] << " " << I->p[i][j] << endl;
                }
         out << endl;
        }

out.close();

cout << "done!" << endl;
}
//***************************************************************************
int Matrix::isBoundaryPixel(int x, int y)
{
if (p[x][y] == 0.0)
        return 0;

int count = 0;
for (int i = -1; i <= 1; i++)
        for (int j = -1; j <= 1; j++)
                {
                 int xn = x + i;
                 int yn = y + j;
                 if (xn >= 0 && xn < nx && yn >=0 && yn < ny &&
                        !(xn == x && yn == y) && p[xn][yn])
                        ++count;
                }

if (count == 1)
        return 1;
else
        return 0;
}
//***************************************************************************