// This file contains functions for image processing


#include "../const.h"

#ifdef AIBO_YES
#include <OPENR/OFbkImage.h>
#include <OPENR/OSyslog.h>
#include <OPENR/OPENRAPI.h>
#include <OPENR/core_macro.h>
#endif

#include <stdlib.h>
#include <sys/stat.h>
#include <math.h>
#include "Image.h"
#include "BMP.h"

using namespace std;


// The morphological operator of "closing" is applied to the input image
bool
Close(const Image& image, Image& closed)
{

int x, y;
bool result;
Image dilated;

#ifdef FUNC_TRACK
PRN(("Starting Close\n"));
#endif


result = Dilate(image, dilated);
if (!result)
   return false;
result = Erode(dilated, closed);
if (!result)
   return false;


// Copying first and last row, first and last column of the original image into the closed image
for (x = 0; x < image.width; x++) {
    closed.data[x][image.height - 1] = image.data[x][image.height - 1];
    closed.data[x][0] = image.data[x][0];
}
for (y = image.height - 2; y > 0; y--) {
    closed.data[0][y] = image.data[0][y];
    closed.data[image.width - 1][y] = image.data[image.width - 1][y];
}

FreeImg(dilated);


#ifdef FUNC_TRACK
PRN(("Ending close\n"));
#endif

return true;

}

// The morphological operator of "opening" is applied to the input image
bool
Open(Image& image, Image& opened)
{

int x, y;
bool result;
Image eroded;


#ifdef FUNC_TRACK
PRN(("Starting Open\n"));
#endif

result = Erode(image, eroded);
if (!result)
   return false;
result = Dilate(eroded, opened);
if (!result)
   return false;

// Copying first and last row, first and last column of the original image into the opened image
for (x = 0; x < image.width; x++) {
    opened.data[x][image.height - 1] = image.data[x][image.height - 1];
    opened.data[x][0] = image.data[x][0];
}
for (y = image.height - 2; y > 0; y--) {
    opened.data[0][y] = image.data[0][y];
    opened.data[image.width - 1][y] = image.data[image.width - 1][y];
}

FreeImg(eroded);

#ifdef FUNC_TRACK
PRN(("Ending Open\n"));
#endif

return true;

}

// The morphological operator of "erosion" is applied to the input image (a 3x3 mask is used)
bool
Erode(const Image& image, Image& eroded)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Erosion\n"));
#endif

eroded = Image(image.width, image.height, image.offsetX, image.offsetY);
if (eroded.data == NULL)
   return false;

// Main cycle (first and last row, first and last column not considered)
for (x = 1; x < image.width - 1; x++) {
     for (y = 1; y < image.height - 1; y++) {
            if (image.data[x-1][y-1]     && image.data[x][y-1] && image.data[x+1][y-1]
                && image.data[x-1][y]    && image.data[x][y]   && image.data[x+1][y]
                && image.data[x-1][y+1]  && image.data[x][y+1] && image.data[x+1][y+1]) {
                   
                                eroded.data[x][y] = WHITE;
            }

        }
}


#ifdef FUNC_TRACK
PRN(("Ending Erosion\n"));
#endif

return true;

}

// The morphological operator of "dilation" is applied to the input image (a 3x3 mask is used)
bool
Dilate(const Image& image, Image& dilated)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Dilation\n"));
#endif

dilated = Image(image.width, image.height, image.offsetX, image.offsetY);
if (dilated.data == NULL)
   return false;

// Main cycle (first and last row, first and last column not considered)
for (x = 1; x < image.width - 1; x++) {
     for (y = 1; y < image.height - 1; y++) {
            if (image.data[x][y]) {
                     dilated.data[x][y] = dilated.data[x-1][y-1] = dilated.data[x][y-1]     = dilated.data[x+1][y-1]
                                        = dilated.data[x-1][y]   = dilated.data[x+1][y]     = dilated.data[x-1][y+1]
                                        = dilated.data[x][y+1]   = dilated.data[x+1][y+1]   = WHITE;
            }

        }
}


#ifdef FUNC_TRACK
PRN(("Ending Dilation\n"));
#endif

return true;

}

// This function erases isolated pixels from the input image
bool
EraseSinglePixels(const Image& image, Image& result)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting EraseSinglePixels\n"));
#endif


result =  Image(image.width, image.height, image.offsetX, image.offsetY);
if (result.data == NULL)
   return false;

// Main cycle (first and last row, first and last column not considered)
for (x = 1; x < image.width - 1; x++) {
     for (y = 0; y < image.height - 1; y++) {
            if (image.data[x][y]) {
                     if (!image.data[x-1][y-1] && !image.data[x][y-1] && !image.data[x+1][y-1] &&
                         !image.data[x-1][y] && !image.data[x+1][y] &&
                         !image.data[x-1][y+1] && !image.data[x][y+1] && !image.data[x+1][y+1])
                                               result.data[x][y] = BLACK;
                     else
                                               result.data[x][y] = WHITE;
            }

        }
}



#ifdef FUNC_TRACK
PRN(("Ending EraseSinglePixels\n"));
#endif

return true;

}


