//***************************************************************************
#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 "utils.h"
//***************************************************************************
#define BUF_SIZE 256
//***************************************************************************
GrayImage::GrayImage ( )        // Constructor
{
width = height = numpixels = 0;
for (int i = 0; i < 10; i++)
        MaxHist[i] = MinHist[i] = 0;
numMaxHist = numMinHist = 0;
hist = NULL;
mincount = maxcount = 0;
numComponents = 0;
allocateFlag = 0;
}
//***************************************************************************
GrayImage::GrayImage( int h, int w )
{
GrayImage(h, w, 255);
}
//***************************************************************************
GrayImage::GrayImage( int h, int w, int numLev )
{
height = h;
width = w;
numLevels = numLev;
numBits = (int) (log((double) (numLevels+2))/log(2.0));
numpixels = height*width;
for (int i = 0; i < 10; i++)
        MaxHist[i] = MinHist[i] = 0;
numMaxHist = numMinHist = 0;
mincount = maxcount = 0;
numComponents = 0;
allocateFlag = 0;
allocateMem( );
}
//***************************************************************************
GrayImage::~GrayImage ( )       // Destructor
{
deAllocateMem( );
}
//***************************************************************************
void GrayImage::deAllocateMem ( )
{
if (!allocateFlag)
        return;

if (numpixels)
        {
         for (int i = 0; i < height; i++)
                delete p[i];
         delete p;
        }
if (hist)
        {
         delete hist;
         hist = NULL;
        }
if (numComponents)
        {
         for (int i = 0; i < numComponents; i++)
                delete Phi[i];
         delete Phi;
         numComponents = 0;
        }
allocateFlag = 0;
}
//***************************************************************************
void GrayImage::allocateMem ( )
{
if (allocateFlag)
        {
         cerr << "allocateMem: Memory already allocated!!" << endl;
         return;
        }
allocateFlag = 1;

p = new unsigned char *[height];
for (int i = 0; i < height; i++)
        {
         p[i] = new unsigned char [width];
         if (!p[i])
                {
                 cerr << "GrayImage: Not enough memory for allocateMem!!" 
                        << endl;
                 exit(-1);
                }
         for (int j = 0; j < width; j++)
                p[i][j] = 0;
        }
hist = new int [numLevels+1];
for (int i = 0; i < numLevels+1; i++)
        hist[i] = 0;
}
//***************************************************************************
void GrayImage::Init ( )
{
if (!allocateFlag)
        allocateMem( );

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                p[i][j] = 0;
for (int i = 0; i < numLevels+1; i++)
        hist[i] = 0;
}
//***************************************************************************
void GrayImage::Init (int h, int w, int numLev )
{
if (allocateFlag)
        deAllocateMem( );

height = h;
width = w;
numpixels = height*width;
numLevels = numLev;

allocateMem( );

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                p[i][j] = 0;
for (int i = 0; i < numLevels+1; i++)
        hist[i] = 0;
}
//***************************************************************************
int GrayImage::isRGB (char *FileName)
{
ifstream inp(FileName);

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

char a;
inp.read((char *) &a,1);

if (a != 'P')   // not a PBM/PGM/PPM image, hence assuming RGB
        return 1;

inp.getline(magicNumber,BUF_SIZE);
int flag = 1;

if (!(magicNumber[1] == '6' || magicNumber[1] == '3'))  // not RGB image
        flag = 0;

return flag;
}
//***************************************************************************
// Read from file in any valid image format recognized by ImageMagick
int GrayImage::ReadFromFile2 (char *FileName)
{
cout << "Reading image from file \"" << FileName << "\":" << endl;

char *tmpfile = convert_image_format(FileName, "ppm");
ReadFromFile(tmpfile);
deleteFile(tmpfile);

return 0;
}
//***************************************************************************
// Read from file in PBM/PGM/PPM format
int GrayImage::ReadFromFile (char *FileName)
{
int i,j;
char buf[BUF_SIZE];
ifstream inp(FileName);

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

if (!strstr(FileName, "mytmp"))
        cout << "Reading image from file \"" << FileName << "\":" << endl;

inp.getline(magicNumber,BUF_SIZE);

inp.getline(buf,BUF_SIZE);

int buflen = strlen(buf);
while (buf[0] == '#' || buflen < 2)
        {
         inp.getline(buf,BUF_SIZE);
         buflen = strlen(buf);
        }

width = atoi(buf);
height = atoi(strpbrk(buf, " \t"));

numpixels = height*width;

if (magicNumber[1] == '1' || magicNumber[1] == '4')
        numLevels = 1;
else
        {
         inp.getline(buf,BUF_SIZE);
         numLevels = atoi(buf);
        }

if (!allocateFlag)
        allocateMem( );

numBits = (int) (log((double) (numLevels+2))/log(2.0));
cout << "[" << numBits << "-bit ";

if (magicNumber[1] == '5' || magicNumber[1] == '2')     // GrayScale image
        {
         if (magicNumber[1] == '5')             // RAWBITs
                {
                 cout << "GrayScale RAWBITs PGM format]";
                 for (i = 0; i < height; i++)
                        inp.read((char *) p[i],width);
                }
         else                                   // ASCII
                {
                 cout << "GrayScale ASCII PGM format]";
                 for (i = 0; i < height; i++)
                        for (j = 0; j < width; j++)
                                {
                                 int pix;
                                 inp >> pix;
                                 p[i][j] = (unsigned char) 
                                                (pix / pow( 2, numBits-8 ));
                                }
                }
        }
else if (magicNumber[1] == '6' || magicNumber[1] == '3')        // RGB image
        {
         unsigned char rgb[3];
         if (magicNumber[1] == '6')             // RAWBITs
                {
                 cout << "RGB RAWBITs PPM format]";
                 for (i = 0; i < height; i++)
                        for (j = 0; j < width; j++)
                                {
                                 inp.read((char *) rgb,3);
                                 p[i][j] = (unsigned char) (sqrt((double) 
                                        (rgb[0]*rgb[0] + rgb[1]*rgb[1] 
                                        + rgb[2]*rgb[2])) /sqrt(3.0));
                                }
                }
         else                                   // ASCII
                {
                 cout << "RGB ASCII PPM format]";
                 for (i = 0; i < height; i++)
                        for (j = 0; j < width; j++)
                                {
                                 int pix;
                                 for (int k = 0; k < 3; k++)
                                        {
                                         inp >> pix;
                                         rgb[k] = (unsigned char) pix;
                                        }
                                 p[i][j] = (unsigned char) (sqrt((double) 
                                        (rgb[0]*rgb[0] + rgb[1]*rgb[1] 
                                        + rgb[2]*rgb[2])) /sqrt(3.0));
                                }
                }
        }
else if (magicNumber[1] == '4' || magicNumber[1] == '1')   // Binary image
        {
         if (magicNumber[1] == '4')             // RAWBITs
                {
                 cout << "Binary RAWBITs PBM format]";
                 for (i = 0; i < height; i++)
                        for (j = 0; j < width; j += 8)
                                {
                                 unsigned char pix[1];
                                 inp.read((char *) pix,1);
                                 unsigned int x = (int) pix[0];
                                 for (int k = 0; k < 8; k++)
                                        {
                                         int y = x/((int) pow(2,7-k));
                                         if (y)
                                                p[i][j+k] = 0;
                                         else
                                                p[i][j+k] = 255;
                                         x -= y*((int)pow(2,7-k));
                                        }
                                }
                }
         else
                {
                 cout << "Binary ASCII PBM format]";
                 for (i = 0; i < height; i++)
                        for (j = 0; j < width; j++)
                                {
                                 int pix;
                                 inp >> pix;
                                 if (pix == 0)
                                        p[i][j] = 0;
                                 else
                                        p[i][j] = 255;
                                }
                }
        }

inp.close();

cout << endl;
cout << "Width = " << height << endl;
cout << "Height = " << width << endl;
cout << "Numpixels = " << numpixels << endl;

numLevels = 255;
numBits = 8;
magicNumber[1] = '5';

return 0;               // No error encountered....successful completion
}
//***************************************************************************
int GrayImage::ReadFromMatrix(Matrix &M)
{
height = M.nx;
width = M.ny;

numpixels = height*width;
numLevels = 255;
numBits = 8;

double minp = 1.0e+300, maxp = -1.0e+300;
for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                {
                 if (M.p[i][j] > maxp)
                        maxp = M.p[i][j]; 
                 if (M.p[i][j] < minp)
                        minp = M.p[i][j]; 
                }
cout << endl;
cout << "Matrix2Image: [Min, Max] for " << height << "x" << width << " matrix = [" 
        << minp << "," << maxp << "]" << endl;
cout.flush();

if (!allocateFlag)
        allocateMem( );

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                p[i][j] = (unsigned char) (255.0*(M.p[i][j]-minp)/(maxp-minp));

magicNumber[1] = '5';

return 0;               // No error encountered....successful completion
}
//***************************************************************************
int GrayImage::WriteToMatrix (Matrix *M)
{
if ((M->nx < height) || (M->ny < width))
        {
         cerr << "Error: Matrix dimensions [" << M->nx << "x"
                << M->ny 
                << "] is not compatible with Image dimension ["
                << height << "x" << width << "] (HxW)!!" << endl;
         exit(-1);
        }

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                M->p[i][j] = (double) p[i][j];

return 0;               // No error encountered....successful completion
}
//***************************************************************************
int GrayImage::WriteToFile(char *FileName, char *fmt)
{
WriteToFile(FileName);
char *tmpfile = convert_image_format(FileName, fmt);
moveFile(tmpfile, FileName);

return 0;
}
//***************************************************************************
int GrayImage::WriteToFile (char *FileName)
{
char buf[BUF_SIZE];
ofstream out(FileName);

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

cout << "Writing " << height << "x" << width << " (HxW) image to file \"" << 
        FileName << "\"...";
cout.flush();

if (magicNumber[0] == 'P')
        out << magicNumber[0] << magicNumber[1] << endl; // For 8-bit GrayScale
else
        out << "P5" << endl;
out << "# Created by GrayImage::WriteToFile" << endl;
out << width << " " << height << endl;
out << numLevels << endl;

for (int i = 0; i < height; i++)
        out.write((char *) p[i],width);

out.close();

cout << "done!" << endl;

return 0;               // No error encountered....successful completion
}
//***************************************************************************
// all pixels above the specified threshold "T" become 255 (white)
void GrayImage::ThresholdWhite (unsigned char T)
{
cout << "Applying threshold for gray level " << (int) T << "....";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] >= T)
                        p[i][j] = 255;

