//**************************************************************************
#include <iostream.h>
#include <fstream.h>
#include <math.h>
#include "image.h"
#include "region.h"
//**************************************************************************
#define MAX_NUM_OF_REGIONS 500
#define JOIN_DIST_FACTOR 0.5
static Region R[MAX_NUM_OF_REGIONS];
static int numRegions;
static GrayImage tagImage;
//**************************************************************************
void makeRegions(Matrix &Im, int thresholdsize)
{
Matrix img;
numRegions = 0;

img = Im;
tagImage.Init(img.nx, img.ny, 255);

cout << "Making regions.....";
cout.flush();

for (int i = 0; i < img.nx; i++)
        for (int j = 0; j < img.ny; j++)
                {
                 if (img.isBoundaryPixel(i,j))
                        {
                         R[numRegions].reInitialize();
                         R[numRegions].Insert(i,j, img.p[i][j]);
                         img.p[i][j] = 0.0;
                         tagImage.p[i][j] = (unsigned char) numRegions;
                         followAndfillRegion(&R[numRegions], &img);
                         if (R[numRegions].size() >= thresholdsize)
                                ++numRegions;
                         if (numRegions == MAX_NUM_OF_REGIONS)
                                {
                                 cout << "Max number of regions reached!!" << endl;
                                 // break;
                                 i = img.nx;
                                 j = img.ny;
                                }
                        }
                }

cout << "done (" << numRegions << " regions created)" << endl;
cout.flush();

//printRegionLengths();
}
//**************************************************************************
void followAndfillRegion(Region *R, Matrix *img)
{
Pixel P = (*R)();
int x = P.x;
int y = P.y;
float strength = P.strength;
double err = 1.0e+30;
int im = -1, jm = -1;
int height = img->nx;
int width = img->ny;

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) && img->p[xn][yn])
                        {
                         float errlocal = fabs(strength - img->p[xn][yn]);
                         if (err > errlocal)    // find min err
                                {
                                 err = errlocal;
                                 im = xn;
                                 jm = yn;
                                }
                        }
                }