// Implementation of the plain algorithm for labelling connected components (8-neighborhood)
// The function returns the labels image, the number of found connected components and 
// the array of used labels.
bool
LabelComponents(const Image& image, Image& image_label, int& number, int*& labels)
{

byte neighbors_label[4];
int dim, label, x, y, i, l;
int* equiv;
bool* found_labels;


#ifdef FUNC_TRACK
PRN(("Starting LabelComponents\n"));
#endif

// Initialization...
image_label = Image(image.width, image.height, image.offsetX, image.offsetY);
if (image_label.data == NULL)
   return false;

dim = NUM_COL;
label = 1;


// Equivalence of labels
equiv  = (int *) calloc(dim, sizeof(int));
for (i = 1; i < dim; i++)
    equiv[i] = i;


// Labels are "propagated" from neighbors and equivalences are stored in the 'equiv' array
for (x = 0; x < image.width; x++) {
     for (y = 0; y < image.height; y++) {
             if (image.data[x][y]) {
                if (y == 0) {

                  if (x == 0)  {

                       image_label.data[x][y] = label;
                       label++;
                   }

                  else {

                       neighbors_label[1] = neighbors_label[4] = 0;
                       neighbors_label[2] = image_label.data[x-1][y];
                       neighbors_label[3] = image_label.data[x-1][y+1];

                       if (neighbors_label[2])  {
                          image_label.data[x][y] = neighbors_label[2];
                          if (neighbors_label[3]) {
                             if (equiv[neighbors_label[3]] < equiv[neighbors_label[2]]) {
                                equiv[equiv[neighbors_label[2]]] = equiv[neighbors_label[2]] = equiv[neighbors_label[3]];
                             }
                             if (equiv[neighbors_label[3]] > equiv[neighbors_label[2]])  {
                                equiv[equiv[neighbors_label[3]]] = equiv[neighbors_label[3]] = equiv[neighbors_label[2]];
                             }
                          }
                       }
                       else if (neighbors_label[3]){
                          image_label.data[x][y] = neighbors_label[3];
                       }
                       else {
                          image_label.data[x][y] = label;
                          label++;
                       }


                   }
              }

              else if (y == (image.height - 1))  {

                  if (x == 0)  {
                               	
                        image_label.data[x][y] = image_label.data[x][y-1];
                        if (! image_label.data[x][y]) {
                            image_label.data[x][y] = label;
                            label++;
                        }

                   }
                   
                  // if (x == (image.width - 1)) {
                   else {                              	

                       neighbors_label[3] = 0;
                       neighbors_label[1] = image_label.data[x-1][y-1];
                       neighbors_label[2] = image_label.data[x-1][y];
                       neighbors_label[4] = image_label.data[x][y-1];

                       if (neighbors_label[1])  {
                          image_label.data[x][y] = neighbors_label[1];
                          if (neighbors_label[2]) {
                             if (equiv[neighbors_label[2]] < equiv[neighbors_label[1]])
                                equiv[equiv[neighbors_label[1]]] = equiv[neighbors_label[1]] = equiv[neighbors_label[2]];
                             if (equiv[neighbors_label[2]] > equiv[neighbors_label[1]])
                                equiv[equiv[neighbors_label[2]]] = equiv[neighbors_label[2]] = equiv[neighbors_label[1]];
                          }
                          if (neighbors_label[4]) {
                             if (equiv[neighbors_label[4]] < equiv[neighbors_label[1]])
                                equiv[equiv[neighbors_label[1]]] = equiv[neighbors_label[1]] = equiv[neighbors_label[4]];
                             if (equiv[neighbors_label[4]] > equiv[neighbors_label[1]])
                                equiv[equiv[neighbors_label[4]]] = equiv[neighbors_label[4]] = equiv[neighbors_label[1]];
                          }
                       }
                       else if (neighbors_label[2]){
                          image_label.data[x][y] = neighbors_label[2];
                          if (neighbors_label[4]) {
                             if (equiv[neighbors_label[4]] < equiv[neighbors_label[2]])
                                equiv[equiv[neighbors_label[2]]] = equiv[neighbors_label[2]] = equiv[neighbors_label[4]];
                             if (equiv[neighbors_label[4]] > equiv[neighbors_label[2]])
                                equiv[equiv[neighbors_label[4]]] = equiv[neighbors_label[4]] = equiv[neighbors_label[2]];
                          }
                       }
                       else if (neighbors_label[4]){
                          image_label.data[x][y] = neighbors_label[4];
                       }
                       else {
                          image_label.data[x][y]  = label;
                          label++;
                       }

                   }


              }
              
              else  {
                     if (x == 0)  {

                       image_label.data[x][y] = image_label.data[x][y-1];
                       if (! image_label.data[x][y]) {
                           image_label.data[x][y] = label;
                           label++;
                       }
                      }
                       else {

                            neighbors_label[1] = image_label.data[x-1][y-1];
                            neighbors_label[2] = image_label.data[x-1][y];
                            neighbors_label[3] = image_label.data[x-1][y+1];
                            neighbors_label[4] = image_label.data[x][y-1];

                           if (neighbors_label[1])  {
                              image_label.data[x][y] = neighbors_label[1];
                              if (neighbors_label[2]) {
                                 if(equiv[neighbors_label[2]] < equiv[neighbors_label[1]])
                                    equiv[equiv[neighbors_label[1]]] = equiv[neighbors_label[1]] = equiv[neighbors_label[2]];
                                 if(equiv[neighbors_label[2]] > equiv[neighbors_label[1]])
                                    equiv[equiv[neighbors_label[2]]] = equiv[neighbors_label[2]] = equiv[neighbors_label[1]];
                              }
                              if (neighbors_label[3]) {
                                 if (equiv[neighbors_label[3]] < equiv[neighbors_label[1]])
                                    equiv[equiv[neighbors_label[1]]] = equiv[neighbors_label[1]] = equiv[neighbors_label[3]];
                                 if (equiv[neighbors_label[3]] > equiv[neighbors_label[1]])
                                    equiv[equiv[neighbors_label[3]]] = equiv[neighbors_label[3]] = equiv[neighbors_label[1]];
                              }
                              if (neighbors_label[4])  {
                                 if (equiv[neighbors_label[4]] < equiv[neighbors_label[1]])
                                    equiv[equiv[neighbors_label[1]]] = equiv[neighbors_label[1]] = equiv[neighbors_label[4]];
                                 if (equiv[neighbors_label[4]] > equiv[neighbors_label[1]])
                                    equiv[equiv[neighbors_label[4]]] = equiv[neighbors_label[4]] = equiv[neighbors_label[1]];
                              }

                           }
                           else if (neighbors_label[2]){
                              image_label.data[x][y] = neighbors_label[2];
                              if (neighbors_label[3]) {
                                 if (equiv[neighbors_label[3]] < equiv[neighbors_label[2]])
                                    equiv[equiv[neighbors_label[2]]] = equiv[neighbors_label[2]] = equiv[neighbors_label[3]];
                                 if (equiv[neighbors_label[3]] > equiv[neighbors_label[2]])
                                    equiv[equiv[neighbors_label[3]]] = equiv[neighbors_label[3]] = equiv[neighbors_label[2]];
                              }
                              if (neighbors_label[4]) {
                                 if (equiv[neighbors_label[4]] < equiv[neighbors_label[2]])
                                    equiv[equiv[neighbors_label[2]]] = equiv[neighbors_label[2]] = equiv[neighbors_label[4]];
                                 if (equiv[neighbors_label[4]] > equiv[neighbors_label[2]])
                                    equiv[equiv[neighbors_label[4]]] = equiv[neighbors_label[4]] = equiv[neighbors_label[2]];
                              }
    
                           }
                           else if (neighbors_label[3]){
                              image_label.data[x][y] = neighbors_label[3];
                              if (neighbors_label[4]) {
                                 if (equiv[neighbors_label[4]] < equiv[neighbors_label[3]])
                                    equiv[equiv[neighbors_label[3]]] = equiv[neighbors_label[3]] = equiv[neighbors_label[4]];
                                 if (equiv[neighbors_label[4]] > equiv[neighbors_label[3]])
                                    equiv[equiv[neighbors_label[4]]] = equiv[neighbors_label[4]] = equiv[neighbors_label[3]];
                              }
    
                           }
                           else if (neighbors_label[4]){
                              image_label.data[x][y] = neighbors_label[4];
                           }
                           else {
                              image_label.data[x][y] = label;
                              label++;
                           }
                   }
              }
          }
          if (label >= dim) {         // Our image of labels cannot store more than 256 labels!
             PRN(("Error: too many connected components in input image!\n Quitting...\n"));
             return false;
          }
    }
}


equiv = (int *) realloc(equiv, label * sizeof(int));
found_labels = (bool*) calloc(label, sizeof(bool));
labels = (int*) calloc(label, sizeof(int));


// Replacing all equivalent labels with the same representative,
// getting labels and the number of found connected components
number = 0;
i = 0;


for (x = 0; x < image.width; x++) {
     for (y = 0; y < image.height; y++) {
             l = image_label.data[x][y];
             if (l) {
                while (equiv[l] != l) {
                     l = equiv[l];
                }
                image_label.data[x][y] = l;
                if (!found_labels[l]) {
                   found_labels[l] = true;
                   labels[i] = l;
                   number++;
                   i++;
                }

             }
    }
}

labels = (int *) realloc(labels, number * sizeof(int));

free(equiv);
free(found_labels);


#ifdef FUNC_TRACK
PRN(("Ending LabelComponents\n"));
#endif

return true;

}

// This function returns a binary image where only the component labelled with the given label is shown.
bool
GetComponent(const Image& image_label, int label, Image& component)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting GetComponent\n"));
#endif

component = Image(image_label.width, image_label.height, image_label.offsetX, image_label.offsetY);
if (component.data == NULL || label == 0)
   return false;


for (x = 0; x < image_label.width; x++) {
     for (y = 0; y < image_label.height; y++) {
               if (image_label.data[x][y] == label)
                   component.data[x][y] = WHITE;
               else
                   component.data[x][y] = BLACK;
        }
}

#ifdef FUNC_TRACK
PRN(("Ending GetComponent\n"));
#endif

return true;

}

// This function returns a graylevel image where only the component labelled with the given label is shown (with original graylevels)
bool
GetComponentGL(const Image& image_label, const Image& original, Image& component, int label)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting GetComponentGL\n"));
#endif

component = Image(image_label.width, image_label.height, image_label.offsetX, image_label.offsetY);
if (component.data == NULL || label == 0)
   return false;

for (x = 0; x < image_label.width; x++) {
     for (y = 0; y < image_label.height; y++) {
               if (image_label.data[x][y] == label)
                   component.data[x][y] = original.data[x][y];
               else
                   component.data[x][y] = BLACK;
        }
}


#ifdef FUNC_TRACK
PRN(("Ending GetComponentGL\n"));
#endif

return true;

}

// This function returns a binary image where only the greatest connected component in the label image is shown.
int
GetMajorComp(const Image& label_image)
{

int x, y, l, max, label_max, dim;
int* comp;

#ifdef FUNC_TRACK
PRN(("Starting GetMajorComp\n"));
#endif

dim = NUM_COL;

// Array counting the number of pixels for each label
comp = (int*) calloc(dim, sizeof(int));
if (comp == NULL)
   return 0;


for (x = 0; x < label_image.width; x++) {
       for (y = 0; y < label_image.height; y++) {
               l = label_image.data[x][y];
               if (l) {
                  comp[l]++;
                }
         }
}

max = 0;
for (l = 0; l < dim; l++) {
       if (comp[l] > max) {
           max = comp[l];
           label_max = l;
       }
}  

free(comp);

#ifdef FUNC_TRACK
PRN(("Ending GetMajorComp\n"));
#endif

return label_max;

}