cout << "done!" << endl;
}
//***************************************************************************
// all pixels below the specified threshold "T" become 0 (black)
void GrayImage::Threshold (unsigned char T)
{
cout << "Applying threshold for gray level " << (int) T << "....";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] < T)
                        p[i][j] = 0;

cout << "done!" << endl;
}
//***************************************************************************
// all pixels below the specified threshold "T1" and above "T2" become 0 (black)
void GrayImage::Threshold (unsigned char T1, unsigned char T2)
{
if (T2 == 0)
        T2 = 255;
if (T1 > T2)
        { 
         unsigned char temp = T1;
         T1 = T2;
         T2 = temp; 
        }

cout << "Applying threshold for gray levels between " << (int) T1 <<
        " and " << (int) T2 << "....";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] < T1 || p[i][j] > T2)
                        p[i][j] = 0;

cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::ThresholdBinary (unsigned char T1, unsigned char T2)
{
if (T2 == 0)
        T2 = 255;
if (T1 > T2)
        { 
         unsigned char temp = T1;
         T1 = T2;
         T2 = temp; 
        }

cout << "Applying binary threshold for gray levels between " << (int) T1 <<
        " and " << (int) T2 << "....";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] < T1 || p[i][j] > T2)
                        p[i][j] = 0;
                else
                        p[i][j] = 255;

cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::ThresholdBinary (unsigned char T)
{
for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] < T)
                        p[i][j] = 0;
                else
                        p[i][j] = 255;
}
//***************************************************************************
// all pixels above the specified threshold "T" become 0 (black)
void GrayImage::reverseThreshold (unsigned char T)
{
for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] >= T)
                        p[i][j] = 0;
}
//***************************************************************************
// all pixels above the specified threshold "T1" and below "T2" become 0 (black)
void GrayImage::reverseThreshold (unsigned char T1, unsigned char T2)
{
if (T2 == 0)
        T2 = 255;
if (T1 > T2)
        { 
         unsigned char temp = T1;
         T1 = T2;
         T2 = temp; 
        }

cout << "Applying threshold for all gray levels except those between " 
        << (int) T1 << " and " << (int) T2 << "....";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] >= T1 && p[i][j] <= T2)
                        p[i][j] = 0;

cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::Histogram ( )
{
cout << "Calculating the histogram.....";

for (int i = 0; i <= numLevels; i++)
        hist[i] = 0;
for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                ++hist[p[i][j]];

cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::printHistogram (char *FileName)
{
ofstream out(FileName);

cout << "Writing histogram data for the image to file \"" 
        << FileName << "\"...";

for (int i = 0; i <= numLevels; i++)
        out << i << "\t" << hist[i] << endl;

cout << "done!" << endl;
}
//***************************************************************************
/*
 *              BEZIER CURVES
 *              Ref : Mathematical Elements for Computer Graphics
 *                      2nd Ed. - Rogers and Adams
 *                      pp 295-298.
 *
 *                                                      -Anirudh Modi
 */
//***************************************************************************
void GrayImage::smoothenHistogram ( )
{
int kn, nn, nkn;
double mu, blend, muk, munk;
double *h = new double [numLevels+1];

cout << "Applying Bezier Smoothing to the Histogram...";
cout.flush();

for (int i = 0; i <= numLevels; i++)
        {
         h[i] = 0.0;
         mu = i*1.0/(numLevels+1);
         muk = 1.0;
         munk = pow( 1.0-mu, numLevels );
         for (int k = 0; k <= numLevels; k++) 
                {
                 nn = numLevels;
                 kn = k;
                 nkn = numLevels - k;
                 blend = muk * munk;
                 muk *= mu;
                 munk /= 1.0 - mu;
                 while (nn >= 1) 
                        {
                         blend *= nn;
                         nn--;
                         if (kn > 1) 
                                {
                                 blend /= (double) kn;
                                 kn--;
                                }
                         if (nkn > 1) 
                                {
                                 blend /= (double) nkn;
                                 nkn--;
                                }
                        }
                 h[i] += hist[k] * blend;
                }
         hist[i] = (int) h[i];
        }

delete h;

cout << "done!" << endl;
cout.flush();
}
//***************************************************************************
void GrayImage::calcHistogramPeaks ( )
{
for (int i = 0; i < 10; i++)
        MaxHist[i] = MinHist[i] = 0;

numMaxHist = numMinHist = 0;
mincount = maxcount = 0;

if (hist[0] < hist[1])
        MinHist[++numMinHist] = 0;

#define THR 10 
for (int i = 1; i < numLevels; i++)
        {
         if ((hist[i-1] <= hist[i]) && (hist[i+1] < hist[i]))
                {
                 MaxHist[++numMaxHist] = i;
                 if (numMinHist && 
                        (MaxHist[numMaxHist] - MinHist[numMinHist] < THR))
                        --numMinHist;
                 if (numMaxHist > 1 &&
                        (MaxHist[numMaxHist] - MaxHist[numMaxHist-1] < THR))
                         --numMaxHist;
                }
         else if ((hist[i-1] >= hist[i]) && (hist[i+1] > hist[i]))
                {
                 MinHist[++numMinHist] = i;
                 if (numMaxHist && 
                        (MinHist[numMinHist] - MaxHist[numMaxHist] < THR))
                        --numMaxHist;
                 if (numMinHist > 1 &&
                        (MinHist[numMinHist] - MinHist[numMinHist-1] < THR))
                         --numMinHist;
                }
        }
if (numLevels - MinHist[numMinHist] > 25)
        MinHist[++numMinHist] = numLevels;

int List[20], minmax[20], n = 0, l = 1, m = 1;

while (l <= numMinHist && m <= numMaxHist)
        if (MinHist[l] < MaxHist[m])
                {
                 List[n] = MinHist[l++];
                 minmax[n++] = -1;
                }
        else
                {
                 List[n] = MaxHist[m++];
                 minmax[n++] = 1;
                }
if (l > numMinHist)
        while (m <= numMaxHist)
                {
                 List[n] = MaxHist[m++];
                 minmax[n++] = 1;
                }
else if (m > numMaxHist)
        while (l <= numMinHist)
                {
                 List[n] = MinHist[l++];
                 minmax[n++] = -1;
                }

for (int i = 0; i < n-1; i++)
        {
         if (minmax[i] == minmax[i+1])
                {
                 if (minmax[i] < 0)
                        List[i+1] = (hist[List[i]] < hist[List[i+1]]) 
                                        ? List[i] : List[i+1];
                 else
                        List[i+1] = (hist[List[i]] > hist[List[i+1]]) 
                                        ? List[i] : List[i+1];
                 minmax[i] = 0;
                }
        }

numMinHist = numMaxHist = 0;
for (int i = 0; i < n; i++)
        if (minmax[i] < 0)
                MinHist[++numMinHist] = List[i];
        else if (minmax[i] > 0)
                MaxHist[++numMaxHist] = List[i];

cout << endl;
for (int i = 1; i <= numMaxHist; i++)
         cout << "Maxima " << i << " for histogram is at " 
                << MaxHist[i] << " with value of "
                        << hist[MaxHist[i]] << endl;
for (int i = 1; i <= numMinHist; i++)
         cout << "Minima " << i << " for histogram is at " 
                << MinHist[i] << " with value of "
                        << hist[MinHist[i]] << endl;
cout << endl;
}
//***************************************************************************
void GrayImage::printHistogramPeaks (char *FileName)
{
ofstream out(FileName);

cout << "Writing histogram peaks for the image to file \""
        << FileName << "\"...";

for (int i = 1; i <= numMinHist; i++)
        out << MinHist[i] << "\t" << hist[MinHist[i]] << endl;
out << endl;
for (int i = 1; i <= numMaxHist; i++)
        out << MaxHist[i] << "\t" << hist[MaxHist[i]] << endl;

cout << "done!" << endl;
}
//***************************************************************************
int GrayImage::getHistMinima ( )
{
int n;

if (mincount < numMinHist)
        n = MinHist[++mincount];
else
        n = -1;

return n;
}
//***************************************************************************
int GrayImage::getHistMaxima ( )
{
int n;

if (maxcount < numMaxHist)
        n = MaxHist[++maxcount];
else
        n = -1;

return n;
}
//***************************************************************************
int GrayImage::getHistMinima ( int i )
{
int n;

if (i <= numMinHist)
        n = MinHist[i];
else
        n = -1;

return n;
}
//***************************************************************************
int GrayImage::getHistMaxima ( int i )
{
int n;

if (i <= numMaxHist)
        n = MaxHist[i];
else
        n = -1;

return n;
}
//***************************************************************************
GrayImage &GrayImage::operator= (const GrayImage &I)
{
if (this == &I)
        return *this;

height = I.height;
width = I.width;
numpixels = I.numpixels;
numBits = I.numBits;
numLevels = I.numLevels;
mincount = I.mincount;
maxcount = I.maxcount;
numMaxHist = I.numMaxHist;
numMinHist = I.numMinHist;
allocatePhi(I.numComponents);
numComponents = I.numComponents;

for (int i = 0; i < 2; i++)
        magicNumber[i] = I.magicNumber[i];
for (int i = 0; i < 10; i++)
        {
         MaxHist[i] = I.MaxHist[i];
         MinHist[i] = I.MinHist[i];
        }

if (!allocateFlag)
        allocateMem( );

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                p[i][j] = I.p[i][j];
for (int i = 0; i <= numLevels; i++)
        hist[i] = I.hist[i];

for (int i = 0; i < numComponents; i++)
        for (int j = 0; j < 7; j++)
                Phi[i][j] = I.Phi[i][j];

return *this;
}
//***************************************************************************
GrayImage &GrayImage::operator+ (GrayImage &I1)
{
GrayImage *I = new GrayImage;

if ((I1.width == width) && (I1.height == height))
        {
         *I = I1;
         for (int i = 0; i < height; i++)
                 for (int j = 0; j < width; j++)
                        I->p[i][j] += p[i][j];
         for (int i = 0; i <= numLevels; i++)
                I->hist[i] += hist[i];
        }
else
        {
         cerr << "operator +: Both images are not of the same size!!" << endl;
         exit(-1);
        }
return *I;
}
//***************************************************************************
GrayImage &GrayImage::operator- (GrayImage &I1)
{
GrayImage *I = new GrayImage;

if ((I1.width == width) && (I1.height == height))
        {
         *I = I1;
         for (int i = 0; i < height; i++)
                 for (int j = 0; j < width; j++)
                        I->p[i][j] = p[i][j] - I->p[i][j];
         for (int i = 0; i <= numLevels; i++)
                {
                 I->hist[i] = hist[i] - I->hist[i];
                 if (I->hist[i] < 0)
                        I->hist[i] = 0;
                }
        }
else
        {
         cerr << "operator -: Both images are not of the same size!!" << endl;
         exit(-1);
        }
return *I;
}
//***************************************************************************
int GrayImage::Perimeter (int F)
//***************************************************************************
// Assumptions:
// Background is 255, foreground is F
// Image does not have any regions of height 1
//***************************************************************************
{
int perimeter = 0;

for (int i = 0; i < height; i++)         // checks horizontal borders for
        {
         if (p[i][0] == F)               // foreground pixels
                perimeter++;
         if (p[i][width-1] == F)
                perimeter++;
        }

for (int j = 1; j < width-1; j++)      // checks vertical borders for
        {
         if (p[0][j] == F)               // foreground pixels
                perimeter++;            // doesn't double-count corners
         if (p[height-1][j] == F)
                perimeter++;
        }

for (int i = 1; i < height-1; i++)               // checks rest of pixels
        for (int j = 1; j < width-1; j++)      // for borders
                {
                 if ((p[i][j] == F)
                        && ((p[i-1][j] == 255) || (p[i+1][j] == 255)
                                || (p[i][j-1] == 255) || (p[i][j+1] == 255)))
                        perimeter++;
                }

return perimeter;
}
//***************************************************************************
int GrayImage::LabelComponents ( )
//***************************************************************************
// Assumptions:
// Background is 255, foreground is 0
// Only one region in image
// Image does not have any regions of height 1
//***************************************************************************
{
#define AREA_THRESHOLD 200      // Min area reqd to classify as a component

int num = 0;

cout << "Labeling components....";

// Find connected components
for (int i = 0; i < width; i++)
        for (int j = 0; j < height; j++)
                if (p[i][j] == 0)               // Is this a new component?
                        {
                         num++;
                         p[i][j] = num;
                                // Assign new label to pixel
                         FindNeighbors (i,j);
                                // Recursively find other connected pixels
                         if (Area(num) < AREA_THRESHOLD)
                                {
                                 for (int k = 0; k < height; k++)
                                        for (int l = 0; l < width; l++)
                                                if (p[k][l] == num)
                                                        p[k][l] = 255;
                                 --num;
                                }
                        }

cout << num << " component(s) found....done!" << endl;

for (int k = 1; k <= num; k++)
        cout << "   Component " << k << ": Area = " << Area(k) << 
                ", Perimeter = " << Perimeter(k) << endl; 

allocatePhi(num);
numComponents = num;

return numComponents;
}
//***************************************************************************
int GrayImage::FindNeighbors (int i, int j)
{
// Check top neighbor
if ((i != 0) && (p[i-1][j] == 0))   // Is top neighbor connected and unmarked?
        {
         p[i-1][j] = p[i][j];           // If so, mark it
         FindNeighbors (i-1,j);         // and find its unmarked neighbors
        }

// Check bottom neighbor
if ((i != height-1) && (p[i+1][j] == 0))        // Is bottom neighbor 
        {                                       // connected and unmarked?
         p[i+1][j] = p[i][j];           // If so, mark it
         FindNeighbors (i+1,j);         // and find its unmarked neighbors
        }

// Check left neighbor
if ((j != 0) && (p[i][j-1] == 0))  // Is left neighbor connected and unmarked?
        {
         p[i][j-1] = p[i][j];           // If so, mark it
         FindNeighbors (i,j-1);         // and find its unmarked neighbors
        }

// Check right neighbor
if ((j != width-1) && (p[i][j+1] == 0)) // Is right neighbor connected 
        {                                       // and unmarked?
         p[i][j+1] = p[i][j];           // If so, mark it
         FindNeighbors (i,j+1);         // and find its unmarked neighbors
        }

return 0;
}
//***************************************************************************
int GrayImage::Area (int P)
{
int area = 0;

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] == P)
                        ++area;

