//***************************************************************************
#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"
//***************************************************************************
typedef struct mydiff_ {
	int value;
	int id;
	} mydiff;
static int findknearest(int list[], int N, int k);
//***************************************************************************
//static int compareIntegers (int *e1, int *e2)
static int compareIntegers (const void *e1, const void *e2)
{
return (*((int *) e1) - *((int *) e2));
}
//***************************************************************************
//static int compareIntegers2 (mydiff *e1, mydiff *e2)
static int compareIntegers2 (const void *e1, const void *e2)
{
return (((mydiff *) e1)->value - ((mydiff *) e2)->value);
}
//***************************************************************************
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;
}
//***************************************************************************
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;
}
//***************************************************************************
void GrayImage::kNearestNeighboringFilterNxN (int N)	// odd N
{
kNearestNeighboringFilterNxN(N, N*N/2+1);
}
//***************************************************************************
void GrayImage::kNearestNeighboringFilterNxN (int N, int k2)	// odd N
{

cout << "Applying " << N << "x" << N << " k-Nearest (k = " << k2 
		<< ") Neigboring filter...";
cout.flush();

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];
		 dummy[i][j] = findknearest(list, m, k2);
		}

for (int i = 0; i < width; ++i)
	for (int j = 0; j < height; ++j)
		p[i][j] = (unsigned char) dummy[i][j];

cout << "done!" << endl;
}
//***************************************************************************
static int findknearest(int list[], int m, int k)
{
int center = list[m/2];	// odd N only
mydiff diff[m];
for (int i = 0; i < m; ++i)
	{
	 diff[i].id = i;
	 diff[i].value = abs(list[i]-center);
	}

qsort(diff,m,sizeof(mydiff),compareIntegers2);

int value = 0;
if (k > m) 
	k = m;
for (int i = 0; i < k; ++i)
	{
	 value += list[diff[i].id];
	}
value /= k;

return value;
}
//***************************************************************************
// A NxN Gaussian filter from
// Machine Vision text, pages 132-137
//***************************************************************************
void GrayImage::GaussianFilterNxN (int N)               // odd N
{
cout << "Applying " << N << "x" << N << " Gaussian filter...";
cout.flush();

double temp[N][N];
int mask[N][N];
double out[width][height];
int i, j, k, l;
int N2 = N/2;
int normfactor = 0;
double sd = 1.0;

for (i = -N2; i <= N2; i++)
	for (j = -N2; j <= N2; j++)
          {
            temp[i+N2][j+N2] = exp(-1.0*((i*i+j*j)/(2.0*(sd*sd))));
            mask[i+N2][j+N2] = (int) rint(1.0/temp[0][0]*temp[i+N2][j+N2]);
	    normfactor += mask[i+N2][j+N2];
          }

int printmask = 1;
if (printmask)
	{
	 printf("\n");
	 for (i = -N2; i <= N2; i++)
		{
	 	 for (j = -N2; j <= N2; j++)
			printf("%g ",mask[i+N2][j+N2]*1.0/normfactor);
 	 	 printf("\n");
		}
	 printf("normfactor = %d\n",normfactor);
	}

for (i = 0; i < width; i++)
	for (j = 0; j < height; j++)
		out[i][j] = 0.0;

for (i = N2; i < width - N2; i++)
	for (j = N2; j < height - N2; j++)
		for (k = -N2; k <= N2; k++)
			for (l = -N2; l <= N2; l++)
				out[i][j] += p[i+k][j+l]*mask[k+N2][l+N2]
					*1.0/normfactor;

for (i = N2; i < width - N2; i++)
	for (j = N2; j < height - N2; j++)
		p[i][j] = (unsigned char) out[i][j];

cout << "done!" << endl;
}
//***************************************************************************