// This function returns a binary image (whose dimension is the same as 'rif') where the two input components are shown together
bool
MergeComp(const Image& comp1, const Image& comp2, const Image& rif, Image& merge)
{

int x, y, ox, oy;

#ifdef FUNC_TRACK
PRN(("Starting MergeComp\n"));
#endif

merge = Image(rif.width, rif.height, rif.offsetX, rif.offsetY);
if (merge.data == NULL)
   return false;


ox = comp1.offsetX - rif.offsetX;
oy = comp1.offsetY - rif.offsetY;

for (x = 0; x < comp1.width; x++) {
       for (y = 0 ; y < comp1.height; y++) {
           if (comp1.data[x][y])
               merge.data[x + ox][y + oy] = WHITE;

      }
}

ox = comp2.offsetX - rif.offsetX;
oy = comp2.offsetY - rif.offsetY;

for (x = 0; x < comp2.width; x++) {
       for (y = 0 ; y < comp2.height; y++) {
            if (comp2.data[x][y])
               merge.data[x + ox][y + oy] = WHITE;

      }
}

#ifdef FUNC_TRACK
PRN(("Ending MergeComp\n"));
#endif

return true;

}


// Computing the area (i.e., number of pixels) of the connected component represented in the input binary (or graylevel) image
int
Area(const Image& image)
{

int x, y, area;

#ifdef FUNC_TRACK
PRN(("Starting Area\n"));
#endif

area = 0;

for (x = 0; x < image.width; x++) {
       for (y = 0 ; y < image.height; y++) {
               if (image.data[x][y])
                   area++;
        }
}

#ifdef FUNC_TRACK
PRN(("Ending Area\n"));
#endif

return area;

}

// Computing the x coordinate of the centroid of the connected component represented in the input binary image
bool
CentroidX(const Image& image, int& cx)
{

int x, y, area;

#ifdef FUNC_TRACK
PRN(("Starting CentroidX\n"));
#endif

cx = 0;

area = Area(image);
if (!area)
   return false;

for (x = 0; x < image.width; x++) {
       for (y = 0 ; y < image.height; y++) {
               if (image.data[x][y])
                   cx += x;
        }
}

cx = cx / area;

#ifdef FUNC_TRACK
PRN(("Ending CentroidX\n"));
#endif

return true;

}

// Computing the y coordinate of the centroid of the connected component represented in the input binary image
bool
CentroidY(const Image& image, int& cy)
{

int x, y, area;

#ifdef FUNC_TRACK
PRN(("Starting CentroidY\n"));
#endif

cy = 0;

area = Area(image);
if (!area)
   return false;

for (x = 0; x < image.width; x++) {
       for (y = 0 ; y < image.height; y++) {
               if (image.data[x][y])
                   cy += y;
        }
}

cy = cy / area;

#ifdef FUNC_TRACK
PRN(("Ending CentroidY\n"));
#endif

return true;

}


// Computing the x and y coordinate of the centroid of the connected component represented in the input binary image
bool
Centroid(const Image& image, int*& centroid)
{

int x, y, cx, cy, area;


#ifdef FUNC_TRACK
PRN(("Starting Centroid\n"));
#endif

cx = 0;
cy = 0;

area = Area(image);
if (!area)
   return false;
   
centroid = (int*) calloc(2, sizeof(int));
if (centroid == NULL)
   return false;   

for (x = 0; x < image.width; x++) {
       for (y = 0 ; y < image.height; y++) {
               if (image.data[x][y]) {
                   cx += x;
                   cy += y;
               }
        }
}

cx = cx / area;
cy = cy / area;

// we want the centroid coordinates to refer to the original image the considered component belongs to
centroid[0] = cx + image.offsetX;
centroid[1] = cy + image.offsetY;


#ifdef FUNC_TRACK
PRN(("Ending Centroid\n"));
#endif

return true;

}


// Computing the x and y coordinate of the weighted (by graylevels) centroid 
// of the connected component represented in the input graylevel image.
bool
WeightedCentroid(const Image& image, int*& centroid)
{

int x, y, cx, cy, total;


#ifdef FUNC_TRACK
PRN(("Starting Centroid\n"));
#endif

cx = 0;
cy = 0;
total = 0;

centroid = (int*) calloc(2, sizeof(int));
if (centroid == NULL)
   return false;

for (x = 0; x < image.width; x++) {
       for (y = 0 ; y < image.height; y++) {
               if (image.data[x][y]) {
                   cx += x * image.data[x][y];
                   cy += y * image.data[x][y];
                   total += image.data[x][y];
               }
        }
}

if (!total)
   return false;

cx = cx / (total);
cy = cy / (total);

// we want the centroid coordinates to refer to the original image the considered component belongs to
centroid[0] = cx + image.offsetX;
centroid[1] = cy + image.offsetY;


#ifdef FUNC_TRACK
PRN(("Ending Centroid\n"));
#endif

return true;

}


// Computing the bounding box of the connected component represented in the input binary image.
// The function returns an array of 4 integers: the maximum and minimum x and y coordinates of the component.
// If mode = 0, coordinates refer to the current image, if mode = 1, they refer to the original image
bool
BoundingBox(const Image& image, int mode, int*& extremes)
{

int x, y, i, x_max, x_min, y_max, y_min;

#ifdef FUNC_TRACK
PRN(("Starting BoundingBox\n"));
#endif

// Initializing...
extremes = (int*) calloc(4, sizeof(int));
if (extremes == NULL)
   return false;

y_max = x_max = 0;
x_min = image.width - 1;
y_min = image.height - 1;

// Computing bounding box
for (x = 0; x < image.width; x++) {
       for (y = 0 ; y < image.height; y++) {
               if (image.data[x][y]) {
                     if (x < x_min)
                        x_min = x;
                     if (x > x_max)
                        x_max = x;
                     if (y < y_min)
                        y_min = y;
                     if (y > y_max)
                        y_max = y;
                }
   }
}

extremes[0] = x_min;
extremes[1] = x_max;
extremes[2] = y_min;
extremes[3] = y_max;

// referrring to the original image...
if (mode) {
   extremes[0] += image.offsetX;
   extremes[1] += image.offsetX;
   extremes[2] += image.offsetY;
   extremes[3] += image.offsetY;
}


#ifdef FUNC_TRACK
PRN(("Ending BoundingBox\n"));
#endif

return true;

}


// This function extracts the background from the input image
bool
Background(const Image& image, Image& background)
{

int x, y, l, dim, dummy_number;
int* dummy_labels;
bool result;
bool* labels;
Image image_label, comp, temp;


#ifdef FUNC_TRACK
PRN(("Starting Background\n"));
#endif

// Initializing...
background = Image(image.width, image.height, image.offsetX, image.offsetY);
if (background.data == NULL)
   return false;

dim = (image.height * 2) + (image.width - 2) * 2 + 1;    // number of pixels on the borders
labels = (bool*) calloc(dim, sizeof(bool));

result = LabelComponents(image, image_label, dummy_number, dummy_labels);
if (!result) {
   return false;
}

// The background is the sum of all connected components adjacent to the border:
// we're looking for components labelled by the same labels we found for pixels
// placed in the first and last row, first and last column.
for (y = image_label.height - 1 ; y >= 0; y--) {
         l = image_label.data[0][y];
         if (l){
            if (!labels[l]) {
               labels[l] = true;
               result = GetComponent(image_label, l, comp);
               if (!result)
                  return false;

               result = CopyImage(background, temp);
               FreeImg(background);
               MergeComp(temp, comp, image, background);
               FreeImg(comp);
               FreeImg(temp);
            }
         }
         l = image_label.data[image_label.width - 1][y];
         if (l){
            if (!labels[l]) {
               labels[l] = true;
               result = GetComponent(image_label, l, comp);
               if (!result)
                  return false;
               result = CopyImage(background, temp);
               FreeImg(background);
               MergeComp(temp, comp, image, background);
               FreeImg(comp);
               FreeImg(temp);
            }
         }
}

for (x = 1 ; x < image_label.width - 1; x++) {
         l = image_label.data[x][0];
         if (l){
            if (!labels[l]) {
               labels[l] = true;
               result = GetComponent(image_label, l, comp);
               if (!result)
                  return false;

               result = CopyImage(background, temp);
               FreeImg(background);
               MergeComp(temp, comp, image, background);
               FreeImg(comp);
               FreeImg(temp);
            }
         }
         l = image_label.data[x][image_label.height - 1];
         if (l){
            if (!labels[l]) {
               labels[l] = true;
               result = GetComponent(image_label, l, comp);
               if (!result)
                  return false;

               result = CopyImage(background, temp);
               FreeImg(background);
               MergeComp(temp, comp, image, background);
               FreeImg(comp);
               FreeImg(temp);
            }
         }
}

// Freeing memory...
free(dummy_labels);
free(labels);
FreeImg(image_label);

#ifdef FUNC_TRACK
PRN(("Ending Background\n"));
#endif

return true;

}

// Binary And: only white pixels in both input images will be white in the output image
bool
And(const Image& im1, const Image& im2, Image& result)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting And\n"));
#endif

