edgeDetect.cpp

Go to the documentation of this file.
00001 //==============================================
00002 //  copyright            : (C) 2003-2005 by Will Stokes
00003 //==============================================
00004 //  This program is free software; you can redistribute it 
00005 //  and/or modify it under the terms of the GNU General 
00006 //  Public License as published by the Free Software 
00007 //  Foundation; either version 2 of the License, or  
00008 //  (at your option) any later version.         
00009 //==============================================
00010 
00011 //Systemwide includes
00012 #include <qpainter.h>
00013 #include <qimage.h>
00014 #include <math.h>
00015 
00016 //Projectwide includes
00017 #include "edgeDetect.h"
00018 #include "blur.h"
00019 #include "../enhancements/contrast.h"
00020 
00021 //----------------------------------------------
00022 // Inputs:
00023 // -------
00024 // QImage* image - image to find edges in
00025 //
00026 // Outputs:
00027 // --------
00028 // QImage* getEdgeImage - returns the produced grayscale edge image
00029 // 
00030 // Other information such as pixel groups etc is also availble through 
00031 // various accesor method.
00032 //
00033 // Description:
00034 // ------------
00035 // This class is the first known publically available implementation of 
00036 // Kim, Lee, and Kweon's 2003 paper titled:
00037 // "Automatic edge detection using 3x3 ideal binary 
00038 // pixel patterns and fuzzy-based edge thresholding"
00039 // http://rcv.kaist.ac.kr/publication/file/foreign_journal/28_DongSuKim_PRL2004.pdf
00040 //
00041 // Edge detection is an old problem, an while many use the edge detectors by 
00042 // Canny, Sobel, etc, they all suffer from a common problem: the user must 
00043 // tweak a series of poorly understood input parameters to get the ideal edge image.
00044 //
00045 // Album Shaper was in need of an automatic edge detector for use when 
00046 // blurring and sharpening images. Having the user manually tweak such 
00047 // paramters first would not only be annoying but error-prone. In an 
00048 // effort to make edge detection automatic I took a stab at implementing 
00049 // this paper and am quite happy with the resulsts...
00050 //
00051 // http://albumshaper.sourceforge.net/images/teasers/peaksAndValleys.jpg
00052 //
00053 // Algorithm:
00054 // ----------
00055 // While complex, the algorithm can be broken up into a series of 
00056 // fairly straightforward tasks:
00057 //
00058 // 1.) "allocateAndInitObjects()" is called to allocate and fill a 
00059 //     few data structures that will be used when finding image edges.
00060 //
00061 // 2.) "fillLumMapAndLumHistogram"() is called, during which an 
00062 //     luminance map is constructed and luminance histogram populated. 
00063 //     For an m x n image, the luminance map will be a m x n integer array.
00064 //
00065 // 3.) The luminance histogram is smoothed using "smoothLumHistogram" to 
00066 //     make peak finding less sensative to noise.
00067 //
00068 // 4.) The fourth step is a little complicated. The edge magnitude and 
00069 //     GSLC (Grey level similitude code) value is computed at each pixel.
00070 //     The paper takes an interesting approach to edge detection by 
00071 //     classifying pixels into one of 9 groups by first computing the
00072 //     average luminance for a 3x3 group centered about a pixel. Pixels are
00073 //     then separated into one of two groups, those that have a luminance 
00074 //     greater than or less than the 3x3 average luminance. For example:
00075 //                                   
00076 //                                                     X
00077 //      ---------          ---------        ---------X
00078 //     | 7 15 18 |        | 0  1  1 |      | 0  1  X |   
00079 //     | 5 17 20 |  -->   | 0  1  1 |  --> | 0  X  1 |
00080 //     | 9  8  3 |        | 0  0  0 |      | X  0  0 |
00081 //      ---------          ---------       X---------
00082 //                                       X
00083 //
00084 //     Here the average luminance is 11.333, placing the top right 4 
00085 //     pixels in one group and the other remaining pixels in another. 
00086 //     The dominant edge diretion is from the bottom left the top right. 
00087 //     The GLSC code is computing by considering the 1/0 values
00088 //     (1 = pixel in same group as central pixel). The central value 
00089 //     is ignored (it's always 1) leaving us with an 8bit = 2^8 = 256 code.
00090 //
00091 //     In this case the GSLC is:
00092 //     0*2^0 + 1*2^1 + 1*2^2 +
00093 //     0*2^3 + XXXXX + 1*2^4 +
00094 //     0*2^5 + 0*2^6 + 0*2^7 = 22.
00095 //
00096 //     The authors of the paper found pixels with one of five GSLC 
00097 //     codes (15,31,7,47,11) were most often associated with edge pixels 
00098 //     when producing an edge image using various competitive techinques. 
00099 //     By looking up the GSLC for a given pixel later one we can suppress 
00100 //     edges where they most likely do not belong.
00101 //
00102 // 5.) The fifth step involves grouping pixels by luminance using the 
00103 //     smoothed luminance histogram. This complicated step is brushed off 
00104 //     as being trivial in the paper. Since I struggled a bit with developing 
00105 //     an algorithm for this step, I'll explain my approach in detail 
00106 //     to avoid others suffering.
00107 //
00108 //     Using the smoothed histogram, we first compute the JND or just 
00109 //     noticible differnce using the maximum count that was found. I'm 
00110 //     not sure how appropriate this is, but 2% is the usual quantity
00111 //     used in other contexts, and it works well here, so 2% it is.
00112 // 
00113 //     Next we walk through the smoothed luminance histogram and find 
00114 //     the midpoint of the valleys. We accomplish this by updating an 
00115 //     index of the deepend last known valley. As the valley slopes
00116 //     down and we move across it we update this last best known index. 
00117 //     Once the valley slopes up one JND above the last deepend location 
00118 //     found we mark that valley midpoint and move on.
00119 //
00120 //     Once all valley midpoints have been marked we can quickly deduce 
00121 //     how many peaks must exist. We pass across the smoothed luminance 
00122 //     histogram again finding the peak index for each pixel between
00123 //     valley midpoints. For all future work pixels one JND+- the peak 
00124 //     center will be used.
00125 //
00126 // 6.) The sixth step, "computeClusterStatistics()", computes various 
00127 //     cluster-specific statsitics that will be used to determine cluster 
00128 //     thresholds. The paper was rather vague in this area, but after 
00129 //     experimenting with various interpreatations of what they were trying 
00130 //     to say I think I finally got it right.
00131 // 
00132 //     First, we iterate over all image pixels, determine which pixel group 
00133 //     they belong by comparing luminance endpoints for all pixelclusters, 
00134 //     and update total edge mag, num pixels, and an edgeMagHistogram for 
00135 //     the given pixel group they belong to.
00136 //
00137 //     Next we compute the average edge meganitude and most frequent edge 
00138 //     magnitude observed for each pixel cluster, in addition to normalizing 
00139 //     the cluster pixel count variable to [0,1]
00140 //
00141 // 7.) The seventh step is quite complicated and encompases computing the 
00142 //     edge thresholds for each pixel cluster using the 18-rule fuzzy 
00143 //     logic approach put forward by the paper. There is nothing ground 
00144 //     breaking here, just a lot of complicated fuzzy logic, although most 
00145 //     of the effort is put into computing the centroid at the end. I had 
00146 //     never touched fuzzy logic before, but found this article more than
00147 //     helpful getting myself up to speed:
00148 //
00149 //     http://www.doc.ic.ac.uk/~nd/surprise_96/journal/vol2/sbaa/article2.html
00150 //
00151 // 8.) The eigth and final step, actually constructing the edge image, is fairly
00152 //     straight forward and employs techinques (such as non-maximum suppression or NMS) 
00153 //     anyone who has implemented Canny shoudl be familiar with.
00154 //
00155 //     First, one looks up the ESF (edge shape factor) for a given pixel from a look-up 
00156 //     table. These values were computed by generating a ton of edge images by carefully
00157 //     setting Canny, Sobel, etc params, then identifying how often a pixel
00158 //     with a given GSLC code was judged to be an edge. Hence ESF's will fall in
00159 //     the [0,1] range and help to suppress edges where no clear edge really
00160 //     can be claimed to exist, e.g.
00161 //
00162 //     . # .
00163 //     # # #   <- If there is an edge here, care to explain what the
00164 //     . # .      edge direction is?!!
00165 //
00166 //
00167 //     If the ESF for a give pixel is 0 we skip it, it's not an edge.
00168 //
00169 //     Next, we look up a pixels edge magnitude threshold by identifying the
00170 //     pixel cluster is belongs to using the pixels luminance using the
00171 //     luminance map.
00172 // 
00173 //     If the pixels edge magnitude is less than the threshold we skip the pixel, 
00174 //     this filtersout low lying noise.
00175 //
00176 //     Finally, the direction of the pixel is looked up using its GSLC and
00177 //     NMS (non-maximum suppression is applied). If the pixel has a greater
00178 //     egdge magnitude than either of its neighbors along the edge direction then
00179 //     the edge is marked.
00180 //
00181 // Final Remarks:
00182 // --------------
00183 // Despite the involved complexity, the implementation appears to work
00184 // really really well. I consider this one of the secret gems of Album Shaper
00185 // and hope to make good use of it in the future for things other than 
00186 // just sharpening and bluring. The only caveat is that edge detection 
00187 // does take a few CPU cycles.
00188 //----------------------------------------------
00189 
00190 
00191 //==============================================
00192 EdgeDetect::EdgeDetect( QImage* image )
00193 {                     
00194   //load image
00195   this->image = image;
00196 
00197   //allocate and initialize objects used for edge detection
00198   allocateAndInitObjects();
00199 
00200   //fill lum map and lum histogram
00201   fillLumMapAndLumHistogram();
00202   
00203   //fill smoothed histogram
00204   smoothLumHistogram();
00205   
00206   //compute edge magnitude and GSLC maps
00207   computeEdgeMagAndGSLCmaps();
00208   
00209   //determine pixel clusters
00210   findPixelClusters();
00211   
00212   computeClusterStatistics();  
00213  
00214   computeClusterThresholds();
00215   
00216   constructEdgeImage();
00217 }
00218 //==============================================
00219 EdgeDetect::~EdgeDetect()
00220 {
00221   deallocateObjects();
00222 }
00223 //==============================================
00224 int EdgeDetect::getNumClusters()
00225 { return numClusters; }
00226 //==============================================
00227 PixelCluster* EdgeDetect::getClusters()
00228 { return clusters; }
00229 //==============================================
00230 int* EdgeDetect::getPeaks()
00231 { return clusterPeaks; }
00232 //==============================================
00233 int* EdgeDetect::getSmoothHist()
00234 { return smoothLumHist; }
00235 //==============================================
00236 QImage* EdgeDetect::getEdgeImage()
00237 {
00238   return image; 
00239 }
00240 //==============================================
00241 int* EdgeDetect::getClusterMap()
00242 {
00243   //construct map
00244   int* clusterMap = new int[image->width() * image->height()];
00245   
00246   //iterate over all pixels, determine cluster each pixel belongs to
00247   int i, cluster;
00248   for(i=0; i<image->width()*image->height(); i++)
00249   {
00250     for(cluster=0; cluster<numClusters; cluster++)
00251     {
00252       if( lumMap[i] >= clusters[cluster].minLuminance &&
00253           lumMap[i] <= clusters[cluster].maxLuminance )
00254       {
00255         clusterMap[i] = cluster;
00256         break;
00257       }
00258     } //cluster
00259   } //pixel
00260 
00261   return clusterMap;
00262 }
00263 //==============================================
00264 void EdgeDetect::allocateAndInitObjects()
00265 {
00266   //initialize: 
00267   //-luminosity histogram
00268   //-smoothed luminosity histogram
00269   //-identified peak regions
00270   int i;
00271   for(i=0; i<256; i++)
00272   { 
00273     lumHist[i] = 0; 
00274     smoothLumHist[i] = 0;
00275     clusterPeaks[i] = 0;
00276   }
00277   
00278   //allocate luminance map
00279   lumMap = new int[image->width() * image->height()];
00280   
00281   //allocate edge magnitude map
00282   edgeMagMap = new float[image->width() * image->height()];
00283   
00284   //allocate GSLC map
00285   GSLCmap = new int[image->width() * image->height()];
00286   
00287   //construct LUT
00288   constructGSLClut();
00289 }
00290 //==============================================
00291 void EdgeDetect::fillLumMapAndLumHistogram()
00292 {
00293   int x, y;
00294   QRgb* rgb;
00295   uchar* scanLine;
00296   int lumVal;
00297   for( y=0; y<image->height(); y++)
00298   {   
00299     scanLine = image->scanLine(y);
00300     for( x=0; x<image->width(); x++)
00301     {
00302       //get lum value for this pixel
00303       rgb = ((QRgb*)scanLine+x);
00304       lumVal = qGray(*rgb);
00305       
00306       //store in lum map
00307       lumMap[x + y*image->width()] = lumVal;
00308       
00309       //update lum histogram
00310       lumHist[ lumVal ]++;
00311     }
00312   }
00313 }
00314 //==============================================
00315 void EdgeDetect::smoothLumHistogram()
00316 {
00317   #define FILTER_SIZE 5
00318   int filter[FILTER_SIZE] = {2, 5, 8, 5, 2};
00319   
00320   int i,j;
00321   int filterIndex, sum, total;
00322   for(i = 0; i<256; i++)
00323   {
00324     sum = 0;
00325     total = 0;
00326     
00327     for( j= -FILTER_SIZE/2; j <= FILTER_SIZE/2; j++)
00328     {
00329       if( i+j > 0 && i+j < 256 )
00330       {
00331         filterIndex = j+ FILTER_SIZE/2;
00332         total+= filter[filterIndex] * lumHist[i+j];
00333         sum  += filter[filterIndex];
00334       }
00335     }
00336     
00337     smoothLumHist[i] = total / sum;
00338   }
00339 }
00340 //==============================================
00341 void EdgeDetect::computeEdgeMagAndGSLCmaps()
00342 {
00343   int x, y;
00344   int idealPattern[9];
00345   int pixelLums[9];
00346   
00347   //-------  
00348   //iterate over all pixels
00349   for( y=0; y<image->height(); y++)
00350   {   
00351     for( x=0; x<image->width(); x++)
00352     {
00353       //compute pixel luminances for entire grid
00354       pixelLums[0] = pixelLum(x-1,y-1);
00355       pixelLums[1] = pixelLum(x  ,y-1);
00356       pixelLums[2] = pixelLum(x+1,y-1);
00357       pixelLums[3] = pixelLum(x-1,y  );
00358       pixelLums[4] = pixelLum(x  ,y  );
00359       pixelLums[5] = pixelLum(x+1,y  );
00360       pixelLums[6] = pixelLum(x-1,y+1);
00361       pixelLums[7] = pixelLum(x  ,y+1);
00362       pixelLums[8] = pixelLum(x+1,y+1);
00363       
00364       //compute average
00365       float avg = 0;
00366       int i;
00367       for(i=0; i<=8; i++)
00368       {
00369         avg+= pixelLums[i];
00370       }
00371       avg = avg / 9;
00372       
00373       //determine ideal pattern and I0 and I1 averages
00374       int centerPixelLum = pixelLums[4];
00375       float centerDiff = centerPixelLum - avg;
00376       
00377       float I0avg = 0;
00378       int I0count = 0;
00379       
00380       float I1avg = 0;
00381       int I1count = 0;
00382       
00383       for(i=0; i<=8; i++)
00384       {
00385         if( centerDiff * (pixelLums[i]-avg) >=0 )
00386         { 
00387           I1avg+=pixelLums[i];
00388           I1count++;
00389           idealPattern[i] = 1; 
00390         }
00391         else 
00392         { 
00393           I0avg+=pixelLums[i];
00394           I0count++;
00395           idealPattern[i] = 0; 
00396         }
00397       }
00398 
00399       //compute and store edge magnitude
00400       if(I0count > 0) I0avg = I0avg/I0count;
00401       if(I1count > 0) I1avg = I1avg/I1count;     
00402       edgeMagMap[x + y*image->width()] = QABS( I1avg - I0avg );
00403       
00404       //compute and store GSLC
00405       int GSLC=0;
00406       int weight = 1;
00407       for(i=0; i<9; i++)
00408       {
00409         //skip center
00410         if(i == 4) continue;
00411         
00412         if(idealPattern[i] == 1)
00413         { GSLC+=weight; }
00414         
00415         weight = weight*2;
00416       }
00417       GSLCmap[x + y*image->width()] = GSLC;
00418     } //x
00419   } //y
00420 }
00421 //==============================================
00422 int EdgeDetect::pixelLum(int x, int y)
00423 {
00424   int clampedX = QMAX( QMIN( x, image->width()-1), 0);
00425   int clampedY = QMAX( QMIN( y, image->height()-1), 0);
00426   return lumMap[ clampedX + clampedY * image->width() ];
00427 }
00428 //==============================================
00429 void EdgeDetect::findPixelClusters()
00430 {
00431   //find max count
00432   int maxCount = 0;
00433   int i;
00434   for(i=0; i<256; i++)
00435   {
00436     if(smoothLumHist[i] > maxCount)
00437       maxCount = smoothLumHist[i];
00438   }
00439 
00440   //compute JND for histogram (2% of total spread)
00441   int histJND = maxCount/50;
00442 
00443   //construct temporary array for valley locations
00444   //1's will indicate a valley midpoint
00445   int tmpValleyArray[256];
00446   for(i=0; i<256; i++) { tmpValleyArray[i] = 0; }
00447   
00448   //move across histogram finding valley midpoints
00449   int curTrackedMin = smoothLumHist[0];
00450   
00451   //first and last indices tracked min was observed
00452   int firstMinIndex = 0;
00453   int lastMinIndex = 0;
00454   
00455   //only add valley midpoint if finished tracking a descent
00456   bool slopeNeg = false;
00457   
00458   for(i = 1; i<256; i++ )
00459   {
00460     if( smoothLumHist[i] < curTrackedMin - histJND )
00461     {
00462       //found a descent!
00463       slopeNeg = true;
00464       curTrackedMin = smoothLumHist[i];
00465       firstMinIndex = i;
00466     }
00467     //starting to go up again, add last min to list
00468     else if( smoothLumHist[i] > curTrackedMin + histJND )
00469     {
00470       //if finished tracing a negative slope find midpoint and set location to true
00471       if(slopeNeg)
00472       {
00473         tmpValleyArray[ (firstMinIndex + lastMinIndex)/2 ] = 1;
00474       }
00475       
00476       curTrackedMin = smoothLumHist[i];
00477       slopeNeg = false;
00478     }
00479     else
00480     {
00481       //still tracking a min, update the right 
00482       //hand index. center of valley is found
00483       //by averaging first and last min index
00484       lastMinIndex = i;
00485     }
00486   }
00487   
00488   //count valleys
00489   int numValleys = 0;
00490   for(i=0; i<256; i++)
00491   {
00492     if(tmpValleyArray[i] == 1 ) numValleys++;
00493   }
00494 
00495   //determine number of clusters
00496   numClusters = numValleys-1;
00497   if(tmpValleyArray[0] != 1)
00498     numClusters++;
00499   if(tmpValleyArray[255] != 1)
00500     numClusters++;
00501   
00502   //allocate clusters
00503   clusters = new PixelCluster[numClusters];
00504   
00505   //automatically start first cluster
00506   int cluster=0;
00507   clusters[cluster].minLuminance = 0;
00508   
00509   //initialize left and right boundaries of all clusters
00510   for(i=1; i<256; i++)
00511   {
00512     //reached next valley, end cluster
00513     if( tmpValleyArray[i] == 1)
00514     {
00515       clusters[cluster].maxLuminance = i-1;
00516       cluster++;
00517       clusters[cluster].minLuminance = i;
00518     }
00519     //end last cluster automatically at end
00520     else if(i == 255)
00521     {
00522       clusters[cluster].maxLuminance = i;
00523     }
00524   }
00525   
00526   //determine cluster peaks
00527   for(cluster=0; cluster<numClusters; cluster++)
00528   {
00529     //find max for current cluster
00530     int maxIndex = clusters[cluster].minLuminance;
00531     for(i=clusters[cluster].minLuminance; i<=clusters[cluster].maxLuminance; i++)
00532     {
00533       if(smoothLumHist[i] > smoothLumHist[maxIndex])
00534          maxIndex = i;
00535     }
00536     
00537     //mark peaks  
00538     int lumJND = 255/50;
00539     for(i=QMAX(0, maxIndex-lumJND); i<QMIN(256, maxIndex+lumJND); i++)
00540     { 
00541       clusterPeaks[i] = 1; 
00542     }
00543   }
00544 }
00545 //==============================================
00546 void EdgeDetect::computeClusterStatistics()
00547 {
00548   //initialize cluster stats
00549   int cluster;
00550   for(cluster=0; cluster<numClusters; cluster++)
00551   {
00552     int i;
00553     for(i=0; i<256; i++)
00554     {
00555       clusters[cluster].edgeMagHistogram[i] = 0;
00556     }
00557     clusters[cluster].totalEdgeMagnitude=0.0f;
00558     clusters[cluster].numPixels = 0;
00559   }
00560   
00561   //iterate over all pixels
00562   int i;
00563   for(i=0; i<image->width()*image->height(); i++)
00564   {
00565     //skip pixels that don't belong to peaks
00566     if( clusterPeaks[ lumMap[i] ] != 1)
00567       continue;
00568     
00569     //determine cluster pixel belongs to
00570     int cluster;
00571     for(cluster=0; cluster<numClusters; cluster++)
00572     {
00573       if( lumMap[i] >= clusters[cluster].minLuminance &&
00574           lumMap[i] <= clusters[cluster].maxLuminance )
00575       {      
00576         clusters[cluster].totalEdgeMagnitude+= edgeMagMap[i]; 
00577         clusters[cluster].numPixels++;
00578         clusters[cluster].edgeMagHistogram[ QMIN( QMAX( (int)edgeMagMap[i], 0), 255) ]++;
00579         break;
00580       }
00581     } //cluster
00582   } //pixel i
00583   
00584   //iterate over clusters to determine min and max peak cluster sizes
00585   minClusterSize = clusters[0].numPixels;
00586   maxClusterSize = clusters[0].numPixels;
00587   for(cluster=1; cluster<numClusters; cluster++)
00588   {
00589     if(clusters[cluster].numPixels < minClusterSize)
00590       minClusterSize = clusters[cluster].numPixels;
00591 
00592     if(clusters[cluster].numPixels > maxClusterSize)
00593       maxClusterSize = clusters[cluster].numPixels;
00594   }
00595 
00596   //iterate over clusters one final time to deduce normalized inputs to fuzzy logic process
00597   int JND = 255/50;
00598   for(cluster=0; cluster<numClusters; cluster++)
00599   {
00600     clusters[cluster].meanMode = QMIN( clusters[cluster].totalEdgeMagnitude / clusters[cluster].numPixels,
00601                                        3*JND );
00602     
00603     int i;
00604     int mode = 0;
00605     for(i=1; i<256; i++)
00606     {
00607       if( clusters[cluster].edgeMagHistogram[i] > clusters[cluster].edgeMagHistogram[ mode ] )
00608         mode = i;
00609     }
00610     clusters[cluster].mode = QMIN( mode, 2*JND );
00611         
00612     clusters[cluster].pixelCount = ((float)(clusters[cluster].numPixels - minClusterSize)) / 
00613                                    (maxClusterSize - minClusterSize);  
00614   }
00615 }
00616 //==============================================
00617 //compute edge thresholds for each cluster using 18-rule fuzzy logic approach
00618 void EdgeDetect::computeClusterThresholds()
00619 {
00620   //iterate over each cluster
00621   int cluster;
00622   float S1,M1,L1;
00623   float S2,M2,L2;
00624   float S3,L3;
00625   float outS, outM, outL;
00626   
00627   int JND = 255/50;
00628   
00629   for(cluster=0; cluster<numClusters; cluster++)
00630   {
00631     //----
00632     //compute S,M, and L values for each input
00633     //----
00634     S1 = QMAX( 1.0f - ((clusters[cluster].meanMode/JND) / 1.5f), 0 );
00635 
00636     if( (clusters[cluster].meanMode/JND) <= 1.5f )
00637       M1 = QMAX( (clusters[cluster].meanMode/JND) - 0.5f, 0 );
00638     else
00639       M1 = QMAX( 2.5f - (clusters[cluster].meanMode/JND), 0 );
00640     
00641     L1 = QMAX( ((clusters[cluster].meanMode/JND) - 1.5f) / 1.5f, 0 );
00642     //----
00643     S2 = QMAX( 1.0f - (clusters[cluster].mode/JND), 0 );
00644     
00645     if( (clusters[cluster].mode/JND) <= 1.0f )
00646       M2 = QMAX( -1.0f + 2*(clusters[cluster].mode/JND), 0 );
00647     else
00648       M2 = QMAX( 3.0f - 2*(clusters[cluster].mode/JND), 0 );
00649     
00650     L2 = QMAX( (clusters[cluster].mode/JND) - 1.0, 0 );
00651     //----
00652     S3 = QMAX( 1.0f - 2*clusters[cluster].pixelCount, 0 );
00653     L3 = QMAX( -1.0f + 2*clusters[cluster].pixelCount, 0 );
00654     //----
00655     
00656     //Compute M,L for outputs using set of 18 rules.
00657     //outS is inherantly S given the ruleset provided
00658     outS = 0.0f;
00659     outM = 0.0f;
00660     outL = 0.0f;
00661     //Out 1
00662     if( numClusters > 2 )
00663     {
00664       outM += S1*S2*S3;   //rule 1
00665       
00666       //rule 2
00667       if( clusters[cluster].meanMode < clusters[cluster].mode )
00668         outS += S1*S2*L3;   
00669       else
00670         outM += S1*S2*L3;
00671 
00672       outM += S1*M2*S3;   //rule 3
00673       outM += S1*M2*L3;   //rule 4
00674       outM += S1*L2*S3;   //rule 5
00675       outM += S1*L2*L3;   //rule 6
00676       outM += M1*S2*S3;   //rule 7
00677       outM += M1*S2*L3;   //rule 8
00678       outM += M1*M2*S3;   //rule 9
00679       outL += M1*M2*L3;   //rule 10
00680       outM += M1*L2*S3;   //rule 11
00681       outL += M1*L2*L3;   //rule 12
00682       outM += L1*S2*S3;   //rule 13
00683       outL += L1*S2*L3;   //rule 14
00684       outM += L1*M2*S3;   //rule 15
00685       outL += L1*M2*L3;   //rule 16
00686       outL += L1*L2*S3;   //rule 17
00687       outL += L1*L2*L3;   //rule 18
00688     }
00689     //Out 2
00690     else
00691     {
00692       outL += S1*S2*S3;   //rule 1
00693       outL += S1*S2*L3;   //rule 2
00694       outM += S1*M2*S3;   //rule 3
00695       outL += S1*M2*L3;   //rule 4
00696       outM += S1*L2*S3;   //rule 5
00697       outM += S1*L2*L3;   //rule 6
00698       outL += M1*S2*S3;   //rule 7
00699       outL += M1*S2*L3;   //rule 8
00700       outL += M1*M2*S3;   //rule 9
00701       outL += M1*M2*L3;   //rule 10
00702       outL += M1*L2*S3;   //rule 11
00703       outL += M1*L2*L3;   //rule 12
00704       outL += L1*S2*S3;   //rule 13
00705       outL += L1*S2*L3;   //rule 14
00706       outL += L1*M2*S3;   //rule 15
00707       outL += L1*M2*L3;   //rule 16
00708       outL += L1*L2*S3;   //rule 17
00709       outL += L1*L2*L3;   //rule 18
00710     }
00711 
00712     //find centroid - Beta[k]
00713     float A = outM + 0.5f;
00714     float B = 2.5f - outM;
00715     float C = 1.5f * (outL + 1);
00716     float D = 1.5f * (outM + 1);
00717     float E = 2.5f - outL;
00718 
00719     //---------------------------------------------------------------
00720     //Case 1: Both outM and outL are above intersection point of diagonals
00721     if( outM > 0.5f && outL > 0.5f )
00722     {
00723       //find area of 7 subregions
00724       float area1 = ((A-0.5f)*outM)/2;
00725       float area2 = outM * (B-A);
00726       float area3 = ((2.1f-B) * (outM - 0.5)) / 2;
00727       float area4 = (2.1 - B) * 0.5f;
00728       float area5 = ((C - 2.1f) * (outL - 0.5)) / 2;
00729       float area6 = (C - 2.1f) * 0.5f;
00730       float area7 = (3.0f - C) * outL;
00731      
00732       //find half of total area
00733       float halfArea = (area1 + area2 + area3 + area4 + area5 + area6 + area7) / 2;
00734       
00735       //determine which region split will be within and resulting horizontal midpoint
00736       
00737       //Within area 1
00738       if( area1 > halfArea )
00739       {
00740         clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);       
00741       }
00742       //Within area 2
00743       else if( area1 + area2 > halfArea )
00744       {
00745         clusters[cluster].beta = ((halfArea - area1) / outM) + A;        
00746       }
00747       //Within area 3-4
00748       else if( area1 + area2 + area3 + area4 > halfArea )
00749       {
00750         float a = -0.5f;
00751         float b = 2.8f;
00752         float c = area1 + area2 + area3 - halfArea - B/2 - 2.625f;
00753         clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
00754       }
00755       //Within area 5-6
00756       else if( area1 + area2 + area3 + area4 + area5 + area6 > halfArea )
00757       {
00758         float a = 1.0f/3;
00759         float b = -0.7f;
00760         float c = area1 + area2 + area3 + area4 - halfArea;
00761         clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
00762       }
00763       //Within area 7
00764       else
00765       {
00766         clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4 + area5 + area6) ) / outL) + C;       
00767       }
00768     } //end case 1
00769     //---------------------------------------------------------------
00770     //Case 2
00771     else if ( outM < 0.5f && outL > outM )
00772     {
00773       //find area of 5 subregions
00774       float area1 = (outM*(A-0.5f)) / 2;
00775       float area2 = (D-A) * outM;
00776       float area3 = ((C-D) * (outL - outM)) / 2;
00777       float area4 = (C-D) * outM;
00778       float area5 = (3.0f - C) * outL;
00779       
00780       //find half of total area
00781       float halfArea = (area1 + area2 + area3 + area4 + area5) / 2;
00782       
00783       //determine which region split will be within and resulting horizontal midpoint
00784 
00785       //Within area 1
00786       if( area1 > halfArea )
00787       {
00788         clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
00789       }
00790       //Within area 2
00791       else if( area1 + area2 > halfArea )
00792       {
00793         clusters[cluster].beta = ((halfArea - area1) / outM) + A;
00794       }
00795       //Within area 3-4
00796       else if( area1 + area2 + area3 + area4 > halfArea )
00797       {
00798         float a = 1.0f/3.0f;
00799         float b = outM - 0.5f - D/3;
00800         float c = area1 + area2 - D*outM + D/2 - halfArea;
00801         clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
00802       }
00803       //Within area5
00804       else
00805       {
00806         clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4) ) / outL) + C;
00807       }
00808     } //end case 2
00809     //---------------------------------------------------------------
00810     //Case 3
00811     else
00812     {
00813       //find area of 5 subregions
00814       float area1 = (outM*(A-0.5f)) / 2;
00815       float area2 = (B-A) * outM;
00816       float area3 = ((E-B) * (outM - outL)) / 2;
00817       float area4 = (E-B) * outL;
00818       float area5 = (3.0f - E) * outL;
00819       
00820       //find half of total area
00821       float halfArea = (area1 + area2 + area3 + area4 + area5) / 2;
00822       
00823       //determine which region split will be within and resulting horizontal midpoint
00824       
00825       //Within area 1
00826       if( area1 > halfArea )
00827       {
00828         clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
00829       }
00830       //Within area 2
00831       else if( area1 + area2 > halfArea )
00832       {
00833         clusters[cluster].beta = ((halfArea - area1) / outM) + A;
00834       }
00835       //Within area 3-4
00836       else if( area1 + area2 + area3 + area4 > halfArea )
00837       {
00838         float a = -0.5f;
00839         float b = E/2 + 2.5f/2; 
00840         float c = area3 - 2.5f*E/2;
00841         clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
00842       }
00843       //Within area5
00844       else
00845       {
00846         clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4) ) / outL) + E;
00847       }
00848     } //end case 3
00849     //---------------------------------------------------------------
00850     
00851     //Compute edge threshold
00852     int lumJND = 255/50;
00853     clusters[cluster].edgeThreshold = clusters[cluster].mode + clusters[cluster].beta*lumJND;
00854 
00855   } //end for cluster
00856   
00857 }
00858 //==============================================
00859 void EdgeDetect::constructEdgeImage()
00860 {
00861   int x, y;
00862   QRgb* rgb;
00863   
00864   uchar* scanLine;
00865   for( y=0; y<image->height(); y++)
00866   {   
00867     scanLine = image->scanLine(y);
00868     for( x=0; x<image->width(); x++)
00869     {
00870       //initialize pixel to black 
00871       rgb = ((QRgb*)scanLine+x);
00872       *rgb = qRgb( 0, 0, 0 );
00873       
00874       //lookup ESF for this pixel
00875       float ESF = LUT[ GSLCmap[x + y*image->width()] ].ESF;
00876 
00877       //If ESF value for this pixel is 0 skip
00878       if( ESF == 0.0f ) continue;
00879 
00880       //lookup edge magnitude threshold
00881       float lum = lumMap[x + y*image->width()];
00882       float edgeMagThresh = -1.0f;
00883       int cluster;
00884       for(cluster=0; cluster<numClusters; cluster++)
00885       {
00886         if(lum >= clusters[cluster].minLuminance &&
00887            lum <= clusters[cluster].maxLuminance)
00888         {
00889           edgeMagThresh = clusters[cluster].edgeThreshold;
00890           break;
00891         }
00892       }
00893 
00894       //if cluster not found bail
00895       if( cluster >= numClusters )
00896       {
00897 //        cout << "Error! Could not find cluster pixel belonged to!\n";
00898         continue;
00899       }
00900       
00901       //if edge mag below thresh then skip
00902       if( edgeMagMap[x + y*image->width()] < edgeMagThresh ) continue;
00903         
00904       //ok, last checks implement NMS (non-maximum supression)
00905       int direction = LUT[ GSLCmap[x + y*image->width()] ].direction;
00906       int neighborIndex1 = -1;
00907       int neighborIndex2 = -1;
00908 
00909       if( direction == 0)
00910       {
00911         if( x > 0) 
00912           neighborIndex1 = x-1 + y*image->width();
00913         if( x < image->width() - 1 ) 
00914           neighborIndex2 = x+1 + y*image->width();
00915       }
00916       else if(direction == 1)
00917       {
00918         if( x > 0 && y < image->height() - 1 )
00919           neighborIndex1 = x-1 + (y+1)*image->width();
00920         if( x < image->width() - 1 && y > 0 )
00921           neighborIndex2 = x+1 + (y-1)*image->width();
00922       }
00923       else if(direction == 2)
00924       {
00925         if( y < image->height() - 1 ) 
00926           neighborIndex1 = x + (y+1)*image->width();
00927         if( y > 0) 
00928           neighborIndex2 = x + (y-1)*image->width();
00929       }
00930       else if(direction == 3)
00931       {
00932         if( x < image->width() - 1 && y < image->height() - 1 )
00933           neighborIndex1 = x+1 + (y+1)*image->width();
00934         if( x > 0 && y > 0 )
00935           neighborIndex2 = x-1 + (y-1)*image->width();
00936       }
00937 
00938       //neighbor 1 has higher confidence, skip!
00939       if( neighborIndex1 != -1 &&
00940           LUT[ GSLCmap[neighborIndex1] ].ESF * edgeMagMap[neighborIndex1] >
00941           ESF * edgeMagMap[x + y*image->width()] )
00942         continue;
00943       
00944       //neighbor 2 has higher confidence, skip!
00945       if( neighborIndex2 != -1 &&
00946           LUT[ GSLCmap[neighborIndex2] ].ESF * edgeMagMap[neighborIndex2] >
00947           ESF * edgeMagMap[x + y*image->width()] )
00948         continue;
00949       
00950       //All tests passed! Mark edge!
00951       *rgb = qRgb( 255, 255, 255 );
00952     } //x
00953   } //y
00954   
00955   //blur image - all of it
00956   blurImage( *image, 2.0f );
00957 
00958   //normalize image
00959   enhanceImageContrast( image );
00960   
00961 }
00962 //==============================================
00963 void EdgeDetect::deallocateObjects()
00964 {
00965   delete[] lumMap;      lumMap = NULL;
00966   delete[] edgeMagMap;  edgeMagMap = NULL;
00967   delete[] GSLCmap;     GSLCmap = NULL;
00968   delete[] clusters;    clusters = NULL;
00969 }
00970 //==============================================
00971 void EdgeDetect::constructGSLClut()
00972 {
00973   //----------------------
00974   //First fill entire table with 0 ESF's and invalid directions
00975   int i;
00976   for(i=0; i<256; i++)
00977   {
00978     LUT[i].ESF = 0.0f;
00979     LUT[i].direction = -1;
00980   }
00981   //----------------------
00982   //Next code in all pattern that are highly 
00983   //likely to be on edges as described in the paper
00984   //----------------------
00985   //Pattern (a)
00986 
00987   // ###
00988   // ##.
00989   // ...
00990   LUT[15].ESF = 0.179f;
00991   LUT[15].direction = 3;
00992 
00993   // ...
00994   // .##
00995   // ###
00996   LUT[240].ESF = 0.179f;
00997   LUT[240].direction = 3;
00998 
00999   // ###
01000   // .##
01001   // ...
01002   LUT[23].ESF = 0.179f;
01003   LUT[23].direction = 1;
01004   
01005   // ...
01006   // ##.
01007   // ###
01008   LUT[232].ESF = 0.179f;
01009   LUT[232].direction = 1;
01010    
01011   // ##.
01012   // ##.
01013   // #..
01014   LUT[43].ESF = 0.179f;
01015   LUT[43].direction = 3;
01016   
01017   // ..#
01018   // .##
01019   // .##
01020   LUT[212].ESF = 0.179f;
01021   LUT[212].direction = 3;
01022 
01023   // #..
01024   // ##.
01025   // ##.
01026   LUT[105].ESF = 0.179f;
01027   LUT[105].direction = 1;
01028 
01029   // .##
01030   // .##
01031   // ..#
01032   LUT[150].ESF = 0.179f;
01033   LUT[150].direction = 1;
01034   //----------------------
01035   //Pattern (b)
01036 
01037   // ###
01038   // ###
01039   // ...
01040   LUT[31].ESF = 0.137f;
01041   LUT[31].direction = 2;
01042 
01043   // ...
01044   // ###
01045   // ###
01046   LUT[248].ESF = 0.137f;
01047   LUT[248].direction = 2;
01048   
01049   // ##.
01050   // ##.
01051   // ##.
01052   LUT[107].ESF = 0.137f;
01053   LUT[107].direction = 0;
01054   
01055   // .##
01056   // .##
01057   // .##
01058   LUT[214].ESF = 0.137f;
01059   LUT[214].direction = 0;
01060   //----------------------
01061   //Pattern (c)
01062   
01063   // ###
01064   // .#.
01065   // ...
01066   LUT[7].ESF = 0.126f;
01067   LUT[7].direction = 2;
01068   
01069   // ...
01070   // .#.
01071   // ###
01072   LUT[224].ESF = 0.126f;
01073   LUT[224].direction = 2;
01074 
01075   // #..
01076   // ##.
01077   // #..
01078   LUT[41].ESF = 0.126f;
01079   LUT[41].direction = 0; 
01080   
01081   // ..#
01082   // .##
01083   // ..#
01084   LUT[148].ESF = 0.126f;
01085   LUT[148].direction = 0;
01086   //----------------------
01087   //Pattern (d)
01088   
01089   // ###
01090   // ##.
01091   // #..
01092   LUT[47].ESF = 0.10f;
01093   LUT[47].direction = 3;
01094    
01095   // ..#
01096   // .##
01097   // ###
01098   LUT[244].ESF = 0.10f;
01099   LUT[244].direction = 3;
01100 
01101   // ###
01102   // .##
01103   // ..#
01104   LUT[151].ESF = 0.10f;
01105   LUT[151].direction = 1;
01106 
01107   // #..
01108   // ##.
01109   // ###
01110   LUT[233].ESF = 0.10f;
01111   LUT[233].direction = 1;
01112   //----------------------
01113   //Pattern (e)
01114   
01115   // ##.
01116   // ##.
01117   // ...
01118   LUT[11].ESF = 0.10f;
01119   LUT[11].direction = 3;
01120   
01121   // ...
01122   // .##
01123   // .##
01124   LUT[208].ESF = 0.10f;
01125   LUT[208].direction = 3;
01126 
01127   // .##
01128   // .##
01129   // ...
01130   LUT[22].ESF = 0.10f;
01131   LUT[22].direction = 1;
01132   
01133   // ...
01134   // ##.
01135   // ##.
01136   LUT[104].ESF = 0.10f;
01137   LUT[104].direction = 1; 
01138   //----------------------    
01139 }
01140 //==============================================

Generated on Thu Jan 3 10:52:46 2008 for AlbumShaper by  doxygen 1.5.4