//==============================================
//  copyright            : (C) 2003-2005 by Will Stokes
//==============================================
//  This program is free software; you can redistribute it 
//  and/or modify it under the terms of the GNU General 
//  Public License as published by the Free Software 
//  Foundation; either version 2 of the License, or  
//  (at your option) any later version.         
//==============================================

//Systemwide includes
#include <qpainter.h>
#include <qimage.h>
#include <math.h>

//Projectwide includes
#include "edgeDetect.h"
#include "blur.h"
#include "../enhancements/contrast.h"

//----------------------------------------------
// Inputs:
// -------
// QImage* image - image to find edges in
//
// Outputs:
// --------
// QImage* getEdgeImage - returns the produced grayscale edge image
// 
// Other information such as pixel groups etc is also availble through 
// various accesor method.
//
// Description:
// ------------
// This class is the first known publically available implementation of 
// Kim, Lee, and Kweon's 2003 paper titled:
// "Automatic edge detection using 3x3 ideal binary 
// pixel patterns and fuzzy-based edge thresholding"
// http://rcv.kaist.ac.kr/publication/file/foreign_journal/28_DongSuKim_PRL2004.pdf
//
// Edge detection is an old problem, an while many use the edge detectors by 
// Canny, Sobel, etc, they all suffer from a common problem: the user must 
// tweak a series of poorly understood input parameters to get the ideal edge image.
//
// Album Shaper was in need of an automatic edge detector for use when 
// blurring and sharpening images. Having the user manually tweak such 
// paramters first would not only be annoying but error-prone. In an 
// effort to make edge detection automatic I took a stab at implementing 
// this paper and am quite happy with the resulsts...
//
// http://albumshaper.sourceforge.net/images/teasers/peaksAndValleys.jpg
//
// Algorithm:
// ----------
// While complex, the algorithm can be broken up into a series of 
// fairly straightforward tasks:
//
// 1.) "allocateAndInitObjects()" is called to allocate and fill a 
//     few data structures that will be used when finding image edges.
//
// 2.) "fillLumMapAndLumHistogram"() is called, during which an 
//     luminance map is constructed and luminance histogram populated. 
//     For an m x n image, the luminance map will be a m x n integer array.
//
// 3.) The luminance histogram is smoothed using "smoothLumHistogram" to 
//     make peak finding less sensative to noise.
//
// 4.) The fourth step is a little complicated. The edge magnitude and 
//     GSLC (Grey level similitude code) value is computed at each pixel.
//     The paper takes an interesting approach to edge detection by 
//     classifying pixels into one of 9 groups by first computing the
//     average luminance for a 3x3 group centered about a pixel. Pixels are
//     then separated into one of two groups, those that have a luminance 
//     greater than or less than the 3x3 average luminance. For example:
//                                   
//                                                     X
//      ---------          ---------        ---------X
//     | 7 15 18 |        | 0  1  1 |      | 0  1  X |   
//     | 5 17 20 |  -->   | 0  1  1 |  --> | 0  X  1 |
//     | 9  8  3 |        | 0  0  0 |      | X  0  0 |
//      ---------          ---------       X---------
//                                       X
//
//     Here the average luminance is 11.333, placing the top right 4 
//     pixels in one group and the other remaining pixels in another. 
//     The dominant edge diretion is from the bottom left the top right. 
//     The GLSC code is computing by considering the 1/0 values
//     (1 = pixel in same group as central pixel). The central value 
//     is ignored (it's always 1) leaving us with an 8bit = 2^8 = 256 code.
//
//     In this case the GSLC is:
//     0*2^0 + 1*2^1 + 1*2^2 +
//     0*2^3 + XXXXX + 1*2^4 +
//     0*2^5 + 0*2^6 + 0*2^7 = 22.
//
//     The authors of the paper found pixels with one of five GSLC 
//     codes (15,31,7,47,11) were most often associated with edge pixels 
//     when producing an edge image using various competitive techinques. 
//     By looking up the GSLC for a given pixel later one we can suppress 
//     edges where they most likely do not belong.
//
// 5.) The fifth step involves grouping pixels by luminance using the 
//     smoothed luminance histogram. This complicated step is brushed off 
//     as being trivial in the paper. Since I struggled a bit with developing 
//     an algorithm for this step, I'll explain my approach in detail 
//     to avoid others suffering.
//
//     Using the smoothed histogram, we first compute the JND or just 
//     noticible differnce using the maximum count that was found. I'm 
//     not sure how appropriate this is, but 2% is the usual quantity
//     used in other contexts, and it works well here, so 2% it is.
// 
//     Next we walk through the smoothed luminance histogram and find 
//     the midpoint of the valleys. We accomplish this by updating an 
//     index of the deepend last known valley. As the valley slopes
//     down and we move across it we update this last best known index. 
//     Once the valley slopes up one JND above the last deepend location 
//     found we mark that valley midpoint and move on.
//
//     Once all valley midpoints have been marked we can quickly deduce 
//     how many peaks must exist. We pass across the smoothed luminance 
//     histogram again finding the peak index for each pixel between
//     valley midpoints. For all future work pixels one JND+- the peak 
//     center will be used.
//
// 6.) The sixth step, "computeClusterStatistics()", computes various 
//     cluster-specific statsitics that will be used to determine cluster 
//     thresholds. The paper was rather vague in this area, but after 
//     experimenting with various interpreatations of what they were trying 
//     to say I think I finally got it right.
// 
//     First, we iterate over all image pixels, determine which pixel group 
//     they belong by comparing luminance endpoints for all pixelclusters, 
//     and update total edge mag, num pixels, and an edgeMagHistogram for 
//     the given pixel group they belong to.
//
//     Next we compute the average edge meganitude and most frequent edge 
//     magnitude observed for each pixel cluster, in addition to normalizing 
//     the cluster pixel count variable to [0,1]
//
// 7.) The seventh step is quite complicated and encompases computing the 
//     edge thresholds for each pixel cluster using the 18-rule fuzzy 
//     logic approach put forward by the paper. There is nothing ground 
//     breaking here, just a lot of complicated fuzzy logic, although most 
//     of the effort is put into computing the centroid at the end. I had 
//     never touched fuzzy logic before, but found this article more than
//     helpful getting myself up to speed:
//
//     http://www.doc.ic.ac.uk/~nd/surprise_96/journal/vol2/sbaa/article2.html
//
// 8.) The eigth and final step, actually constructing the edge image, is fairly
//     straight forward and employs techinques (such as non-maximum suppression or NMS) 
//     anyone who has implemented Canny shoudl be familiar with.
//
//     First, one looks up the ESF (edge shape factor) for a given pixel from a look-up 
//     table. These values were computed by generating a ton of edge images by carefully
//     setting Canny, Sobel, etc params, then identifying how often a pixel
//     with a given GSLC code was judged to be an edge. Hence ESF's will fall in
//     the [0,1] range and help to suppress edges where no clear edge really
//     can be claimed to exist, e.g.
//
//     . # .
//     # # #   <- If there is an edge here, care to explain what the
//     . # .      edge direction is?!!
//
//
//     If the ESF for a give pixel is 0 we skip it, it's not an edge.
//
//     Next, we look up a pixels edge magnitude threshold by identifying the
//     pixel cluster is belongs to using the pixels luminance using the
//     luminance map.
// 
//     If the pixels edge magnitude is less than the threshold we skip the pixel, 
//     this filtersout low lying noise.
//
//     Finally, the direction of the pixel is looked up using its GSLC and
//     NMS (non-maximum suppression is applied). If the pixel has a greater
//     egdge magnitude than either of its neighbors along the edge direction then
//     the edge is marked.
//
// Final Remarks:
// --------------
// Despite the involved complexity, the implementation appears to work
// really really well. I consider this one of the secret gems of Album Shaper
// and hope to make good use of it in the future for things other than 
// just sharpening and bluring. The only caveat is that edge detection 
// does take a few CPU cycles.
//----------------------------------------------