result = Image(im1.width, im1.height, im1.offsetX, im1.offsetY);
if (result.data == NULL)
   return false;

for (x = 0; x < result.width; x++) {
      for (y = 0 ; y < result.height; y++) {
                 if (im1.data[x][y] && im2.data[x][y])
                        result.data[x][y] = WHITE;
        }
}


#ifdef FUNC_TRACK
PRN(("Ending And\n"));
#endif

return true;

}

// This function returns the negative of the input binary image
bool
Negative(const Image& image, Image& neg)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Negative\n"));
#endif

neg = Image(image.width, image.height, image.offsetX, image.offsetY);
if (neg.data == NULL)
   return false;


for (x = 0; x < neg.width; x++) {
      for (y = 0; y  < neg.height; y++) {
              neg.data[x][y] = WHITE - image.data[x][y];
       }
}


#ifdef FUNC_TRACK
PRN(("Ending Negative\n"));
#endif

return true;

}

// This function implements the mean filter with a 3x3 mask
bool
Mean(const Image& image, Image& filtered)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Mean Filter\n"));
#endif

filtered = Image(image.width, image.height, image.offsetX, image.offsetY);
if (filtered.data == NULL)
   return false;


for (x = 1; x < image.width - 1; x++) {
       for (y = 1; y < image.height - 1; y++) {
               filtered.data[x][y] = ((image.data[x-1][y-1]) + (image.data[x][y-1]) + (image.data[x+1][y-1]) +
                                     (image.data[x-1][y])   + (image.data[x][y])   + (image.data[x+1][y]) +
                                     (image.data[x-1][y+1]) + (image.data[x][y+1]) + (image.data[x+1][y+1])) / 9;
        }
}

// reconstructing border...
for (y = filtered.height - 1 ; y >= 0; y--) {
    filtered.data[0][y] = image.data[0][y];
    filtered.data[filtered.width - 1][y] = image.data[filtered.width - 1][y];
}
for (x = 1; x < filtered.width - 1; x++) {
    filtered.data[x][0] = image.data[x][0];
    filtered.data[x][filtered.height - 1] = image.data[x][filtered.height - 1];
}


#ifdef FUNC_TRACK
PRN(("Ending Mean Filter\n"));
#endif

return true;

}


// This function performs edge detection in the y direction, applying the 3x3 Sobel Filter
bool
SobelX(const Image& image, Image& sobel_image)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting SobelX\n"));
#endif

// Filtered values are scaled to be in the [0, 255] interval.
// Negative values are included in the [0, 127] interval, whereas positive values are represented in [128, 255]

sobel_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (sobel_image.data == NULL)
   return false;


for (x = 1; x < image.width - 1; x++) {
      for (y = 1; y < image.height - 1; y++) {
               sobel_image.data[x][y] = (-1 * (image.data[x-1][y-1]) + 1 * (image.data[x+1][y-1]) +
                                        -2 * (image.data[x-1][y])   + 2 * (image.data[x+1][y]) +
                                        -1 * (image.data[x-1][y+1]) + 1 * (image.data[x+1][y+1]) + SOBEL_SUM) / SOBEL_DIV;
        }
}


#ifdef FUNC_TRACK
PRN(("Ending SobelX\n"));
#endif

return true;

}


// This function performs edge detection in the x direction, applying the 3x3 Sobel Filter
bool
SobelY(const Image& image, Image& sobel_image)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting SobelY\n"));
#endif

// Filtered values are scaled to be in the [0, 255] interval.
// Negative values are included in the [0, 127] interval, whereas positive values are represented in [128, 255]

sobel_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (sobel_image.data == NULL)
   return false;

for (x = 1; x < image.width - 1; x++) {
      for (y = 1; y < image.height - 1; y++) {
               sobel_image.data[x][y] = (1 * (image.data[x-1][y-1]) +   2 * (image.data[x][y-1]) +   1 * (image.data[x+1][y-1]) +
                                        -1 * (image.data[x-1][y+1]) +  -2 * (image.data[x][y+1]) +  -1 * (image.data[x+1][y+1]) + SOBEL_SUM) / SOBEL_DIV;

        }
}


#ifdef FUNC_TRACK
PRN(("Ending SobelY\n"));
#endif

return true;

}


// This function performs edge detection in all directions, using the gradient module (|SobelX| + |SobelY|)
bool
GradientImage(const Image& image, Image& gradient_image)
{

int x, y, value;
Image sx, sy;

#ifdef FUNC_TRACK
PRN(("Starting GradientImage\n"));
#endif

gradient_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (gradient_image.data == NULL)
   return false;

SobelX(image, sx);
SobelY(image, sy);

// Note: we need to transform the [0, 255] interval into [-127, 128] in order to consider only
// the strength of the edge and not its direction
for (x = 1; x < image.width - 1; x++) {
      for (y = 1; y < image.height - 1; y++) {
               value = abs(sx.data[x][y] - 127) + abs(sy.data[x][y] - 127);
               if (value > WHITE)
                  gradient_image.data[x][y] = WHITE;
               else
                  gradient_image.data[x][y] = value;
        }
}

FreeImg(sx);
FreeImg(sy);

#ifdef FUNC_TRACK
PRN(("Ending GradientImage\n"));
#endif

return true;

}

// This function performs edge detection in all directions, using the gradient module (|SobelX| + |SobelY|),
// and returns the gradient direction for all pixels.
bool
GradientImageWithDirection(const Image& image, Image& gradient_image, double**& directions)
{

int x, y, i, value;
Image sx, sy;

#ifdef FUNC_TRACK
PRN(("Starting GradientImageWithDirection\n"));
#endif

// Initialization...

gradient_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (gradient_image.data == NULL)
   return false;
   
directions = (double**) calloc(image.width, sizeof(double*));
if (directions == NULL)
   return false;
for (i = 0; i < image.width; i++) {
   directions[i] = (double*) calloc(image.height, sizeof(double));
   if (directions[i] == NULL)
      return false;
}


SobelX(image, sx);
SobelY(image, sy);

// Note: we need to transform the [0, 255] interval into [-127, 128] in order to consider only
// the strength of the edge and not its direction
for (x = 1; x < image.width - 1; x++) {
      for (y = 1; y < image.height - 1; y++) {
               value = abs(sx.data[x][y] - 127) + abs(sy.data[x][y] - 127);
               if (value > WHITE)
                  gradient_image.data[x][y] = WHITE;
               else
                  gradient_image.data[x][y] = value;
               directions[x][y] = atan2((sy.data[x][y] - 127), (sx.data[x][y] - 127));
               if (!ceilf(directions[x][y]) && (sy.data[x][y] - 127) < 0)
                    directions[x][y] = PI;
        }
}

FreeImg(sx);
FreeImg(sy);

#ifdef FUNC_TRACK
PRN(("Ending GradientImageWithDirection\n"));
#endif

return true;

}


// This function thresholds the input image (usually an edge image) showing only
// those pixels whose graylevel is in the interval [th_min, th_max].
bool
Threshold(const Image& image, int th_min, int th_max, Image& th_image)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Threshold\n"));
#endif

th_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (th_image.data == NULL)
   return false;

for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
               if ((image.data[x][y] >= th_min) && (image.data[x][y] <= th_max))
                  th_image.data[x][y] = WHITE;
               else
                  th_image.data[x][y] = BLACK;
        }
}

#ifdef FUNC_TRACK
PRN(("Ending Threshold\n"));
#endif

return true;

}


// This function thresholds the input image (usually an edge image) showing only
// the 'perc'% pixels having the highest graylevels.
bool
ThresholdPerc(const Image& image, double perc, Image& th_image)
{

int x, y, i, min, total, count;
int* hist;

#ifdef FUNC_TRACK
PRN(("Starting ThresholdPerc\n"));
#endif

th_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (th_image.data == NULL)
   return false;

if (!Histogram(image, hist))
   return false;

total = (int) (perc * (th_image.width * th_image.height));

count = 0;
i = WHITE;
while (count < total) {
      count += hist[i];
      i--;
}
min = i + 1;

for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
               if (image.data[x][y] >= min)
                  th_image.data[x][y] = WHITE;
               else
                  th_image.data[x][y] = BLACK;
        }
}

free(hist);

#ifdef FUNC_TRACK
PRN(("Ending ThresholdPerc\n"));
#endif

return true;

}

