//***************************************************************************
#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"
//***************************************************************************
int **allocateIntegerMatrix(int m, int n);
double **allocateDoubleMatrix(int m, int n);
void deallocateMatrix(int **p, int m, int n);
void deallocateMatrix(double **p, int m, int n);
//***************************************************************************
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 = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++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 >= height || j+l < 0 
                                        || j+l >= width)
                                        && p[i+k][j+l] == num)
                                        {
                                         dummy[i][j] = num;
                                         break;
                                        }
                                }
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
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 = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++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 >= height || j+l < 0 
                                        || j+l >= width)
                                        && p[i+k][j+l] != num)
                                        {
                                         dummy[i][j] = 0;
                                         break;
                                        }
                                }
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
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(height, width, 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(height, width, 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(height, width, 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 = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++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 >= height || j+l < 0 
                                        || j+l >= width)
                                        && p[i+k][j+l] < dummy[i][j])
                                        dummy[i][j] = p[i+k][j+l];
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
void GrayImage::DilationNxN (int N)             // odd N
{

cout << "   Applying " << N << "x" << N << " Local Maximum (Dilation) filter...";

int **dummy = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++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 >= height || j+l < 0 
                                        || j+l >= width)
                                        && p[i+k][j+l] > dummy[i][j])
                                        dummy[i][j] = p[i+k][j+l];
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
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 *[height];
dummyk = new int *[height];
dummyk1 = new int *[height];
for (int i = 0; i < height; i++)
        {
         dummy[i] = new int [width];
         dummyk[i] = new int [width];
         dummyk1[i] = new int [width];
        }

ptrk = dummyk;
ptrk1 = dummyk1;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++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 < height-1; ++i)
                for (int j = 1; j < width-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 < height; ++i)
                for (int j = 0; j < width; ++j)
                        if (ptrk[i][j] == k)
                                ptrk1[i][j] = k;
        } while (flag == 0);

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

for (int i = 1; i < height-1; ++i)
        for (int j = 1; j < width-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 < height; ++i)
        dummy[i][0] = dummy[i][width-1] = 0;
for (int j = 0; j < width; ++j)
        dummy[0][j] = dummy[height-1][j] = 0;

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

cout << "done!" << endl;

for (int i = 0; i < height; 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 = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++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 >= height || j+l < 0 
                                        || j+l >= width))
                                        {
                                         dummy[i][j] += p[i+k][j+l];
                                         ++m;
                                        }
                 dummy[i][j] /= m;
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
void GrayImage::MedianFilterNxN (int N)         // odd N
{

cout << "Applying " << N << "x" << N << " median filter...";

int **dummy = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++j)
                {
                 int *list = new int [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 >= height || j+l < 0 
                                        || j+l >= width))
                                        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];

                 delete list;
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
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 = allocateIntegerMatrix(height, width);

int n2 = N/2;

for (int i = 0; i < height; ++i)
        for (int j = 0; j < width; ++j)
                {
                 int *list = new int [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 >= height || j+l < 0 
                                        || j+l >= width))
                                        list[m++] = p[i+k][j+l];
                 dummy[i][j] = findknearest(list, m, k2);

                 delete list;
                }

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

cout << "done!" << endl;

deallocateMatrix(dummy, height, width);
}
//***************************************************************************
static int findknearest(int list[], int m, int k)
{
int center = list[m/2]; // odd N only
mydiff *diff = new mydiff [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;

delete diff;
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();

int **temp = allocateIntegerMatrix(N, N);
int **mask = allocateIntegerMatrix(N, N);
double **out = allocateDoubleMatrix(height, width);
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 < height; i++)
        for (j = 0; j < width; j++)
                out[i][j] = 0.0;

for (i = 0; i < height; i++)
        for (j = 0; j < width; j++)
                {
                 for (k = -N2; k <= N2; k++)
                        for (int l = -N2; l <= N2; l++)
                                {
                                 int ik = i+k;
                                 int jl = j+l;
                                 if (ik >= height)
                                        ik = height-1-(ik-(height-1));
                                 else if (ik < 0)
                                        ik = -ik;
                                 if (jl >= width)
                                        jl = width-1-(jl-(width-1));
                                 else if (jl < 0)
                                        jl = -jl;
                                 out[i][j] += p[ik][jl]*mask[k+N2][l+N2]
                                                *1.0/normfactor;
                                }
                }

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

cout << "done!" << endl;

deallocateMatrix(temp, N, N);
deallocateMatrix(mask, N, N);
deallocateMatrix(out, height, width);
}
//***************************************************************************
int **allocateIntegerMatrix(int m, int n)
{
int **p;

p = new int *[m];
for (int i = 0; i < m; i++)
        {
         p[i] = new int [n];
         if (!p[i])
                {
                 cerr << "allocateIntegerMatrix: Not enough memory!!"
                        << endl;
                 exit(-1);
                }
         for (int j = 0; j < n; j++)
                p[i][j] = 0;
        }

return p;
}
//***************************************************************************
double **allocateDoubleMatrix(int m, int n)
{
double **p;

p = new double *[m];
for (int i = 0; i < m; i++)
        {
         p[i] = new double [n];
         if (!p[i])
                {
                 cerr << "allocateIntegerMatrix: Not enough memory!!"
                        << endl;
                 exit(-1);
                }
         for (int j = 0; j < n; j++)
                p[i][j] = 0.0;
        }

return p;
}
//***************************************************************************
void deallocateMatrix(int **p, int m, int n)
{
for (int i = 0; i < m; i++)
        {
         delete p[i];
        }

delete p;
}
//***************************************************************************
void deallocateMatrix(double **p, int m, int n)
{
for (int i = 0; i < m; i++)
        {
         delete p[i];
        }

delete p;
}
//***************************************************************************