return area;
}
//***************************************************************************
void GrayImage::ConvertToBinary ( )
{
cout << "Converting image to Binary...";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j])
                        p[i][j] = 0;            // Black foreground
                else
                        p[i][j] = 255;          // White background

cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::ReverseBinary ( )
{
cout << "Reversing Binary image...";

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] == 0)
                        p[i][j] = 255;          // Black foreground
                else
                        p[i][j] = 0;            // White background

cout << "done!" << endl;
}
//***************************************************************************
int GrayImage::WriteLabelsToFile (char *FileName)
{
char buf[BUF_SIZE];
ofstream out(FileName);

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

unsigned int **p2;

p2 = new unsigned int *[height];
for (int i = 0; i < height; i++)
         p2[i] = new unsigned int [width];

int interval = 16777215 / numComponents;        
                        // Calculate gray scale gap between objects
for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] != 255)             // If pixel is labeled,
                        p2[i][j] = 16777215-p[i][j]*interval;   
                                        // then multiply by interval
                else
                        p2[i][j] = 16777215;

cout << "Writing " << height << "x" << width
        << " image to file \"" << FileName << "\"...";

out << "P6" << endl;    // For 8-bit GrayScale
out << "# Created by GrayImage::WriteToFile" << endl;
out << height << " " << width << endl;
out << numLevels << endl;