// This function thresholds the input image (usually an edge image) showing only
// the 'perc'% pixels having the highest graylevels.
bool
ThresholdPerc(const Image& image, Image& th_image, double perc, int starting_gl, int ending_gl)
{

int x, y, i, min, total, count;
int* hist;

#ifdef FUNC_TRACK
PRN(("Starting ThresholdPerc\n"));
#endif

th_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (th_image.data == NULL)    
   return false;

if (!HistogramWithoutBorder(image, hist))
   return false;

// If we consider 0 as the maximum value for an edge, we have
// to exclude from our computation the border pixel
// (which are manually set to 0)
if (starting_gl == 0)
   total = (int) (perc * ((th_image.width - 2) * (th_image.height - 2)));
else
   total = (int) (perc * (th_image.width * th_image.height));


count = 0;
i = starting_gl;

if (starting_gl > ending_gl) {
      while (count < total) {
            count += hist[i];
            i--;
      
      }

      min = i + 1;

      for (x = 0; x < image.width; x++) {
            for (y = 0; y < image.height; y++) {
                     if (image.data[x][y] >= min)
                        th_image.data[x][y] = WHITE;
                     else
                        th_image.data[x][y] = BLACK;
              }
      }
}
else  {
      while (count < total) {
            count += hist[i];
            i++;
      
      }

      min = i - 1;

      for (x = 0; x < image.width; x++) {
            for (y = 0; y < image.height; y++) {
                     if (image.data[x][y] <= min)
                        th_image.data[x][y] = WHITE;
                     else
                        th_image.data[x][y] = BLACK;
              }
      }
}

free(hist);

#ifdef FUNC_TRACK
PRN(("Ending ThresholdPerc\n"));
#endif

return true;

}

// As the previous function, except it returns also the value used as a threshold
bool
ThresholdPerc(const Image& image, Image& th_image, double perc, int starting_gl, int ending_gl, int& th)
{

int x, y, i, min, total, count;
int* hist;

#ifdef FUNC_TRACK
PRN(("Starting ThresholdPerc\n"));
#endif

th_image = Image(image.width, image.height, image.offsetX, image.offsetY);
if (th_image.data == NULL)    
   return false;

if (!HistogramWithoutBorder(image, hist))
   return false;

// If we consider 0 as the maximum value for an edge, we have
// to exclude from our computation the border pixel
// (which are manually set to 0)
if (starting_gl == 0)
   total = (int) (perc * ((th_image.width - 2) * (th_image.height - 2)));
else
   total = (int) (perc * (th_image.width * th_image.height));


count = 0;
i = starting_gl;

if (starting_gl > ending_gl) {
      while (count < total) {
            count += hist[i];
            i--;
      
      }

      min = i + 1;

      for (x = 0; x < image.width; x++) {
            for (y = 0; y < image.height; y++) {
                     if (image.data[x][y] >= min)
                        th_image.data[x][y] = WHITE;
                     else
                        th_image.data[x][y] = BLACK;
              }
      }
}
else  {
      while (count < total) {
            count += hist[i];
            i++;
      
      }

      min = i - 1;

      for (x = 0; x < image.width; x++) {
            for (y = 0; y < image.height; y++) {
                     if (image.data[x][y] <= min)
                        th_image.data[x][y] = WHITE;
                     else
                        th_image.data[x][y] = BLACK;
              }
      }
}

free(hist);
th = min;

#ifdef FUNC_TRACK
PRN(("Ending ThresholdPerc\n"));
#endif

return true;

}


// This function returns, for each row of the input image, the average graylevel of pixels on that row
bool
HorProjection(const Image& image, int*& proj)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting HorProjection\n"));
#endif

// Array representing horizontal projection for each row
proj = (int*) calloc(image.height, sizeof(int));
if (proj == NULL || image.width == 0)
   return false;


for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
            proj[y] += image.data[x][y];
        }
}


for (y = image.height - 1 ; y >= 0; y--) {
          proj[y] = (int) (proj[y]/image.width);
}

#ifdef FUNC_TRACK
PRN(("Ending HorProjection\n"));
#endif

return true;

}


// This function returns, for each column of the input image, the average graylevel of pixels on that column
bool
VerProjection(const Image& image, int*& proj)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting VerProjection\n"));
#endif

// Array representing vertical projection for each row
proj = (int*) calloc(image.width, sizeof(int));
if (proj == NULL || image.height == 0)
   return false;

for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
                  proj[x] += image.data[x][y];
        }
}


for (x = 0; x < image.width; x++) {
          proj[x] = (int) (proj[x]/image.height);
}

#ifdef FUNC_TRACK
PRN(("Ending VerProjection\n"));
#endif

return true;

}


// This function returns, for each row of the input binary image, the number of white pixels on that row
bool
BinaryHorProjection(const Image& image, int*& proj)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting BinaryHorProjection\n"));
#endif


// Array representing horizontal projection for each row
proj = (int*) calloc(image.height, sizeof(int));
if (proj == NULL)
   return false;


for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
               if (image.data[x][y])
                  proj[y]++;
        }
}

#ifdef FUNC_TRACK
PRN(("Ending BinaryHorProjection\n"));
#endif

return true;

}


// This function returns, for each column of the input binary image, the number of white pixels on that column
bool
BinaryVerProjection(const Image& image, int*& proj)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting BinaryVerProjection\n"));
#endif


// Array representing vertical projection for each row
proj = (int*) calloc(image.width, sizeof(int));
if (proj == NULL)
   return false;


for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
               if (image.data[x][y])
                  proj[x]++;
        }
}

#ifdef FUNC_TRACK
PRN(("Ending BinaryVerProjection\n"));
#endif

return true;

}


// This funcion computes the histogram of input image assuming that pixel values are all in the interval [0, 255]
bool
Histogram(const Image& image, int*& histogram)
{
int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Histogram\n"));
#endif

histogram = (int*) calloc(NUM_COL, sizeof(int));
if (histogram == NULL)
   return false;

for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
              histogram[image.data[x][y]]++;

       }
}

#ifdef FUNC_TRACK
PRN(("Ending Histogram\n"));
#endif

return true;

}

// This funcion computes the histogram of input image assuming that pixel values are all in the interval [0, 255],
// ignoring the image borders (this is useful when considering images whose borders have been manually set to 0)
bool
HistogramWithoutBorder(const Image& image, int*& histogram)
{
int x, y;

#ifdef FUNC_TRACK
PRN(("Starting Histogram\n"));
#endif

histogram = (int*) calloc(NUM_COL, sizeof(int));
if (histogram == NULL)
   return false;

for (x = 1; x < image.width - 1; x++) {
      for (y = 1; y < image.height - 1; y++) {
              histogram[image.data[x][y]]++;

       }
}

#ifdef FUNC_TRACK
PRN(("Ending Histogram\n"));
#endif

return true;

}

// This funcion computes the brigthness of the input image as the average graylevel.
int
Brightness(const Image& image)
{
int x, y, b;

#ifdef FUNC_TRACK
PRN(("Starting Brightness\n"));
#endif

b = 0;
for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
              b += image.data[x][y];
       }
}

b = (int) (b / (image.width * image.height));

#ifdef FUNC_TRACK
PRN(("Ending Brightness\n"));
#endif

return b;

}

// This funcion computes the constrast of the input image.
int
Contrast(const Image& image)
{
int x, y, b, c;

#ifdef FUNC_TRACK
PRN(("Starting Brightness\n"));
#endif

b = Brightness(image);
c = 0;
for (x = 0; x < image.width; x++) {
      for (y = 0; y < image.height; y++) {
              c += (image.data[x][y] - b) * (image.data[x][y] - b);
       }
}

c = (int) sqrt((c / (image.width * image.height)));

#ifdef FUNC_TRACK
PRN(("Ending Brightness\n"));
#endif

return c;

}



// This function computes the bounding box of the face ('ext') in the image and, simultaneously, crops
// the clustered image according to the bounding box, in order to get the face image; 
// also, it returns the skin_map centroid.
bool
FindFace(const Image& clustered, Image& skin_map, int*& ext, int*& centroid, const Init_para& parameters, int serialNumber)
{

int l, dummy_number;
int* dummy_labels;
bool result;
Image im_label, closed, face;

#ifdef FUNC_TRACK
PRN(("Starting FindFace\n"));
#endif

// Closing clustered image
result = Close(clustered, closed);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// Labeling components and getting face region
result = LabelComponents(closed, im_label, dummy_number, dummy_labels);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

l = GetMajorComp(im_label);
result = GetComponent(im_label, l, face);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// We found a face only if the component's area is at most perc_face% of the picture's dimension
if (Area(face) < (parameters.perc_face * clustered.width * clustered.height)) {
   PRN(("Face not found!\n"));
   return false;
}

result = BoundingBox(face, 0, ext);
if (!result)
   return false;
result = Centroid(face, centroid);
if (!result)
   return false;

// Cropping face
result = CropImage(ext, clustered, skin_map, 0);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// Freeing memory...
free(dummy_labels);
FreeImg(closed);
FreeImg(im_label);
FreeImg(face);

#ifdef FUNC_TRACK
PRN(("Ending FindFace\n"));
#endif

return true;

}