//==============================================
EdgeDetect::EdgeDetect( QImage* image )
{                     
  //load image
  this->image = image;

  //allocate and initialize objects used for edge detection
  allocateAndInitObjects();

  //fill lum map and lum histogram
  fillLumMapAndLumHistogram();
  
  //fill smoothed histogram
  smoothLumHistogram();
  
  //compute edge magnitude and GSLC maps
  computeEdgeMagAndGSLCmaps();
  
  //determine pixel clusters
  findPixelClusters();
  
  computeClusterStatistics();  
 
  computeClusterThresholds();
  
  constructEdgeImage();
}
//==============================================
EdgeDetect::~EdgeDetect()
{
  deallocateObjects();
}
//==============================================
int EdgeDetect::getNumClusters()
{ return numClusters; }
//==============================================
PixelCluster* EdgeDetect::getClusters()
{ return clusters; }
//==============================================
int* EdgeDetect::getPeaks()
{ return clusterPeaks; }
//==============================================
int* EdgeDetect::getSmoothHist()
{ return smoothLumHist; }
//==============================================
QImage* EdgeDetect::getEdgeImage()
{
  return image; 
}
//==============================================
int* EdgeDetect::getClusterMap()
{
  //construct map
  int* clusterMap = new int[image->width() * image->height()];
  
  //iterate over all pixels, determine cluster each pixel belongs to
  int i, cluster;
  for(i=0; i<image->width()*image->height(); i++)
  {
    for(cluster=0; cluster<numClusters; cluster++)
    {
      if( lumMap[i] >= clusters[cluster].minLuminance &&
          lumMap[i] <= clusters[cluster].maxLuminance )
      {
        clusterMap[i] = cluster;
        break;
      }
    } //cluster
  } //pixel

  return clusterMap;
}
//==============================================
void EdgeDetect::allocateAndInitObjects()
{
  //initialize: 
  //-luminosity histogram
  //-smoothed luminosity histogram
  //-identified peak regions
  int i;
  for(i=0; i<256; i++)
  { 
    lumHist[i] = 0; 
    smoothLumHist[i] = 0;
    clusterPeaks[i] = 0;
  }
  
  //allocate luminance map
  lumMap = new int[image->width() * image->height()];
  
  //allocate edge magnitude map
  edgeMagMap = new float[image->width() * image->height()];
  
  //allocate GSLC map
  GSLCmap = new int[image->width() * image->height()];
  
  //construct LUT
  constructGSLClut();
}
//==============================================
void EdgeDetect::fillLumMapAndLumHistogram()
{
  int x, y;
  QRgb* rgb;
  uchar* scanLine;
  int lumVal;
  for( y=0; y<image->height(); y++)
  {   
    scanLine = image->scanLine(y);
    for( x=0; x<image->width(); x++)
    {
      //get lum value for this pixel
      rgb = ((QRgb*)scanLine+x);
      lumVal = qGray(*rgb);
      
      //store in lum map
      lumMap[x + y*image->width()] = lumVal;
      
      //update lum histogram
      lumHist[ lumVal ]++;
    }
  }
}
//==============================================
void EdgeDetect::smoothLumHistogram()
{
  #define FILTER_SIZE 5
  int filter[FILTER_SIZE] = {2, 5, 8, 5, 2};
  
  int i,j;
  int filterIndex, sum, total;
  for(i = 0; i<256; i++)
  {
    sum = 0;
    total = 0;
    
    for( j= -FILTER_SIZE/2; j <= FILTER_SIZE/2; j++)
    {
      if( i+j > 0 && i+j < 256 )
      {
        filterIndex = j+ FILTER_SIZE/2;
        total+= filter[filterIndex] * lumHist[i+j];
        sum  += filter[filterIndex];
      }
    }
    
    smoothLumHist[i] = total / sum;
  }
}
//==============================================
void EdgeDetect::computeEdgeMagAndGSLCmaps()
{
  int x, y;
  int idealPattern[9];
  int pixelLums[9];
  
  //-------  
  //iterate over all pixels
  for( y=0; y<image->height(); y++)
  {   
    for( x=0; x<image->width(); x++)
    {
      //compute pixel luminances for entire grid
      pixelLums[0] = pixelLum(x-1,y-1);
      pixelLums[1] = pixelLum(x  ,y-1);
      pixelLums[2] = pixelLum(x+1,y-1);
      pixelLums[3] = pixelLum(x-1,y  );
      pixelLums[4] = pixelLum(x  ,y  );
      pixelLums[5] = pixelLum(x+1,y  );
      pixelLums[6] = pixelLum(x-1,y+1);
      pixelLums[7] = pixelLum(x  ,y+1);
      pixelLums[8] = pixelLum(x+1,y+1);
      
      //compute average
      float avg = 0;
      int i;
      for(i=0; i<=8; i++)
      {
        avg+= pixelLums[i];
      }
      avg = avg / 9;
      
      //determine ideal pattern and I0 and I1 averages
      int centerPixelLum = pixelLums[4];
      float centerDiff = centerPixelLum - avg;
      
      float I0avg = 0;
      int I0count = 0;
      
      float I1avg = 0;
      int I1count = 0;
      
      for(i=0; i<=8; i++)
      {
        if( centerDiff * (pixelLums[i]-avg) >=0 )
        { 
          I1avg+=pixelLums[i];
          I1count++;
          idealPattern[i] = 1; 
        }
        else 
        { 
          I0avg+=pixelLums[i];
          I0count++;
          idealPattern[i] = 0; 
        }
      }

      //compute and store edge magnitude
      if(I0count > 0) I0avg = I0avg/I0count;
      if(I1count > 0) I1avg = I1avg/I1count;     
      edgeMagMap[x + y*image->width()] = QABS( I1avg - I0avg );
      
      //compute and store GSLC
      int GSLC=0;
      int weight = 1;
      for(i=0; i<9; i++)
      {
        //skip center
        if(i == 4) continue;
        
        if(idealPattern[i] == 1)
        { GSLC+=weight; }
        
        weight = weight*2;
      }
      GSLCmap[x + y*image->width()] = GSLC;
    } //x
  } //y
}
//==============================================
int EdgeDetect::pixelLum(int x, int y)
{
  int clampedX = QMAX( QMIN( x, image->width()-1), 0);
  int clampedY = QMAX( QMIN( y, image->height()-1), 0);
  return lumMap[ clampedX + clampedY * image->width() ];
}
//==============================================
void EdgeDetect::findPixelClusters()
{
  //find max count
  int maxCount = 0;
  int i;
  for(i=0; i<256; i++)
  {
    if(smoothLumHist[i] > maxCount)
      maxCount = smoothLumHist[i];
  }

  //compute JND for histogram (2% of total spread)
  int histJND = maxCount/50;

  //construct temporary array for valley locations
  //1's will indicate a valley midpoint
  int tmpValleyArray[256];
  for(i=0; i<256; i++) { tmpValleyArray[i] = 0; }
  
  //move across histogram finding valley midpoints
  int curTrackedMin = smoothLumHist[0];
  
  //first and last indices tracked min was observed
  int firstMinIndex = 0;
  int lastMinIndex = 0;
  
  //only add valley midpoint if finished tracking a descent
  bool slopeNeg = false;
  
  for(i = 1; i<256; i++ )
  {
    if( smoothLumHist[i] < curTrackedMin - histJND )
    {
      //found a descent!
      slopeNeg = true;
      curTrackedMin = smoothLumHist[i];
      firstMinIndex = i;
    }
    //starting to go up again, add last min to list
    else if( smoothLumHist[i] > curTrackedMin + histJND )
    {
      //if finished tracing a negative slope find midpoint and set location to true
      if(slopeNeg)
      {
        tmpValleyArray[ (firstMinIndex + lastMinIndex)/2 ] = 1;
      }
      
      curTrackedMin = smoothLumHist[i];
      slopeNeg = false;
    }
    else
    {
      //still tracking a min, update the right 
      //hand index. center of valley is found
      //by averaging first and last min index
      lastMinIndex = i;
    }
  }
  
  //count valleys
  int numValleys = 0;
  for(i=0; i<256; i++)
  {
    if(tmpValleyArray[i] == 1 ) numValleys++;
  }

  //determine number of clusters
  numClusters = numValleys-1;
  if(tmpValleyArray[0] != 1)
    numClusters++;
  if(tmpValleyArray[255] != 1)
    numClusters++;
  
  //allocate clusters
  clusters = new PixelCluster[numClusters];
  
  //automatically start first cluster
  int cluster=0;
  clusters[cluster].minLuminance = 0;
  
  //initialize left and right boundaries of all clusters
  for(i=1; i<256; i++)
  {
    //reached next valley, end cluster
    if( tmpValleyArray[i] == 1)
    {
      clusters[cluster].maxLuminance = i-1;
      cluster++;
      clusters[cluster].minLuminance = i;
    }
    //end last cluster automatically at end
    else if(i == 255)
    {
      clusters[cluster].maxLuminance = i;
    }
  }
  
  //determine cluster peaks
  for(cluster=0; cluster<numClusters; cluster++)
  {
    //find max for current cluster
    int maxIndex = clusters[cluster].minLuminance;
    for(i=clusters[cluster].minLuminance; i<=clusters[cluster].maxLuminance; i++)
    {
      if(smoothLumHist[i] > smoothLumHist[maxIndex])
         maxIndex = i;
    }
    
    //mark peaks  
    int lumJND = 255/50;
    for(i=QMAX(0, maxIndex-lumJND); i<QMIN(256, maxIndex+lumJND); i++)
    { 
      clusterPeaks[i] = 1; 
    }
  }
}
//==============================================
void EdgeDetect::computeClusterStatistics()
{
  //initialize cluster stats
  int cluster;
  for(cluster=0; cluster<numClusters; cluster++)
  {
    int i;
    for(i=0; i<256; i++)
    {
      clusters[cluster].edgeMagHistogram[i] = 0;
    }
    clusters[cluster].totalEdgeMagnitude=0.0f;
    clusters[cluster].numPixels = 0;
  }
  
  //iterate over all pixels
  int i;
  for(i=0; i<image->width()*image->height(); i++)
  {
    //skip pixels that don't belong to peaks
    if( clusterPeaks[ lumMap[i] ] != 1)
      continue;
    
    //determine cluster pixel belongs to
    int cluster;
    for(cluster=0; cluster<numClusters; cluster++)
    {
      if( lumMap[i] >= clusters[cluster].minLuminance &&
          lumMap[i] <= clusters[cluster].maxLuminance )
      {      
        clusters[cluster].totalEdgeMagnitude+= edgeMagMap[i]; 
        clusters[cluster].numPixels++;
        clusters[cluster].edgeMagHistogram[ QMIN( QMAX( (int)edgeMagMap[i], 0), 255) ]++;
        break;
      }
    } //cluster
  } //pixel i
  
  //iterate over clusters to determine min and max peak cluster sizes
  minClusterSize = clusters[0].numPixels;
  maxClusterSize = clusters[0].numPixels;
  for(cluster=1; cluster<numClusters; cluster++)
  {
    if(clusters[cluster].numPixels < minClusterSize)
      minClusterSize = clusters[cluster].numPixels;

    if(clusters[cluster].numPixels > maxClusterSize)
      maxClusterSize = clusters[cluster].numPixels;
  }

  //iterate over clusters one final time to deduce normalized inputs to fuzzy logic process
  int JND = 255/50;
  for(cluster=0; cluster<numClusters; cluster++)
  {
    clusters[cluster].meanMode = QMIN( clusters[cluster].totalEdgeMagnitude / clusters[cluster].numPixels,
                                       3*JND );
    
    int i;
    int mode = 0;
    for(i=1; i<256; i++)
    {
      if( clusters[cluster].edgeMagHistogram[i] > clusters[cluster].edgeMagHistogram[ mode ] )
        mode = i;
    }
    clusters[cluster].mode = QMIN( mode, 2*JND );
        
    clusters[cluster].pixelCount = ((float)(clusters[cluster].numPixels - minClusterSize)) / 
                                   (maxClusterSize - minClusterSize);  
  }
}
//==============================================
//compute edge thresholds for each cluster using 18-rule fuzzy logic approach
void EdgeDetect::computeClusterThresholds()
{
  //iterate over each cluster
  int cluster;
  float S1,M1,L1;
  float S2,M2,L2;
  float S3,L3;
  float outS, outM, outL;
  
  int JND = 255/50;
  
  for(cluster=0; cluster<numClusters; cluster++)
  {
    //----
    //compute S,M, and L values for each input
    //----
    S1 = QMAX( 1.0f - ((clusters[cluster].meanMode/JND) / 1.5f), 0 );

    if( (clusters[cluster].meanMode/JND) <= 1.5f )
      M1 = QMAX( (clusters[cluster].meanMode/JND) - 0.5f, 0 );
    else
      M1 = QMAX( 2.5f - (clusters[cluster].meanMode/JND), 0 );
    
    L1 = QMAX( ((clusters[cluster].meanMode/JND) - 1.5f) / 1.5f, 0 );
    //----
    S2 = QMAX( 1.0f - (clusters[cluster].mode/JND), 0 );
    
    if( (clusters[cluster].mode/JND) <= 1.0f )
      M2 = QMAX( -1.0f + 2*(clusters[cluster].mode/JND), 0 );
    else
      M2 = QMAX( 3.0f - 2*(clusters[cluster].mode/JND), 0 );
    
    L2 = QMAX( (clusters[cluster].mode/JND) - 1.0, 0 );
    //----
    S3 = QMAX( 1.0f - 2*clusters[cluster].pixelCount, 0 );
    L3 = QMAX( -1.0f + 2*clusters[cluster].pixelCount, 0 );
    //----
    
    //Compute M,L for outputs using set of 18 rules.
    //outS is inherantly S given the ruleset provided
    outS = 0.0f;
    outM = 0.0f;
    outL = 0.0f;
    //Out 1
    if( numClusters > 2 )
    {
      outM += S1*S2*S3;   //rule 1
      
      //rule 2
      if( clusters[cluster].meanMode < clusters[cluster].mode )
        outS += S1*S2*L3;   
      else
        outM += S1*S2*L3;

      outM += S1*M2*S3;   //rule 3
      outM += S1*M2*L3;   //rule 4
      outM += S1*L2*S3;   //rule 5
      outM += S1*L2*L3;   //rule 6
      outM += M1*S2*S3;   //rule 7
      outM += M1*S2*L3;   //rule 8
      outM += M1*M2*S3;   //rule 9
      outL += M1*M2*L3;   //rule 10
      outM += M1*L2*S3;   //rule 11
      outL += M1*L2*L3;   //rule 12
      outM += L1*S2*S3;   //rule 13
      outL += L1*S2*L3;   //rule 14
      outM += L1*M2*S3;   //rule 15
      outL += L1*M2*L3;   //rule 16
      outL += L1*L2*S3;   //rule 17
      outL += L1*L2*L3;   //rule 18
    }
    //Out 2
    else
    {
      outL += S1*S2*S3;   //rule 1
      outL += S1*S2*L3;   //rule 2
      outM += S1*M2*S3;   //rule 3
      outL += S1*M2*L3;   //rule 4
      outM += S1*L2*S3;   //rule 5
      outM += S1*L2*L3;   //rule 6
      outL += M1*S2*S3;   //rule 7
      outL += M1*S2*L3;   //rule 8
      outL += M1*M2*S3;   //rule 9
      outL += M1*M2*L3;   //rule 10
      outL += M1*L2*S3;   //rule 11
      outL += M1*L2*L3;   //rule 12
      outL += L1*S2*S3;   //rule 13
      outL += L1*S2*L3;   //rule 14
      outL += L1*M2*S3;   //rule 15
      outL += L1*M2*L3;   //rule 16
      outL += L1*L2*S3;   //rule 17
      outL += L1*L2*L3;   //rule 18
    }

    //find centroid - Beta[k]
    float A = outM + 0.5f;
    float B = 2.5f - outM;
    float C = 1.5f * (outL + 1);
    float D = 1.5f * (outM + 1);
    float E = 2.5f - outL;

    //---------------------------------------------------------------
    //Case 1: Both outM and outL are above intersection point of diagonals
    if( outM > 0.5f && outL > 0.5f )
    {
      //find area of 7 subregions
      float area1 = ((A-0.5f)*outM)/2;
      float area2 = outM * (B-A);
      float area3 = ((2.1f-B) * (outM - 0.5)) / 2;
      float area4 = (2.1 - B) * 0.5f;
      float area5 = ((C - 2.1f) * (outL - 0.5)) / 2;
      float area6 = (C - 2.1f) * 0.5f;
      float area7 = (3.0f - C) * outL;
     
      //find half of total area
      float halfArea = (area1 + area2 + area3 + area4 + area5 + area6 + area7) / 2;
      
      //determine which region split will be within and resulting horizontal midpoint
      
      //Within area 1
      if( area1 > halfArea )
      {
        clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);       
      }
      //Within area 2
      else if( area1 + area2 > halfArea )
      {
        clusters[cluster].beta = ((halfArea - area1) / outM) + A;        
      }
      //Within area 3-4
      else if( area1 + area2 + area3 + area4 > halfArea )
      {
        float a = -0.5f;
        float b = 2.8f;
        float c = area1 + area2 + area3 - halfArea - B/2 - 2.625f;
        clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
      }
      //Within area 5-6
      else if( area1 + area2 + area3 + area4 + area5 + area6 > halfArea )
      {
        float a = 1.0f/3;
        float b = -0.7f;
        float c = area1 + area2 + area3 + area4 - halfArea;
        clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
      }
      //Within area 7
      else
      {
        clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4 + area5 + area6) ) / outL) + C;       
      }
    } //end case 1
    //---------------------------------------------------------------
    //Case 2
    else if ( outM < 0.5f && outL > outM )
    {
      //find area of 5 subregions
      float area1 = (outM*(A-0.5f)) / 2;
      float area2 = (D-A) * outM;
      float area3 = ((C-D) * (outL - outM)) / 2;
      float area4 = (C-D) * outM;
      float area5 = (3.0f - C) * outL;
      
      //find half of total area
      float halfArea = (area1 + area2 + area3 + area4 + area5) / 2;
      
      //determine which region split will be within and resulting horizontal midpoint

      //Within area 1
      if( area1 > halfArea )
      {
        clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
      }
      //Within area 2
      else if( area1 + area2 > halfArea )
      {
        clusters[cluster].beta = ((halfArea - area1) / outM) + A;
      }
      //Within area 3-4
      else if( area1 + area2 + area3 + area4 > halfArea )
      {
        float a = 1.0f/3.0f;
        float b = outM - 0.5f - D/3;
        float c = area1 + area2 - D*outM + D/2 - halfArea;
        clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
      }
      //Within area5
      else
      {
        clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4) ) / outL) + C;
      }
    } //end case 2
    //---------------------------------------------------------------
    //Case 3
    else
    {
      //find area of 5 subregions
      float area1 = (outM*(A-0.5f)) / 2;
      float area2 = (B-A) * outM;
      float area3 = ((E-B) * (outM - outL)) / 2;
      float area4 = (E-B) * outL;
      float area5 = (3.0f - E) * outL;
      
      //find half of total area
      float halfArea = (area1 + area2 + area3 + area4 + area5) / 2;
      
      //determine which region split will be within and resulting horizontal midpoint
      
      //Within area 1
      if( area1 > halfArea )
      {
        clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
      }
      //Within area 2
      else if( area1 + area2 > halfArea )
      {
        clusters[cluster].beta = ((halfArea - area1) / outM) + A;
      }
      //Within area 3-4
      else if( area1 + area2 + area3 + area4 > halfArea )
      {
        float a = -0.5f;
        float b = E/2 + 2.5f/2; 
        float c = area3 - 2.5f*E/2;
        clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
      }
      //Within area5
      else
      {
        clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4) ) / outL) + E;
      }
    } //end case 3
    //---------------------------------------------------------------
    
    //Compute edge threshold
    int lumJND = 255/50;
    clusters[cluster].edgeThreshold = clusters[cluster].mode + clusters[cluster].beta*lumJND;

  } //end for cluster
  
}
//==============================================
void EdgeDetect::constructEdgeImage()
{
  int x, y;
  QRgb* rgb;
  
  uchar* scanLine;
  for( y=0; y<image->height(); y++)
  {   
    scanLine = image->scanLine(y);
    for( x=0; x<image->width(); x++)
    {
      //initialize pixel to black 
      rgb = ((QRgb*)scanLine+x);
      *rgb = qRgb( 0, 0, 0 );
      
      //lookup ESF for this pixel
      float ESF = LUT[ GSLCmap[x + y*image->width()] ].ESF;

      //If ESF value for this pixel is 0 skip
      if( ESF == 0.0f ) continue;

      //lookup edge magnitude threshold
      float lum = lumMap[x + y*image->width()];
      float edgeMagThresh = -1.0f;
      int cluster;
      for(cluster=0; cluster<numClusters; cluster++)
      {
        if(lum >= clusters[cluster].minLuminance &&
           lum <= clusters[cluster].maxLuminance)
        {
          edgeMagThresh = clusters[cluster].edgeThreshold;
          break;
        }
      }

      //if cluster not found bail
      if( cluster >= numClusters )
      {
//        cout << "Error! Could not find cluster pixel belonged to!\n";
        continue;
      }
      
      //if edge mag below thresh then skip
      if( edgeMagMap[x + y*image->width()] < edgeMagThresh ) continue;
        
      //ok, last checks implement NMS (non-maximum supression)
      int direction = LUT[ GSLCmap[x + y*image->width()] ].direction;
      int neighborIndex1 = -1;
      int neighborIndex2 = -1;

      if( direction == 0)
      {
        if( x > 0) 
          neighborIndex1 = x-1 + y*image->width();
        if( x < image->width() - 1 ) 
          neighborIndex2 = x+1 + y*image->width();
      }
      else if(direction == 1)
      {
        if( x > 0 && y < image->height() - 1 )
          neighborIndex1 = x-1 + (y+1)*image->width();
        if( x < image->width() - 1 && y > 0 )
          neighborIndex2 = x+1 + (y-1)*image->width();
      }
      else if(direction == 2)
      {
        if( y < image->height() - 1 ) 
          neighborIndex1 = x + (y+1)*image->width();
        if( y > 0) 
          neighborIndex2 = x + (y-1)*image->width();
      }
      else if(direction == 3)
      {
        if( x < image->width() - 1 && y < image->height() - 1 )
          neighborIndex1 = x+1 + (y+1)*image->width();
        if( x > 0 && y > 0 )
          neighborIndex2 = x-1 + (y-1)*image->width();
      }

      //neighbor 1 has higher confidence, skip!
      if( neighborIndex1 != -1 &&
          LUT[ GSLCmap[neighborIndex1] ].ESF * edgeMagMap[neighborIndex1] >
          ESF * edgeMagMap[x + y*image->width()] )
        continue;
      
      //neighbor 2 has higher confidence, skip!
      if( neighborIndex2 != -1 &&
          LUT[ GSLCmap[neighborIndex2] ].ESF * edgeMagMap[neighborIndex2] >
          ESF * edgeMagMap[x + y*image->width()] )
        continue;
      
      //All tests passed! Mark edge!
      *rgb = qRgb( 255, 255, 255 );
    } //x
  } //y
  
  //blur image - all of it
  blurImage( *image, 2.0f );

  //normalize image
  enhanceImageContrast( image );
  
}
//==============================================
void EdgeDetect::deallocateObjects()
{
  delete[] lumMap;      lumMap = NULL;
  delete[] edgeMagMap;  edgeMagMap = NULL;
  delete[] GSLCmap;     GSLCmap = NULL;
  delete[] clusters;    clusters = NULL;
}
//==============================================
void EdgeDetect::constructGSLClut()
{
  //----------------------
  //First fill entire table with 0 ESF's and invalid directions
  int i;
  for(i=0; i<256; i++)
  {
    LUT[i].ESF = 0.0f;
    LUT[i].direction = -1;
  }
  //----------------------
  //Next code in all pattern that are highly 
  //likely to be on edges as described in the paper
  //----------------------
  //Pattern (a)

  // ###
  // ##.
  // ...
  LUT[15].ESF = 0.179f;
  LUT[15].direction = 3;

  // ...
  // .##
  // ###
  LUT[240].ESF = 0.179f;
  LUT[240].direction = 3;

  // ###
  // .##
  // ...
  LUT[23].ESF = 0.179f;
  LUT[23].direction = 1;
  
  // ...
  // ##.
  // ###
  LUT[232].ESF = 0.179f;
  LUT[232].direction = 1;
   
  // ##.
  // ##.
  // #..
  LUT[43].ESF = 0.179f;
  LUT[43].direction = 3;
  
  // ..#
  // .##
  // .##
  LUT[212].ESF = 0.179f;
  LUT[212].direction = 3;

  // #..
  // ##.
  // ##.
  LUT[105].ESF = 0.179f;
  LUT[105].direction = 1;

  // .##
  // .##
  // ..#
  LUT[150].ESF = 0.179f;
  LUT[150].direction = 1;
  //----------------------
  //Pattern (b)

  // ###
  // ###
  // ...
  LUT[31].ESF = 0.137f;
  LUT[31].direction = 2;

  // ...
  // ###
  // ###
  LUT[248].ESF = 0.137f;
  LUT[248].direction = 2;
  
  // ##.
  // ##.
  // ##.
  LUT[107].ESF = 0.137f;
  LUT[107].direction = 0;
  
  // .##
  // .##
  // .##
  LUT[214].ESF = 0.137f;
  LUT[214].direction = 0;
  //----------------------
  //Pattern (c)
  
  // ###
  // .#.
  // ...
  LUT[7].ESF = 0.126f;
  LUT[7].direction = 2;
  
  // ...
  // .#.
  // ###
  LUT[224].ESF = 0.126f;
  LUT[224].direction = 2;

  // #..
  // ##.
  // #..
  LUT[41].ESF = 0.126f;
  LUT[41].direction = 0; 
  
  // ..#
  // .##
  // ..#
  LUT[148].ESF = 0.126f;
  LUT[148].direction = 0;
  //----------------------
  //Pattern (d)
  
  // ###
  // ##.
  // #..
  LUT[47].ESF = 0.10f;
  LUT[47].direction = 3;
   
  // ..#
  // .##
  // ###
  LUT[244].ESF = 0.10f;
  LUT[244].direction = 3;

  // ###
  // .##
  // ..#
  LUT[151].ESF = 0.10f;
  LUT[151].direction = 1;

  // #..
  // ##.
  // ###
  LUT[233].ESF = 0.10f;
  LUT[233].direction = 1;
  //----------------------
  //Pattern (e)
  
  // ##.
  // ##.
  // ...
  LUT[11].ESF = 0.10f;
  LUT[11].direction = 3;
  
  // ...
  // .##
  // .##
  LUT[208].ESF = 0.10f;
  LUT[208].direction = 3;

  // .##
  // .##
  // ...
  LUT[22].ESF = 0.10f;
  LUT[22].direction = 1;
  
  // ...
  // ##.
  // ##.
  LUT[104].ESF = 0.10f;
  LUT[104].direction = 1; 
  //----------------------    
}
//==============================================