for (int i = 0; i < height; i++)
        {
         for (int j = 0; j < width; j++)
                {
                 unsigned char r,g,b;
                 int P = p2[i][j];
                 r = P % 256;
                 g = (P/256) % 256;
                 b = P / (256*256);
                 out.write((char *) &r,1);
                 out.write((char *) &g,1);
                 out.write((char *) &b,1);
                }
         delete p2[i];
        }
delete p2;

out.close();

cout << "done!" << endl;

return 0;               // No error encountered....successful completion
}
//***************************************************************************
double GrayImage::Moment (int x, int y, int P)
{
double moment = 0.0;

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] == P)
                        moment += pow(i,x)*pow(j,y);

return moment;
}
//***************************************************************************
double GrayImage::CentralizedMoment (int x, int y, int P)
{
double moment = 0.0;

double m00 = Moment(0,0,P);
double m10 = Moment(1,0,P);
double m01 = Moment(0,1,P);
int i_c = (int) (m10/m00);
int j_c = (int) (m01/m00);

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] == P)
                        moment += pow(i-i_c,x)*pow(j-j_c,y);

return moment;
}
//***************************************************************************
double GrayImage::NormalizedCentralMoment (int x, int y, int P)
{

if ((x+y) < 2)
        {
         cerr << "Normalized Central Moments cannot be computed for p+q < 2!"
                << endl;
         exit(-1);
        }

double gamma = (x+y)/2.0+1.0;

return CentralizedMoment(x,y,P)/pow(CentralizedMoment(0,0,P),gamma);
}
//***************************************************************************
void GrayImage::InvariantMoments (int P)
{
double n20 = CentralizedMoment(2,0,P);
double n02 = CentralizedMoment(0,2,P);
double n11 = CentralizedMoment(1,1,P);
double n12 = CentralizedMoment(1,2,P);
double n21 = CentralizedMoment(2,1,P);
double n30 = CentralizedMoment(3,0,P);
double n03 = CentralizedMoment(0,3,P);

double phi[8];

phi[1] = n20 + n02;
phi[2] = pow(n20-n02,2) + 4*pow(n11,2);
phi[3] = pow(n30-3*n12,2) + pow(3*n21-n03,2);
phi[4] = pow(n30+n12,2) + pow(n21+n03,2);
phi[5] = (n30-3*n12) *  (n30+n12) * (pow(n30+n12,2) - 3*pow(n21+n03,2))
          + (3*n12-n03) * (n21+n03) * (3*pow(n30+n12,2) - pow(n21+n03,2));
phi[6] = (n20-n02) * (pow(n30+n12,2) - pow(n21+n03,2)) 
          + 4*n11*(n30+n12)*(n21+n03);
phi[7] = (3*n21-n30) * (n30+n12) * (pow(n30+n12,2) - 3*pow(n21+n03,2))
          + (3*n12-n30)*(n21+n03)*(3*pow(n30+n12,2) - pow(n21+n03,2));

for (int i = 1; i <= 7; i++)
        phi[i] = log10(fabs(phi[i])+1.0);

cout << "------------------------------------------------" << endl;
cout << "Moment invariants for component " << P << " are:" << endl;
cout << "------------------------------------------------" << endl;
for (int i = 1; i <= 7; i++)
        cout << "   phi[" << i << "] = " << phi[i] << endl;
cout << "------------------------------------------------" << endl;

for (int i = 0; i < 7; i++)
        Phi[P-1][i] = phi[i+1];
}
//***************************************************************************
double GrayImage::getInvariantMoment (int P, int i)
{
return Phi[P-1][i-1];
}
//***************************************************************************
void GrayImage::MinBoundRectangle (int P)
{
int i_min = height, i_max = -1;
int j_min = width, j_max = -1;

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] == P)
                        {
                         if (i < i_min) i_min = i;
                         if (i > i_max) i_max = i;
                         if (j < j_min) j_min = j;
                         if (j > j_max) j_max = j;
                        }