// This function returns the centroid of the face in the image
bool
FindFaceCentroid(const Image& clustered, int*& centroid, const Init_para& parameters, int serialNumber)
{

int l, dummy_number;
int* dummy_labels;
bool result;
Image temp, closed, im_label;

#ifdef FUNC_TRACK
PRN(("Starting FindFaceCentroid\n"));
#endif

// Closing clustered image
result = Close(clustered, closed);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// Labeling components and getting face region
result = LabelComponents(closed, im_label, dummy_number, dummy_labels);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

l = GetMajorComp(im_label);
result = GetComponent(im_label, l, temp);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// We found a face only if the component's area is at most perc_face% of the picture's dimension
if (Area(temp) < (parameters.perc_face * clustered.width * clustered.height)) {
   PRN(("Face not found!\n"));
   return false;
}


if (!Centroid(temp, centroid))
   return false;

free(dummy_labels);

#ifdef FUNC_TRACK
PRN(("Ending FindFaceCentroid\n"));
#endif

return true;

}

// This function returns the centroid of face and mouth in the image
bool
FindFaceAndMouthCentroid(const Image& clustered, const Image& cb, const Image& cr, int*& face_centroid, int*& mouth_centroid, const Init_para& parameters, int serialNumber)
{

int l, dummy_number;
int* dummy_labels;
int* face_extr;
bool result;
Image temp, closed, im_label;

#ifdef FUNC_TRACK
PRN(("Starting FindFaceAndMouthCentroid\n"));
#endif

// Closing clustered image
result = Close(clustered, closed);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// Labeling components and getting face region
result = LabelComponents(closed, im_label, dummy_number, dummy_labels);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

l = GetMajorComp(im_label);
result = GetComponent(im_label, l, temp);
if (!result) {
   PRN(("Face not found!\n"));
   return false;
}

// We found a face only if the component's area is at most perc_face% of the picture's dimension
if (Area(temp) < (parameters.perc_face * clustered.width * clustered.height)) {
   PRN(("Face not found!\n"));
   return false;
}

// Getting face centroid
if (!Centroid(temp, face_centroid))
   return false;

// Finding mouth
BoundingBox(temp, 0, face_extr);
result = FindMouth(face_extr, cb, cr, temp, ((face_extr[1] - face_extr[0] + 1) * (face_extr[3] - face_extr[2] + 1)), parameters, serialNumber);
if (!result)
   return false;

// Getting mouth centroid
if (!Centroid(temp, mouth_centroid))
   return false;

// Freeing memory
free(dummy_labels);
free(face_extr);
FreeImg(temp);
FreeImg(closed);
FreeImg(im_label);

#ifdef FUNC_TRACK
PRN(("Ending FindFaceAndMouthCentroid\n"));
#endif

return true;

}

// This function performs a transformation involving Cb and Cr values to highlight lips pixels.
bool
MouthMap(const Image& cb, const Image& cr, Image& mouth, int dim, const Init_para& parameters, int serialNumber)
{

int x, y, i, total, cont, min;
int* hist;

#ifdef FUNC_TRACK
PRN(("Starting MouthMap\n"));
#endif

// Initializing...
mouth = Image(cr.width, cr.height, cr.offsetX, cr.offsetY);
if (mouth.data == NULL)
   return false;

hist = (int*) calloc(NUM_COL, sizeof(int));
if (hist == NULL)
   return false;

// Applying a transformation to highlight probable lips pixels and building the so obtained image histogram
for (x = 0; x < mouth.width; x++) {
     for (y = 0; y < mouth.height; y++) {
                mouth.data[x][y] = ((WHITE - (cr.data[x][y] - (cb.data[x][y] * cb.data[x][y]))) * cr.data[x][y] * cr.data[x][y])/MOUTH_SCALE;
                hist[mouth.data[x][y]]++;
       }
}

// Debug: writing to file the mouth map
if (parameters.debug_mouth) {
   char name[MAX_LENGTH_FILE];
   BMP bmp;

   sprintf(name, "%s%sMM%d.BMP", DATA_FOLDER, parameters.file_debug_mouth, serialNumber);
   bmp.Write2File(name, mouth);
}

// We set to 1 the PERC_MOUTH% of pixels having the greatest graylevels...
total = (int) (parameters.perc_mouth * dim);   // PERC_MOUTH of total pixels in face region
cont = 0;
min = 0;

i = WHITE;
while (cont < total) {
      cont += hist[i];
      i--;
}
min = i;

// only pixels whose value is above 'min' are shown in the output image
for (x = 0; x < mouth.width; x++) {
     for (y = 0; y < mouth.height; y++) {
                if (mouth.data[x][y] > min)
                   mouth.data[x][y] = WHITE;
                else
                   mouth.data[x][y] = BLACK;
       }
}

free(hist);

#ifdef FUNC_TRACK
PRN(("Ending MouthMap\n"));
#endif

return true;

}

// This function returns the binary image of the mouth
bool
FindMouth(int* ext, const Image& cb, const Image& cr, Image& mouth, int dim, const Init_para& parameters, int serialNumber)
{

int l, dummy_number;
int* dummy_labels;
bool result;
Image cb_cropped, cr_cropped, 
      im_label, mouth_map;

#ifdef FUNC_TRACK
PRN(("Starting FindMouth\n"));
#endif

// Getting mouth map
result = CropImage(ext, cb, cb_cropped, NO_BORDER);
if (!result)
   return false;
result = CropImage(ext, cr, cr_cropped, NO_BORDER);
if (!result)
   return false;

result = MouthMap(cb_cropped, cr_cropped, mouth_map, dim, parameters, serialNumber);
if (!result) {
   PRN(("Mouth not found!\n"));
   return result;
}

// Debug: writing to file the mouth map
if (parameters.debug_mouth) {
   char name[MAX_LENGTH_FILE];
   BMP bmp;

   sprintf(name, "%s%s%d.BMP", DATA_FOLDER, parameters.file_debug_mouth, serialNumber);
   bmp.Write2File(name, mouth_map);
}


// Labeling components and getting mouth region, as the greatest component in the mouth map
result = LabelComponents(mouth_map, im_label, dummy_number, dummy_labels);
if (!result) {
   PRN(("Mouth not found!\n"));
   return result;
}
l = GetMajorComp(im_label);
result = GetComponent(im_label, l, mouth);
if (!result) {
   PRN(("Mouth not found!\n"));
   return result;
}

// Freeing memory
free(dummy_labels);
FreeImg(im_label);
FreeImg(cb_cropped);
FreeImg(cr_cropped);
FreeImg(mouth_map);

#ifdef FUNC_TRACK
PRN(("Ending FindMouth\n"));
#endif

return true;

}


// This function returns a binary image representing only one eye (which = LEFT, RIGHT)
bool
FindEye(const Image& skin_map, const Image& mouth, Image& eye, int which, const Init_para& parameters, int serialNumber)
{

int x, y, i, n, h_min, h_max, h,
    cx_mouth, max, peak, diff, 
    comp, cy, num;

bool result;

int* ext;
int* proj;
int* labels;
int* centroids;

Image cropped, neg_cropped,
      background, neg_background,
      and_image, search_image, 
      image_label, component;


#ifdef FUNC_TRACK
PRN(("Starting FindEye\n"));
#endif

// Limiting search area: we are cropping our image of the skin map
// to consider only the part above the given 'mouth'. Only the central
// part of this region is actually considered as a valid search area.

BoundingBox(mouth, 0, ext);

if (!CentroidX(mouth, cx_mouth)) {
   PRN(("Eye not found!\n"));
   return false;
}


h_min = (int) (parameters.perc_eye_min_y * ext[2]);
h_max = (int) (parameters.perc_eye_max_y * ext[2]);
h = h_max - h_min + 1;

// The search area is further limited according to which eye we're looking for
if (which == LEFT) {
    ext[0] = 0;
    ext[1] = cx_mouth;
}
else if (which == RIGHT) {
    ext[0] = cx_mouth + 1;
    ext[1] = skin_map.width - 1;
}
ext[2] = h_min;
ext[3] = h_max;


result = CropImage(ext, skin_map, cropped, 0);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}


// Extracting holes in the considered portion of the skin map

result = Negative(cropped, neg_cropped);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}
result = Background(neg_cropped, background);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}
result = Negative(background, neg_background);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}
result = And(neg_background, neg_cropped, and_image);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}

