image.cc.html
//***************************************************************************
#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 "image.h"
#include <stdio.h>
//***************************************************************************
#define BUF_SIZE 256
//***************************************************************************
static int compareIntegers (int *e1, int *e2);
//***************************************************************************
GrayImage::GrayImage ( ) // Constructor
{
height = width = numpixels = 0;
for (int i = 0; i < 10; i++)
MaxHist[i] = MinHist[i] = 0;
numMaxHist = numMinHist = 0;
mincount = maxcount = 0;
numComponents = 0;
allocateFlag = 0;
}
//***************************************************************************
GrayImage::GrayImage( int w, int h, int numLev )
{
width = w;
height = h;
numLevels = numLev;
numBits = (int) (log((double) (numLevels+2))/log(2.0));
numpixels = width*height;
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 (numpixels)
{
for (int i = 0; i < width; i++)
delete p[i];
delete p;
}
if (hist)
delete hist;
if (numComponents)
{
for (int i = 0; i < numComponents; i++)
delete Phi[i];
delete Phi;
}
allocateFlag = 0;
}
//***************************************************************************
void GrayImage::allocateMem ( )
{
if (allocateFlag)
{
cerr << "allocateMem: Memory already allocated!!" << endl;
return;
}
allocateFlag = 1;
p = new unsigned char *[width];
for (int i = 0; i < width; i++)
{
p[i] = new unsigned char [height];
if (!p[i])
{
cerr << "GrayImage: Not enough memory for allocateMem!!"
<< endl;
exit(-1);
}
for (int j = 0; j < height; 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 < width; i++)
for (int j = 0; j < height; j++)
p[i][j] = 0;
for (int i = 0; i < numLevels+1; i++)
hist[i] = 0;
}
//***************************************************************************
void GrayImage::Init (int w, int h, int numLev )
{
if (allocateFlag)
deAllocateMem( );
width = w;
height = h;
numpixels = width*height;
numLevels = numLev;
allocateMem( );
for (int i = 0; i < width; i++)
for (int j = 0; j < height; j++)
p[i][j] = 0;
for (int i = 0; i < numLevels+1; i++)
hist[i] = 0;
}
//***************************************************************************
int GrayImage::ReadFromFile (char *FileName)
{
int i,j;
char buf[BUF_SIZE];
ifstream inp(FileName);
if (!inp) // Invalid FileName
{
cerr << "Invalid filname: \"" << FileName << "\" does not exist!!"
<< endl;
exit(-1);
}
cout << "Reading image from file \"" << FileName << "\" " << endl;
inp.getline(magicNumber,BUF_SIZE);
inp.getline(buf,BUF_SIZE);
while (buf[0] == '#')
inp.getline(buf,BUF_SIZE);
width = atoi(buf);
height = atoi(strpbrk(buf, " \t"));
numpixels = width*height;
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 < width; i++)
inp.read(p[i],height);
}
else // ASCII
{
cout << "GrayScale ASCII PGM format]";
for (i = 0; i < width; i++)
for (j = 0; j < height; 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 < width; i++)
for (j = 0; j < height; j++)
{
inp.read(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 < width; i++)
for (j = 0; j < height; 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 < width; i++)
for (j = 0; j < height; j += 8)
{
unsigned char pix[1];
inp.read(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 < width; i++)
for (j = 0; j < height; j++)
{
int pix;
inp >> pix;
if (pix == 0)
p[i][j] = 0;
else
p[i][j] = 255;
}
}
}
inp.close();
cout << endl;
cout << "Width = " << width << endl;
cout << "Height = " << height << endl;
cout << "Numpixels = " << numpixels << endl;
numLevels = 255;
numBits = 8;
magicNumber[1] = '5';
return 0; // No error encountered....successful completion
}
//***************************************************************************
int GrayImage::WriteToFile (char *FileName)
{
char buf[BUF_SIZE];
ofstream out(FileName);
if (!out)
return -1; // Invalid Filename
cout << "Writing " << width << "x" << height << " image to file \"" <<
FileName << "\"...";
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 < width; i++)
out.write(p[i],height);
out.close();
cout << "done!" << endl;
return 0; // No error encountered....successful completion
}
//***************************************************************************
int GrayImage::Threshold (unsigned char T)
{
cout << "Applying threshold for gray level " << (int) T << "....";
for (int i = 0; i < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] < T)
p[i][j] = 0;
cout << "done!" << endl;
}
//***************************************************************************
int 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 < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] < T1 || p[i][j] > T2)
p[i][j] = 0;
cout << "done!" << endl;
}
//***************************************************************************
int 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 < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] < T1 || p[i][j] > T2)
p[i][j] = 0;
else
p[i][j] = 255;
cout << "done!" << endl;
}
//***************************************************************************
int GrayImage::ThresholdBinary (unsigned char T)
{
for (int i = 0; i < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] < T)
p[i][j] = 0;
else
p[i][j] = 255;
}
//***************************************************************************
int GrayImage::reverseThreshold (unsigned char T)
{
for (int i = 0; i < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] >= T)
p[i][j] = 0;
}
//***************************************************************************
int 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 < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] >= T1 && p[i][j] <= T2)
p[i][j] = 0;
cout << "done!" << endl;
}
//***************************************************************************
int GrayImage::Histogram ( )
{
cout << "Calculating the histogram.....";
for (int i = 0; i <= numLevels; i++)
hist[i] = 0;
for (int i = 0; i < width; i++)
for (int j = 0; j < height; j++)
++hist[p[i][j]];
cout << "done!" << endl;
}
//***************************************************************************
int 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
*/
//***************************************************************************
GrayImage::smoothenHistogram ( )
{
int kn, nn, nkn;
double mu, blend, muk, munk;
double p[numLevels+1];
cout << "Applying Bezier Smoothing to the Histogram...";
cout.flush();
for (int i = 0; i <= numLevels; i++)
{
p[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--;
}
}
p[i] += hist[k] * blend;
}
hist[i] = (int) p[i];
}
cout << "done!" << endl;
cout.flush();
}
//***************************************************************************
int GrayImage::calcHistogramPeaks ( )
{
for (int i = 0; i < 10; i++)
MaxHist[i] = MinHist[i] = 0;
numMaxHist = numMinHist = 0;
mincount = maxcount = 0;
//if (hist[0] < numpixels/10000)
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;
}
//***************************************************************************
int 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;
width = I.width;
height = I.height;
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 < width; i++)
for (int j = 0; j < height; 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.height == height) && (I1.width == width))
{
*I = I1;
for (int i = 0; i < width; i++)
for (int j = 0; j < height; 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.height == height) && (I1.width == width))
{
*I = I1;
for (int i = 0; i < width; i++)
for (int j = 0; j < height; 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 width 1
//***************************************************************************
{
int perimeter = 0;
for (int i = 0; i < width; i++) // checks horizontal borders for
{
if (p[i][0] == F) // foreground pixels
perimeter++;
if (p[i][height-1] == F)
perimeter++;
}
for (int j = 1; j < height-1; j++) // checks vertical borders for
{
if (p[0][j] == F) // foreground pixels
perimeter++; // doesn't double-count corners
if (p[width-1][j] == F)
perimeter++;
}
for (int i = 1; i < width-1; i++) // checks rest of pixels
for (int j = 1; j < height-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 width 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 < height; i++)
for (int j = 0; j < width; 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 < width; k++)
for (int l = 0; l < height; 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 != width-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 != height-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 < width; i++)
for (int j = 0; j < height; j++)
if (p[i][j] == P)
++area;
return area;
}
//***************************************************************************
void GrayImage::ConvertToBinary ( )
{
cout << "Converting image to Binary...";
for (int i = 0; i < width; i++)
for (int j = 0; j < height; 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 < width; i++)
for (int j = 0; j < height; 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 *[width];
for (int i = 0; i < width; i++)
p2[i] = new unsigned int [height];
int interval = 16777215 / numComponents;
// Calculate gray scale gap between objects
for (int i = 0; i < width; i++)
for (int j = 0; j < height; 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 " << width << "x" << height
<< " image to file \"" << FileName << "\"...";
out << "P6" << endl; // For 8-bit GrayScale
out << "# Created by GrayImage::WriteToFile" << endl;
out << width << " " << height << endl;
out << numLevels << endl;
for (int i = 0; i < width; i++)
{
for (int j = 0; j < height; j++)
{
unsigned char r,g,b;
int P = p2[i][j];
r = P % 256;
g = (P/256) % 256;
b = P / (256*256);
out.write(&r,1);
out.write(&g,1);
out.write(&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 < width; i++)
for (int j = 0; j < height; 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 < width; i++)
for (int j = 0; j < height; 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 = width, i_max = -1;
int j_min = height, j_max = -1;
for (int i = 0; i < width; i++)
for (int j = 0; j < height; 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;
}
/*
for (int j = j_min; j <= j_max; j++)
{
p[i_min][j] = P;
p[i_max][j] = P;
}
for (int i = i_min; i <= i_max; i++)
{
p[i][j_min] = P;
p[i][j_max] = P;
}
*/
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 = width, i_max = -1;
int j_min = height, j_max = -1;
for (int i = 0; i < width; i++)
for (int j = 0; j < height; 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::MeanFilterNxN (int N) // odd N
{
cout << "Applying " << N << "x" << N << " mean filter...";
int dummy[width][height];
int n2 = N/2;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
{
dummy[i][j] = 0;
int m = 0;
for (int k = -n2; k <= n2; ++k)
for (int l = -n2; l <= n2; ++l)
if (!(i+k < 0 || i+k >= width || j+l < 0
|| j+l >= height))
{
dummy[i][j] += p[i+k][j+l];
++m;
}
dummy[i][j] /= m;
}
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::MedianFilterNxN (int N) // odd N
{
cout << "Applying " << N << "x" << N << " median filter...";
int dummy[width][height];
int n2 = N/2;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
{
int list[N*N];
int m = 0;
for (int k = -n2; k <= n2; ++k)
for (int l = -n2; l <= n2; ++l)
if (!(i+k < 0 || i+k >= width || j+l < 0
|| j+l >= height))
list[m++] = p[i+k][j+l];
qsort(list,m,sizeof(int),compareIntegers);
if (m % 2)
dummy[i][j] = (list[m/2]+list[m/2-1])/2;
else
dummy[i][j] = list[m/2];
}
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
}
//***************************************************************************
static int compareIntegers (int *e1, int *e2)
{
return (*e1 - *e2);
}
//***************************************************************************
void GrayImage::DilationNxN_Binary (int N, int num) // odd N
{ // num -> object with label "num"
cout << " Applying " << N << "x" << N << " Dilation operator to component "
<< num << "...";
int dummy[width][height];
int n2 = N/2;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
{
dummy[i][j] = p[i][j];
for (int k = -n2; k <= n2; ++k)
for (int l = -n2; l <= n2; ++l)
{
if (!(i+k < 0 || i+k >= width || j+l < 0
|| j+l >= height)
&& p[i+k][j+l] == num)
{
dummy[i][j] = num;
break;
}
}
}
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::ErosionNxN_Binary (int N, int num) // odd N
{ // num -> object with label "num"
cout << " Applying " << N << "x" << N << " Erosion operator to component "
<< num << "...";
int dummy[width][height];
int n2 = N/2;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
{
dummy[i][j] = p[i][j];
for (int k = -n2; k <= n2; ++k)
for (int l = -n2; l <= n2; ++l)
{
if (!(i+k < 0 || i+k >= width || j+l < 0
|| j+l >= height)
&& p[i+k][j+l] != num)
{
dummy[i][j] = 0;
break;
}
}
}
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::OpeningNxN_Binary (int N, int num) // odd N, num => component number
{
cout << "Applying " << N << "x" << N << " OPENING operator to component "
<< num << "..." << endl;
ErosionNxN_Binary ( N, num );
DilationNxN_Binary ( N, num );
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::ClosingNxN_Binary (int N, int num) // odd N, num => component number
{
cout << "Applying " << N << "x" << N << " CLOSING operator to component "
<< num << "..." << endl;
DilationNxN_Binary ( N, num );
ErosionNxN_Binary ( N, num );
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::OpeningNxN_Binary (int N) // odd N
{
GrayImage final(width, height, numLevels);
GrayImage original;
original = *this;
for (int i = 1; i <= numComponents; i++)
{
GrayImage dummy;
dummy = original;
dummy.OpeningNxN_Binary ( N, i );
final = final + dummy;
}
*this = final;
}
//***************************************************************************
void GrayImage::ClosingNxN_Binary (int N) // odd N
{
GrayImage final(width, height, numLevels);
GrayImage original;
original = *this;
for (int i = 1; i <= numComponents; i++)
{
GrayImage dummy;
dummy = original;
dummy.ClosingNxN_Binary ( N, i );
final = final + dummy;
}
*this = final;
}
//***************************************************************************
void GrayImage::ErosionNxN_Binary (int N) // odd N
{
GrayImage final(width, height, numLevels);
GrayImage original;
original = *this;
for (int i = 1; i <= numComponents; i++)
{
GrayImage dummy;
dummy = original;
dummy.ErosionNxN_Binary ( N, i );
final = final + dummy;
}
*this = final;
}
//***************************************************************************
void GrayImage::DilationNxN_Binary (int N) // odd N
{
for (int i = 1; i <= numComponents; i++)
DilationNxN_Binary ( N, i );
}
//***************************************************************************
void GrayImage::ErosionNxN (int N) // odd N
{
cout << " Applying " << N << "x" << N << " Local Minimum (Erosion) filter...";
int dummy[width][height];
int n2 = N/2;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
{
dummy[i][j] = INT_MAX;
for (int k = -n2; k <= n2; ++k)
for (int l = -n2; l <= n2; ++l)
if (!(i+k < 0 || i+k >= width || j+l < 0
|| j+l >= height)
&& p[i+k][j+l] < dummy[i][j])
dummy[i][j] = p[i+k][j+l];
}
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::DilationNxN (int N) // odd N
{
cout << " Applying " << N << "x" << N << " Local Maximum (Dilation) filter...";
int dummy[width][height];
int n2 = N/2;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
{
dummy[i][j] = INT_MIN;
for (int k = -n2; k <= n2; ++k)
for (int l = -n2; l <= n2; ++l)
if (!(i+k < 0 || i+k >= width || j+l < 0
|| j+l >= height)
&& p[i+k][j+l] > dummy[i][j])
dummy[i][j] = p[i+k][j+l];
}
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::OpeningNxN (int N) // odd N
{
cout << "Applying " << N << "x" << N << " OPENING operator..." << endl;
ErosionNxN( N );
DilationNxN( N );
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::ClosingNxN (int N) // odd N
{
cout << "Applying " << N << "x" << N << " CLOSING operator..." << endl;
DilationNxN( N );
ErosionNxN( N );
cout << "done!" << endl;
}
//***************************************************************************
void GrayImage::Skeletonization ( ) // Medial Axis Thinning
{
cout << "Applying Skeletonization to binary image...";
int **dummy, **dummyk, **dummyk1, **ptrk, **ptrk1;
dummy = new int *[width];
dummyk = new int *[width];
dummyk1 = new int *[width];
for (int i = 0; i < width; i++)
{
dummy[i] = new int [height];
dummyk[i] = new int [height];
dummyk1[i] = new int [height];
}
ptrk = dummyk;
ptrk1 = dummyk1;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
if (p[i][j] == 0)
dummy[i][j] = ptrk[i][j] = ptrk1[i][j] = 0;
else
dummy[i][j] = ptrk[i][j] = ptrk1[i][j] = 1;
int k = 1;
int flag;
do
{
flag = 1;
ptrk = dummyk;
ptrk1 = dummyk1;
for (int i = 1; i < width-1; ++i)
for (int j = 1; j < height-1; ++j)
if (ptrk1[i][j] == k)
{
int min = ptrk1[i-1][j];
if (min > ptrk1[i+1][j])
min = ptrk1[i+1][j];
if (min > ptrk1[i][j+1])
min = ptrk1[i][j+1];
if (min > ptrk1[i][j-1])
min = ptrk1[i][j-1];
if (min == k)
{
ptrk[i][j] = min + 1;
flag = 0;
}
}
k++; /* Next Iteration */
ptrk = dummyk;
ptrk1 = dummyk1;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
if (ptrk[i][j] == k)
ptrk1[i][j] = k;
} while (flag == 0);
ptrk = dummyk;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = ptrk[i][j];
for (int i = 1; i < width-1; ++i)
for (int j = 1; j < height-1; ++j)
{
unsigned char val = p[i][j];
if (val)
{
unsigned char max = p[i-1][j];
if (max < p[i+1][j])
max = p[i+1][j];
if (max < p[i][j+1])
max = p[i][j+1];
if (max < p[i][j-1])
max = p[i][j-1];
if (val >= max)
dummy[i][j] = 255;
else
dummy[i][j] = 0;
}
}
for (int i = 0; i < width; ++i)
dummy[i][0] = dummy[i][height-1] = 0;
for (int j = 0; j < height; ++j)
dummy[0][j] = dummy[width-1][j] = 0;
for (int i = 0; i < width; ++i)
for (int j = 0; j < height; ++j)
p[i][j] = dummy[i][j];
cout << "done!" << endl;
for (int i = 0; i < width; i++)
{
delete dummy[i];
delete dummyk[i];
delete dummyk1[i];
}
delete dummy;
delete dummyk;
delete dummyk1;
}
//***************************************************************************