cout << "   MBR for component " << P 
        << " is " << i_max-i_min << "x" << j_max-j_min << " box with area " 
        << (i_max-i_min)*(j_max-j_min) << " -> [" 
        << i_min << "," << j_min << "] and ["
        << i_max << "," << j_max << "]" << endl;

}
//***************************************************************************
double GrayImage::Elongation (int P)
{
int i_min = height, i_max = -1;
int j_min = width, j_max = -1;

for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                if (p[i][j] == P)
                        {
                         if (i < i_min) i_min = i;
                         if (i > i_max) i_max = i;
                         if (j < j_min) j_min = j;
                         if (j > j_max) j_max = j;
                        }

double ratio = (i_max - i_min)*1.0/(j_max-j_min);
if (ratio < 1.0)
        ratio = 1.0/ratio;

return ratio;
}
//***************************************************************************
void GrayImage::allocatePhi (int num)
{
if (num == numComponents)
        return;

if (numComponents)
        {
         for (int i = 0; i < numComponents; i++)
                delete Phi[i];
         delete Phi;
        }
Phi = new double *[num];
        for (int i = 0; i < num; i++)
                Phi[i] = new double [7];
}
//***************************************************************************
void GrayImage::setValue (unsigned char T)
{
for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++j)
                p[i][j] = T;
}
//***************************************************************************
void WriteMatrixToImageFile(Matrix &M, char *filename)
{
GrayImage img;
img.ReadFromMatrix(M);
img.WriteToFile(filename);
}
//***************************************************************************
void GrayImage::RemoveLonelyPixels(void)
{
for (int i = 1; i < height-1; ++i)
        for (int j = 1; j < width-1; ++j)
                {
                 // if all the eight neighbors are unconnected
                 if (p[i-1][j-1] == 0 && p[i-1][j] == 0 && p[i-1][j+1] == 0
                    && p[i+1][j-1] == 0 && p[i+1][j] == 0 && p[i+1][j+1] == 0
                    && p[i][j-1] == 0 && p[i][j+1] == 0)
                        p[i][j] = 0;
                }
}
//***************************************************************************
void GrayImage::Zero(void)
{
for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
                p[i][j] = 0;
                
}
//***************************************************************************
int GrayImage::isBoundaryPixel(int x, int y)
{
if (p[x][y] == 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 < height && yn >=0 && yn < width &&
                        !(xn == x && yn == y) && p[xn][yn])
                        ++count;
                }

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