// Clearing the resulting image from single (noisy) pixels
result = EraseSinglePixels(and_image, search_image);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}

// Freeing service images...
FreeImg(cropped);
FreeImg(neg_cropped);
FreeImg(neg_background);
FreeImg(background);
FreeImg(and_image);


// Debug: writing to file the holes image
if (parameters.debug_eyes) {
   char name[MAX_LENGTH_FILE];
   BMP bmp;
   struct stat st;

   sprintf(name, "%s%s%d1.BMP", DATA_FOLDER, parameters.file_debug_eyes, serialNumber);
   if (stat(name, &st) == 0)
      sprintf(name, "%s%s%d2.BMP", DATA_FOLDER, parameters.file_debug_eyes, serialNumber);
   bmp.Write2File(name, search_image);
}

// Labelling Components
result = LabelComponents(search_image, image_label, n, labels);
if (!result) {
   PRN(("Eye not found!\n"));
   return false;
}


// If no components are found in the considered area...
if (n < 1) {
   PRN(("Eye not found!\n"));
   return false;
}

// If only one component is found...
else if (n == 1) {
   comp = labels[0];
   result = GetComponent(image_label, comp, eye);
   return result;
}

// If more than one component are found...
else  {

   // Computing the horizontal projection of our cropped image
   result = BinaryHorProjection(search_image, proj);
   if (!result) {
       PRN(("Eye not found!\n"));
       return false;
   }

  // Looking for the first peak from below (in order to ignore peaks due to the eyebrows)...
     i = h - 2;
     diff = 0;
     peak = -1;
     while (diff >= 0 && i >= 0) {
        diff = proj[i] - proj[i + 1];
        if (diff < 0)
           peak = i + 1;
        i--;
    }


    centroids = (int*) calloc(n, sizeof(int));


    // Centroid Criterion: we keep only those components whose y-centroid is close to the peak
    for (i = 0; i < n; i++) {
          result = GetComponent(image_label, labels[i], component);
          if (!result)
             return result;
          result = CentroidY(component, cy);
          if (!result)
             return result;
          FreeImg(component);
          centroids[i] = cy;
    }

    num = n;

    for (i = 0; i < n; i++) {
         if (labels[i]) {
           if (abs(centroids[i] - peak) > (search_image.height * parameters.eye_offset)) {
                 labels[i]= 0;
                 num--;
            }
         }
    }

    // If no components are left...
    if (num < 1) {
       PRN(("Eye not found!\n"));
       return false;
    }
    
    // If only one component is left...
    else if (num == 1) {
       comp = -1;
       i = 0;
       while (comp < 0 && i < n) {
           if (labels[i])
                comp = labels[i];
           i++;
        }
     }

     // If more than one component are left...
     else {

           // We're keeping the component with lower y-centroid
           max = 0;

           for (i = 0; i < n; i++) {
                      if (labels[i]) {
                            if (centroids[i] > max) {
                                  max = centroids[i];
                                  comp = labels[i];
                            }
                      }
           }
     }


    free(centroids);
    result = GetComponent(image_label, comp, eye);
    if (!result)
        return result;


}

// Freeing memory
free(ext);
free(proj);
free(labels);
FreeImg(search_image);
FreeImg(image_label);


#ifdef FUNC_TRACK
PRN(("Ending FindEye\n"));
#endif

return true;

}


// This function returns a binary image representing an eyebrow
bool
GetEyebrow(const Image& eye, const int* face_bb, const Image& im_exp, Image& eyebrow, const Init_para& parameters, int serialNumber)
{

int x, y, i, dim, min, max, n, 
    area, num, comp, cy, h_max, 
    cx, eye_cx;

bool result;

int* eyeband;
int* labels;
int* centroids;

Image cropped, edges, thres_edges,
      im_label, component;


#ifdef FUNC_TRACK
PRN(("Starting GetEyebrow\n"));
#endif

// Limiting search area: we are cropping our image of the face
// to consider only the part above the given 'eye'
// We add some border on the sides of the bounding box
// in order not to lose details of the eyebrow
result = BoundingBox(eye, 0, eyeband);
if (!result)
   return false;

eyeband[3] = eyeband[2];
eyeband[2] = eyeband[3] + (face_bb[2] - eyeband[3])/2;

result = CropImage(eyeband, im_exp, cropped, LARGE_BORDER);
if (!result)
   return false;

// Searching for horizontal borders...
result = SobelY(cropped, edges);
if (!result)
   return false;
result = ThresholdPerc(edges, thres_edges, parameters.edges_eyebr_perc, parameters.eyebr_th_max, parameters.eyebr_th_min);
if (!result)
   return false;

// Freeing memory
FreeImg(cropped);
FreeImg(edges);


// Debug: writing to file the thresholded image
if (parameters.debug_eyebr) {
   char name[MAX_LENGTH_FILE];
   BMP bmp;
   struct stat st;

   sprintf(name, "%s%s%d1.BMP", DATA_FOLDER, parameters.file_debug_eyebr, serialNumber);
   if (stat(name, &st) == 0)
      sprintf(name, "%s%s%d2.BMP", DATA_FOLDER, parameters.file_debug_eyebr, serialNumber);
   bmp.Write2File(name, thres_edges);
}



// Setting border to 0...
for (y = thres_edges.height - 1 ; y >= 0; y--) {
    thres_edges.data[0][y] = thres_edges.data[thres_edges.width - 1][y] = BLACK;
}
for (x = 1; x < thres_edges.width - 1; x++) {
    thres_edges.data[x][0] = thres_edges.data[x][thres_edges.height - 1] = BLACK;
}

// Looking for connected components
result = LabelComponents(thres_edges, im_label, n, labels);
if (!result) {
   PRN(("NO EYEBROWS FOUND!\n"));
   return false;
}


dim = thres_edges.width * thres_edges.height;
min = (int) (parameters.perc_eyebr_min * dim);
max = (int) (parameters.perc_eyebr_max * dim);


// If there are no components in the eyebrow area...
if (n == 0) {
   PRN(("NO EYEBROWS FOUND!\n"));
   return false;
}

num = n;


// Area Criterion: a component can be an eyebrow only if its area is reasonable.
// Components that don't satisfy this criterion are dismissed.
for (i = 0; i < n; i++) {
     result = GetComponent(im_label, labels[i], component);
     if (!result)
        return false;
     area = Area(component);
     if (! (area >= min && area <= max)) {
          labels[i] = 0;
          num--;
     }
     FreeImg(component);
}

// If no components are left...
if (num == 0) {
   PRN(("NO EYEBROWS FOUND!\n"));
   return false;
}

// If only one component is left...
else if (num == 1) {
      i = 0;
      while(!labels[i]) {
             i++;
      }
      comp = i;
      
      result = GetComponent(im_label, labels[comp], eyebrow);
      if (!result)
         return false;

      #ifdef FUNC_TRACK
      PRN(("Ending GetEyebrow\n"));
      #endif
      
      return true;
}

// If more than one component are left...
else if (num > 1) {

      centroids = (int*) calloc(n, sizeof(int));


      // Centroid Criterion: we're skipping the first 'offset' pixels from the bottom border (in order to avoid che upper edge of the eye)

      h_max = im_label.height - (int)(im_label.height * parameters.eyebr_offset);

      for (i = 0; i < n; i++) {
           if (labels[i]) {
              result = GetComponent(im_label, labels[i], component);
              if (!result)
                 return false;
              result = CentroidY(component, cy);
              if (!result)
                 return false;

              centroids[i] = cy;
              if (! (centroids[i] <= h_max)) {
                    labels[i] = 0;
                    num--;
              }
              FreeImg(component);
           }
      }
      
      free(centroids);

      // If no components are left...
      if (num == 0) {
         PRN(("NO EYEBROWS FOUND!\n"));
         return false;
      }

      // If only one component is left...
      else if (num == 1) {
           i = 0;
           while(!labels[i]) {
                i++;
           }
           comp = i;

           result = GetComponent(im_label, labels[comp], eyebrow);
           if (!result)
              return false;

           #ifdef FUNC_TRACK
           PRN(("Ending GetEyebrow\n"));
           #endif

           return true;
      }

      // If more than one component are left...
      else  {

            centroids = (int*) calloc(n, sizeof(int));
            result = CentroidX(eye, eye_cx);
            if (!result)
               return false;

            // Centroid Criterion: eyebrow x-centroid must be similar to eye x-centroid
            // We choose as eyebrow the component whose x-centroid is closer to the eye's.
            for (i = 0; i < n; i++) {
                if (labels[i]) {
                   result = GetComponent(im_label, labels[i], component);
                   if (!result)
                      return false;
                   result = CentroidX(component, cx);
                   if (!result)
                      return false;
                   centroids[i] = cx + im_label.offsetX;
                   FreeImg(component);
                }
            }


            min = im_label.width;

            for (i = 0; i < n; i++) {
                if (labels[i]) {
                      if (abs(centroids[i] - eye_cx) < min) {
                                     min = abs(centroids[i] - eye_cx);
                                     comp = i;
                      }

                }

            }
      

            free(centroids);
            result = GetComponent(im_label, labels[comp], eyebrow);
            if (!result)
               return false;
        }

}

// Freeing memory...
free(eyeband);
free(labels);
FreeImg(thres_edges);
FreeImg(im_label);


#ifdef FUNC_TRACK
PRN(("Ending GetEyebrow\n"));
#endif

return true;

}