if (im >=0 && jm >= 0)  // neighbouring pixel found
        {
         R->Insert(im, jm, img->p[im][jm]);
         img->p[im][jm] = 0.0;
         tagImage.p[im][jm] = (unsigned char) numRegions;
         followAndfillRegion(R, img); 
        }
}
//**************************************************************************
void makeRegions(GrayImage &Im, int thresholdsize)
{
Matrix img(Im.height, Im.width);

Im.WriteToMatrix(&img);
makeRegions(img, thresholdsize);
}
//**************************************************************************
void writeAllRegions(GrayImage *img)
{
int height = img->height;
int width = img->width;

img->Zero();
for (int i = 0; i < numRegions; i++)
        {
         for (R[i].begin(); !R[i]; ++R[i])
                {
                 Pixel p = R[i]();
                 img->p[p.x][p.y] = i+1;
                }
        }
}
//**************************************************************************
void joinRegions(Region &r1, Region &r2, Region *final, int r_num)
{
Pixel pr1[2];
Pixel pr2[2];

pr1[0] = r1.head();
pr1[1] = r1.tail();
pr2[0] = r2.head();
pr2[1] = r2.tail();

int mindist = 100000;
int min_i = -1;
for (int i = 0; i < 4; i++)
        {
         int mydist = dist(pr1[i/2], pr2[i%2]);
         if (mydist < mindist)
                {
                 mindist = mydist;
                 min_i = i;
                }
        }

final->reInitialize();
if (min_i/2 == 0)       // Join head of r1
        {
         Pixel *pR = new Pixel [r1.size()];
         int k = 0;
         for (r1.begin(); !r1; ++r1)
                {
                 pR[k] = r1();
                 ++k;
                }
         // insert in reverse order from tail to head
         for (int i = r1.size()-1; i >= 0; i--)
                final->Insert(pR[i]);

         delete pR;
        }
else                    // Join tail of r1
        {
         // insert in normal order from head to tail
         for (r1.begin(); !r1; ++r1)
                final->Insert(r1());
        }

// Fill the gap between the two regions (pixels pr1[min_i/2] and pr2[min_i%2])
insertLine(final, pr1[min_i/2], pr2[min_i%2], r_num);

if (min_i%2 == 0)       // Join head of r2
        {
         // insert in normal order from head to tail
         for (r2.begin(); !r2; ++r2)
                final->Insert(r2());
        }
else                    // Join tail of r2
        {
         Pixel *pR = new Pixel [r2.size()];
         int k = 0;
         for (r2.begin(); !r2; ++r2)
                {
                 pR[k] = r2();
                 ++k;
                }
         // insert in reverse order from tail to head
         for (int i = r2.size()-1; i >= 0; i--)
                final->Insert(pR[i]);

         delete pR;
        }
}
//**************************************************************************
// MidPointLine algorithm from Foley & Van Dam pg 78
void insertLine(Region *r, Pixel &p0, Pixel &p1, int value)
{
int n = dist(p0, p1);
for (int i = 1; i < n; i++)
        {
         double t = 1.0*i/n;
         int x = (int) (p0.x + t * (p1.x - p0.x));
         int y = (int) (p0.y + t * (p1.y - p0.y));
         r->Insert(x, y, (double) value); 
         tagImage.p[x][y] = value;
        }

//cout << "Inserting line between [" << p0.x << "," << p0.y << "] and ["
//      << p1.x << "," << p1.y << "] (" << n << " pixels)" << endl;
}
//**************************************************************************
void joinRegions(void)
{
int k = 0;
for (int i = 0; i < numRegions-1; i++)
        {
         if (R[i].isClosed() || R[i].size() == 0)
                continue;
         int n = FindNearestRegion(R[i], i);
         if (n == -1) 
                 continue;

         int newnum = numRegions+k;
/*
         cout << "Creating region " << newnum+1 << " by joining regions "
                << i+1 << " and " << n+1 << endl;
         cout.flush();
*/

         joinRegions(R[i], R[n], &R[newnum], newnum);
         R[i].reInitialize();
         R[n].reInitialize();
         ++k;
        }
numRegions += k;
compactRegions();
}
//**************************************************************************
int FindNearestRegion(Region &r, int regnum)
{
int mindist = 100000;
int minregnum = -1;

for (int i = regnum+1; i < numRegions; i++)
        {
         if (R[i].isClosed())
                continue;
         int mydist = dist(r, R[i]);     
         if (mydist < mindist)
                {
                 mindist = mydist;
                 minregnum = i;
                }
        }

if (mindist < JOIN_DIST_FACTOR*r.size())
        return minregnum;
else
        return -1;
}
//**************************************************************************
int dist(Region &r1, Region &r2)
{
Pixel pr1[2];
Pixel pr2[2];

pr1[0] = r1.head();
pr1[1] = r1.tail();
pr2[0] = r2.head();
pr2[1] = r2.tail();

int mindist = 100000;
for (int i = 0; i < 4; i++)
        {
         int mydist = dist(pr1[i/2], pr2[i%2]);
         if (mydist < mindist)
                 mindist = mydist;
        }

return mindist;
}
//**************************************************************************
void printRegionLengths(void)
{
int n = 0;
for (int i = 0; i < numRegions; i++)
        {
         if (R[i].size())
                {
                 cout << "Length of region " << i+1 << " = " << R[i].size() << endl;
                 ++n;
                }
        }
cout << "Number of non-zero regions = " << n << endl;
cout.flush();
}
//**************************************************************************
int numNonzeroRegions(void)
{
int n = 0;
for (int i = 0; i < numRegions; i++)
        {
         if (R[i].size())
                ++n;
        }

return n;
}
//**************************************************************************
void compactRegions(void)
{
int n = numNonzeroRegions();
Region *Rnew = new Region [n];

int count = 0;
for (int i = 0; i < numRegions; i++)
        {
         if (R[i].size())
                {
                 Rnew[count] = R[i];
                 ++count;
                 R[i].reInitialize();
                }
        }
numRegions = count;
for (int i = 0; i < numRegions; i++)
        {
         R[i] = Rnew[i];
        }

delete Rnew;
}
//**************************************************************************
Region &Region::operator= (Region &Reg)
{
closed = Reg.closed;
length = 0;

for (Reg.begin(); !Reg; ++Reg)
        Insert(Reg());

if (length != Reg.length)
        {
         cerr << "Region: operator= Error!!" << endl;
         exit(-1);
        }
return *this;
}
//**************************************************************************
void refineRegions(int reg_size, int n)
{
for (int i = 0; i < n; i++)
        joinRegions();

GrayImage dummy;
int height = tagImage.height;
int width = tagImage.width;
dummy.Init(height, width, 255);
writeAllRegions(&dummy);

makeRegions(dummy, 2*reg_size);
joinRegions();
}
//**************************************************************************