// This function simultaneously returns both eyebrows images.
bool
GetBothEyebrows(const Image& eye_left, const Image& eye_right, const int* face_bb, const Image& im_exp, 
                Image& eyebrow_left, Image& eyebrow_right, const Init_para& parameters, int serialNumber)
{

bool result;

#ifdef FUNC_TRACK
PRN(("Starting GetBothEyebrows\n"));
#endif

result  = GetEyebrow(eye_left, face_bb, im_exp, eyebrow_left, parameters, serialNumber);
if (!result)
   return false;
result = GetEyebrow(eye_right, face_bb, im_exp, eyebrow_right, parameters, serialNumber);
if (!result)
   return false;


#ifdef FUNC_TRACK
PRN(("Ending GetBothEyebrows\n"));
#endif

return true;

}




// This function finds the mouth/eye corners (in the expanded image) and return their coordinates
bool
FindCorners(const Image& image, int*& ext, const Init_para& parameters, int serialNumber)
{

int x, y, min, max;
bool result;
Image mean, sobel, th_sobel;

#ifdef FUNC_TRACK
PRN(("Starting FindCorners\n"));
#endif


ext = (int*) calloc(4, sizeof(int));
if (ext == NULL)
   return false;

// The image is filtered to show horizontal borders and then thresholded.
// Then we look for the two extreme white pixels.
result = Mean(image, mean);
if (!result)
  return false;
result = SobelY(mean, sobel);
if (!result)
  return false;
result = ThresholdPerc(sobel, th_sobel, parameters.extr_mouth_perc, parameters.extr_mouth_th_min, parameters.extr_mouth_th_max);
if (!result)
  return false;


min = sobel.width;
max = 0;

// not considering image borders...
for (x = 1; x < th_sobel.width - 1; x++) {
     for (y = 1; y < th_sobel.height - 1; y++) {
             if (th_sobel.data[x][y]) {
                if (x < min) {
                   min = x;
                   ext[0] = min;
                   ext[1] = y;
                }
                if (x > max) {
                   max = x;
                   ext[2] = max;
                   ext[3] = y;
                }
             }
       }
}

// debug
if (parameters.debug_mouth) {
   char name[MAX_LENGTH_FILE];
   BMP bmp;

   sprintf(name, "%s%s%d.BMP", DATA_FOLDER, parameters.file_debug_mouth, serialNumber);
   bmp.Write2File(name, th_sobel);

   sprintf(name, "%s%scorners%d.BMP", DATA_FOLDER, parameters.file_debug_mouth, serialNumber);
   ShowPoints(name, image, ext);
}

// Freeing memory...
FreeImg(mean);
FreeImg(sobel);
FreeImg(th_sobel);

#ifdef FUNC_TRACK
PRN(("Ending FindCorners\n"));
#endif

return true;

}



// This function crops the input image according to the bounding box 'ext'.
// If 'border' is not 0, the width of the resulting image is increased
// by 'border' pixels on the left and on the right of the actual bounding box.
bool
CropImage(int* ext, const Image& image, Image& cropped, int border)
{

int i, x, y, x1, y1, x_min, x_max;

#ifdef FUNC_TRACK
PRN(("Starting CropImage\n"));
#endif

cropped.width  = ext[1] - ext[0] + 1;
cropped.height = ext[3] - ext[2] + 1;

x_min = ext[0];
x_max = ext[1];


if (border) {
   if (x_min - border >= 0) {
      x_min         -= border;
      cropped.width += border;
      }
   if (x_max + border < image.width) {
      x_max         += border;
      cropped.width += border;
      }
}

// The offsetX and offsetY fields of the cropped image are set
// so that we can map coordinates in the cropped image into the relative 
// coordinates in the original image
cropped.offsetX = image.offsetX + x_min;
cropped.offsetY = image.offsetY + ext[2];

cropped.data = (byte**) calloc(cropped.width, sizeof(byte*));
if (cropped.data == NULL)
   return false;
for (i = 0; i < cropped.width; i++) {
    cropped.data[i] = (byte*) calloc(cropped.height, sizeof(byte));
    if (cropped.data[i] == NULL)
       return false;
}


for (x = 0, x1 = x_min; x < cropped.width, x1 <= x_max; x++, x1++) {
     for (y = 0, y1 = ext[2]; y < cropped.height, y1 <= ext[3]; y++, y1++) {
                 cropped.data[x][y] = image.data[x1][y1];

        }
}


#ifdef FUNC_TRACK
PRN(("Ending CropImage\n"));
#endif

return true;

}

// Copying an image into a new one
bool
CopyImage(const Image& from, Image& to)
{

int x, y;

#ifdef FUNC_TRACK
PRN(("Starting CopyImage\n"));
#endif

to = Image(from.width, from.height, from.offsetX, from.offsetY);
if (to.data == NULL)
   return false;

for (x = 0; x < to.width; x++) {
       for (y = 0; y < to.height; y++) {
              to.data[x][y] = from.data[x][y];

       }
}

#ifdef FUNC_TRACK
PRN(("Ending CopyImage\n"));
#endif

return true;

}

// Freeing memory used by an Image
void
FreeImg(Image& image)
{

int i;

#ifdef FUNC_TRACK
PRN(("Starting FreeImg\n"));
#endif

if (image.data != NULL) {
    for (i = 0; i < image.width; i++)
        free(image.data[i]);
    free(image.data);
}
image.data = NULL;

#ifdef FUNC_TRACK
PRN(("Ending FreeImg\n"));
#endif

return;

}


// This function shows a pair of points (i.e., mouth corners), highlighting them on the input image
bool
ShowPoints(char* path, const Image& image, int* points)
{

bool result;
Image to_print;
BMP bmp;

#ifdef FUNC_TRACK
PRN(("Starting ShowPoints\n"));
#endif

to_print = Image(image.width, image.height, image.offsetX, image.offsetY);
if (to_print.data == NULL)
   return false;

to_print.data[points[0]][points[1]] =  to_print.data[points[2]][points[3]] = WHITE;

result = bmp.WriteAndHighlight(path, image, to_print);
FreeImg(to_print);

#ifdef FUNC_TRACK
PRN(("Ending ShowPoints\n"));
#endif

return result;

}

// This function shows the bounding box of a region of interest, highlighting them on the input image
bool
ShowBoundingBox(char* path, const Image& image, int* points) {

int i;
bool result;
Image to_print;
BMP bmp;

#ifdef FUNC_TRACK
PRN(("Starting ShowBoundingBox\n"));
#endif

to_print = Image(image.width, image.height, image.offsetX, image.offsetY);
if (to_print.data == NULL)
   return false;

for (i = points[0]; i <= points[1]; i++)
    to_print.data[i][points[2]] = to_print.data[i][points[3]] = WHITE;
for (i = points[2] + 1; i < points[3]; i++)
    to_print.data[points[0]][i] = to_print.data[points[1]][i] = WHITE;

result = bmp.WriteAndHighlight(path, image, to_print);
FreeImg(to_print);

#ifdef FUNC_TRACK
PRN(("Ending ShowBoundingBox\n"));
#endif

return result;

}


// This function shows a given component highlighting it on a given image
bool
ShowComp(char* path, const Image& image, const Image& comp)
{

int x, y, ox, oy;
bool result;
BMP bmp;
Image to_print;

#ifdef FUNC_TRACK
PRN(("Starting ShowComp\n"));
#endif

to_print = Image(image.width, image.height, image.offsetX, image.offsetY);
if (to_print.data == NULL)
   return false;

ox = comp.offsetX;
oy = comp.offsetY;

for (x = 0; x < comp.width; x++) {
      for (y = 0; y < comp.height; y++) {
                if (comp.data[x][y])
                   to_print.data[x + ox][y + oy] = WHITE;

      }
}

result = bmp.WriteAndHighlight(path, image, to_print);
FreeImg(to_print);

#ifdef FUNC_TRACK
PRN(("Ending ShowComp\n"));
#endif

return result;

}



