AlbumShaper  1.0a3
edgeDetect.cpp
Go to the documentation of this file.
1 //==============================================
2 // copyright : (C) 2003-2005 by Will Stokes
3 //==============================================
4 // This program is free software; you can redistribute it
5 // and/or modify it under the terms of the GNU General
6 // Public License as published by the Free Software
7 // Foundation; either version 2 of the License, or
8 // (at your option) any later version.
9 //==============================================
10 
11 //Systemwide includes
12 #include <qpainter.h>
13 #include <qimage.h>
14 #include <math.h>
15 
16 //Projectwide includes
17 #include "edgeDetect.h"
18 #include "blur.h"
19 #include "../enhancements/contrast.h"
20 
21 #define MIN(x,y) ((x) < (y) ? (x) : (y))
22 #define MAX(x,y) ((x) < (y) ? (x) : (y))
23 
24 //----------------------------------------------
25 // Inputs:
26 // -------
27 // QImage* image - image to find edges in
28 //
29 // Outputs:
30 // --------
31 // QImage* getEdgeImage - returns the produced grayscale edge image
32 //
33 // Other information such as pixel groups etc is also availble through
34 // various accesor method.
35 //
36 // Description:
37 // ------------
38 // This class is the first known publically available implementation of
39 // Kim, Lee, and Kweon's 2003 paper titled:
40 // "Automatic edge detection using 3x3 ideal binary
41 // pixel patterns and fuzzy-based edge thresholding"
42 // http://rcv.kaist.ac.kr/publication/file/foreign_journal/28_DongSuKim_PRL2004.pdf
43 //
44 // Edge detection is an old problem, an while many use the edge detectors by
45 // Canny, Sobel, etc, they all suffer from a common problem: the user must
46 // tweak a series of poorly understood input parameters to get the ideal edge image.
47 //
48 // Album Shaper was in need of an automatic edge detector for use when
49 // blurring and sharpening images. Having the user manually tweak such
50 // paramters first would not only be annoying but error-prone. In an
51 // effort to make edge detection automatic I took a stab at implementing
52 // this paper and am quite happy with the resulsts...
53 //
54 // http://albumshaper.sourceforge.net/images/teasers/peaksAndValleys.jpg
55 //
56 // Algorithm:
57 // ----------
58 // While complex, the algorithm can be broken up into a series of
59 // fairly straightforward tasks:
60 //
61 // 1.) "allocateAndInitObjects()" is called to allocate and fill a
62 // few data structures that will be used when finding image edges.
63 //
64 // 2.) "fillLumMapAndLumHistogram"() is called, during which an
65 // luminance map is constructed and luminance histogram populated.
66 // For an m x n image, the luminance map will be a m x n integer array.
67 //
68 // 3.) The luminance histogram is smoothed using "smoothLumHistogram" to
69 // make peak finding less sensative to noise.
70 //
71 // 4.) The fourth step is a little complicated. The edge magnitude and
72 // GSLC (Grey level similitude code) value is computed at each pixel.
73 // The paper takes an interesting approach to edge detection by
74 // classifying pixels into one of 9 groups by first computing the
75 // average luminance for a 3x3 group centered about a pixel. Pixels are
76 // then separated into one of two groups, those that have a luminance
77 // greater than or less than the 3x3 average luminance. For example:
78 //
79 // X
80 // --------- --------- ---------X
81 // | 7 15 18 | | 0 1 1 | | 0 1 X |
82 // | 5 17 20 | --> | 0 1 1 | --> | 0 X 1 |
83 // | 9 8 3 | | 0 0 0 | | X 0 0 |
84 // --------- --------- X---------
85 // X
86 //
87 // Here the average luminance is 11.333, placing the top right 4
88 // pixels in one group and the other remaining pixels in another.
89 // The dominant edge diretion is from the bottom left the top right.
90 // The GLSC code is computing by considering the 1/0 values
91 // (1 = pixel in same group as central pixel). The central value
92 // is ignored (it's always 1) leaving us with an 8bit = 2^8 = 256 code.
93 //
94 // In this case the GSLC is:
95 // 0*2^0 + 1*2^1 + 1*2^2 +
96 // 0*2^3 + XXXXX + 1*2^4 +
97 // 0*2^5 + 0*2^6 + 0*2^7 = 22.
98 //
99 // The authors of the paper found pixels with one of five GSLC
100 // codes (15,31,7,47,11) were most often associated with edge pixels
101 // when producing an edge image using various competitive techinques.
102 // By looking up the GSLC for a given pixel later one we can suppress
103 // edges where they most likely do not belong.
104 //
105 // 5.) The fifth step involves grouping pixels by luminance using the
106 // smoothed luminance histogram. This complicated step is brushed off
107 // as being trivial in the paper. Since I struggled a bit with developing
108 // an algorithm for this step, I'll explain my approach in detail
109 // to avoid others suffering.
110 //
111 // Using the smoothed histogram, we first compute the JND or just
112 // noticible differnce using the maximum count that was found. I'm
113 // not sure how appropriate this is, but 2% is the usual quantity
114 // used in other contexts, and it works well here, so 2% it is.
115 //
116 // Next we walk through the smoothed luminance histogram and find
117 // the midpoint of the valleys. We accomplish this by updating an
118 // index of the deepend last known valley. As the valley slopes
119 // down and we move across it we update this last best known index.
120 // Once the valley slopes up one JND above the last deepend location
121 // found we mark that valley midpoint and move on.
122 //
123 // Once all valley midpoints have been marked we can quickly deduce
124 // how many peaks must exist. We pass across the smoothed luminance
125 // histogram again finding the peak index for each pixel between
126 // valley midpoints. For all future work pixels one JND+- the peak
127 // center will be used.
128 //
129 // 6.) The sixth step, "computeClusterStatistics()", computes various
130 // cluster-specific statsitics that will be used to determine cluster
131 // thresholds. The paper was rather vague in this area, but after
132 // experimenting with various interpreatations of what they were trying
133 // to say I think I finally got it right.
134 //
135 // First, we iterate over all image pixels, determine which pixel group
136 // they belong by comparing luminance endpoints for all pixelclusters,
137 // and update total edge mag, num pixels, and an edgeMagHistogram for
138 // the given pixel group they belong to.
139 //
140 // Next we compute the average edge meganitude and most frequent edge
141 // magnitude observed for each pixel cluster, in addition to normalizing
142 // the cluster pixel count variable to [0,1]
143 //
144 // 7.) The seventh step is quite complicated and encompases computing the
145 // edge thresholds for each pixel cluster using the 18-rule fuzzy
146 // logic approach put forward by the paper. There is nothing ground
147 // breaking here, just a lot of complicated fuzzy logic, although most
148 // of the effort is put into computing the centroid at the end. I had
149 // never touched fuzzy logic before, but found this article more than
150 // helpful getting myself up to speed:
151 //
152 // http://www.doc.ic.ac.uk/~nd/surprise_96/journal/vol2/sbaa/article2.html
153 //
154 // 8.) The eigth and final step, actually constructing the edge image, is fairly
155 // straight forward and employs techinques (such as non-maximum suppression or NMS)
156 // anyone who has implemented Canny shoudl be familiar with.
157 //
158 // First, one looks up the ESF (edge shape factor) for a given pixel from a look-up
159 // table. These values were computed by generating a ton of edge images by carefully
160 // setting Canny, Sobel, etc params, then identifying how often a pixel
161 // with a given GSLC code was judged to be an edge. Hence ESF's will fall in
162 // the [0,1] range and help to suppress edges where no clear edge really
163 // can be claimed to exist, e.g.
164 //
165 // . # .
166 // # # # <- If there is an edge here, care to explain what the
167 // . # . edge direction is?!!
168 //
169 //
170 // If the ESF for a give pixel is 0 we skip it, it's not an edge.
171 //
172 // Next, we look up a pixels edge magnitude threshold by identifying the
173 // pixel cluster is belongs to using the pixels luminance using the
174 // luminance map.
175 //
176 // If the pixels edge magnitude is less than the threshold we skip the pixel,
177 // this filtersout low lying noise.
178 //
179 // Finally, the direction of the pixel is looked up using its GSLC and
180 // NMS (non-maximum suppression is applied). If the pixel has a greater
181 // egdge magnitude than either of its neighbors along the edge direction then
182 // the edge is marked.
183 //
184 // Final Remarks:
185 // --------------
186 // Despite the involved complexity, the implementation appears to work
187 // really really well. I consider this one of the secret gems of Album Shaper
188 // and hope to make good use of it in the future for things other than
189 // just sharpening and bluring. The only caveat is that edge detection
190 // does take a few CPU cycles.
191 //----------------------------------------------
192 
193 
194 //==============================================
195 EdgeDetect::EdgeDetect( QImage* image )
196 {
197  //load image
198  this->image = image;
199 
200  //allocate and initialize objects used for edge detection
202 
203  //fill lum map and lum histogram
205 
206  //fill smoothed histogram
208 
209  //compute edge magnitude and GSLC maps
211 
212  //determine pixel clusters
214 
216 
218 
220 }
221 //==============================================
223 {
225 }
226 //==============================================
228 { return numClusters; }
229 //==============================================
231 { return clusters; }
232 //==============================================
234 { return clusterPeaks; }
235 //==============================================
237 { return smoothLumHist; }
238 //==============================================
240 {
241  return image;
242 }
243 //==============================================
245 {
246  //construct map
247  int* clusterMap = new int[image->width() * image->height()];
248 
249  //iterate over all pixels, determine cluster each pixel belongs to
250  int i, cluster;
251  for(i=0; i<image->width()*image->height(); i++)
252  {
253  for(cluster=0; cluster<numClusters; cluster++)
254  {
255  if( lumMap[i] >= clusters[cluster].minLuminance &&
256  lumMap[i] <= clusters[cluster].maxLuminance )
257  {
258  clusterMap[i] = cluster;
259  break;
260  }
261  } //cluster
262  } //pixel
263 
264  return clusterMap;
265 }
266 //==============================================
268 {
269  //initialize:
270  //-luminosity histogram
271  //-smoothed luminosity histogram
272  //-identified peak regions
273  int i;
274  for(i=0; i<256; i++)
275  {
276  lumHist[i] = 0;
277  smoothLumHist[i] = 0;
278  clusterPeaks[i] = 0;
279  }
280 
281  //allocate luminance map
282  lumMap = new int[image->width() * image->height()];
283 
284  //allocate edge magnitude map
285  edgeMagMap = new float[image->width() * image->height()];
286 
287  //allocate GSLC map
288  GSLCmap = new int[image->width() * image->height()];
289 
290  //construct LUT
292 }
293 //==============================================
295 {
296  int x, y;
297  QRgb* rgb;
298  uchar* scanLine;
299  int lumVal;
300  for( y=0; y<image->height(); y++)
301  {
302  scanLine = image->scanLine(y);
303  for( x=0; x<image->width(); x++)
304  {
305  //get lum value for this pixel
306  rgb = ((QRgb*)scanLine+x);
307  lumVal = qGray(*rgb);
308 
309  //store in lum map
310  lumMap[x + y*image->width()] = lumVal;
311 
312  //update lum histogram
313  lumHist[ lumVal ]++;
314  }
315  }
316 }
317 //==============================================
319 {
320  #define FILTER_SIZE 5
321  int filter[FILTER_SIZE] = {2, 5, 8, 5, 2};
322 
323  int i,j;
324  int filterIndex, sum, total;
325  for(i = 0; i<256; i++)
326  {
327  sum = 0;
328  total = 0;
329 
330  for( j= -FILTER_SIZE/2; j <= FILTER_SIZE/2; j++)
331  {
332  if( i+j > 0 && i+j < 256 )
333  {
334  filterIndex = j+ FILTER_SIZE/2;
335  total+= filter[filterIndex] * lumHist[i+j];
336  sum += filter[filterIndex];
337  }
338  }
339 
340  smoothLumHist[i] = total / sum;
341  }
342 }
343 //==============================================
345 {
346  int x, y;
347  int idealPattern[9];
348  int pixelLums[9];
349 
350  //-------
351  //iterate over all pixels
352  for( y=0; y<image->height(); y++)
353  {
354  for( x=0; x<image->width(); x++)
355  {
356  //compute pixel luminances for entire grid
357  pixelLums[0] = pixelLum(x-1,y-1);
358  pixelLums[1] = pixelLum(x ,y-1);
359  pixelLums[2] = pixelLum(x+1,y-1);
360  pixelLums[3] = pixelLum(x-1,y );
361  pixelLums[4] = pixelLum(x ,y );
362  pixelLums[5] = pixelLum(x+1,y );
363  pixelLums[6] = pixelLum(x-1,y+1);
364  pixelLums[7] = pixelLum(x ,y+1);
365  pixelLums[8] = pixelLum(x+1,y+1);
366 
367  //compute average
368  float avg = 0;
369  int i;
370  for(i=0; i<=8; i++)
371  {
372  avg+= pixelLums[i];
373  }
374  avg = avg / 9;
375 
376  //determine ideal pattern and I0 and I1 averages
377  int centerPixelLum = pixelLums[4];
378  float centerDiff = centerPixelLum - avg;
379 
380  float I0avg = 0;
381  int I0count = 0;
382 
383  float I1avg = 0;
384  int I1count = 0;
385 
386  for(i=0; i<=8; i++)
387  {
388  if( centerDiff * (pixelLums[i]-avg) >=0 )
389  {
390  I1avg+=pixelLums[i];
391  I1count++;
392  idealPattern[i] = 1;
393  }
394  else
395  {
396  I0avg+=pixelLums[i];
397  I0count++;
398  idealPattern[i] = 0;
399  }
400  }
401 
402  //compute and store edge magnitude
403  if(I0count > 0) I0avg = I0avg/I0count;
404  if(I1count > 0) I1avg = I1avg/I1count;
405  edgeMagMap[x + y*image->width()] = QABS( I1avg - I0avg );
406 
407  //compute and store GSLC
408  int GSLC=0;
409  int weight = 1;
410  for(i=0; i<9; i++)
411  {
412  //skip center
413  if(i == 4) continue;
414 
415  if(idealPattern[i] == 1)
416  { GSLC+=weight; }
417 
418  weight = weight*2;
419  }
420  GSLCmap[x + y*image->width()] = GSLC;
421  } //x
422  } //y
423 }
424 //==============================================
425 int EdgeDetect::pixelLum(int x, int y)
426 {
427  int clampedX = MAX( MIN( x, image->width()-1), 0);
428  int clampedY = MAX( MIN( y, image->height()-1), 0);
429  return lumMap[ clampedX + clampedY * image->width() ];
430 }
431 //==============================================
433 {
434  //find max count
435  int maxCount = 0;
436  int i;
437  for(i=0; i<256; i++)
438  {
439  if(smoothLumHist[i] > maxCount)
440  maxCount = smoothLumHist[i];
441  }
442 
443  //compute JND for histogram (2% of total spread)
444  int histJND = maxCount/50;
445 
446  //construct temporary array for valley locations
447  //1's will indicate a valley midpoint
448  int tmpValleyArray[256];
449  for(i=0; i<256; i++) { tmpValleyArray[i] = 0; }
450 
451  //move across histogram finding valley midpoints
452  int curTrackedMin = smoothLumHist[0];
453 
454  //first and last indices tracked min was observed
455  int firstMinIndex = 0;
456  int lastMinIndex = 0;
457 
458  //only add valley midpoint if finished tracking a descent
459  bool slopeNeg = false;
460 
461  for(i = 1; i<256; i++ )
462  {
463  if( smoothLumHist[i] < curTrackedMin - histJND )
464  {
465  //found a descent!
466  slopeNeg = true;
467  curTrackedMin = smoothLumHist[i];
468  firstMinIndex = i;
469  }
470  //starting to go up again, add last min to list
471  else if( smoothLumHist[i] > curTrackedMin + histJND )
472  {
473  //if finished tracing a negative slope find midpoint and set location to true
474  if(slopeNeg)
475  {
476  tmpValleyArray[ (firstMinIndex + lastMinIndex)/2 ] = 1;
477  }
478 
479  curTrackedMin = smoothLumHist[i];
480  slopeNeg = false;
481  }
482  else
483  {
484  //still tracking a min, update the right
485  //hand index. center of valley is found
486  //by averaging first and last min index
487  lastMinIndex = i;
488  }
489  }
490 
491  //count valleys
492  int numValleys = 0;
493  for(i=0; i<256; i++)
494  {
495  if(tmpValleyArray[i] == 1 ) numValleys++;
496  }
497 
498  //determine number of clusters
499  numClusters = numValleys-1;
500  if(tmpValleyArray[0] != 1)
501  numClusters++;
502  if(tmpValleyArray[255] != 1)
503  numClusters++;
504 
505  //allocate clusters
507 
508  //automatically start first cluster
509  int cluster=0;
510  clusters[cluster].minLuminance = 0;
511 
512  //initialize left and right boundaries of all clusters
513  for(i=1; i<256; i++)
514  {
515  //reached next valley, end cluster
516  if( tmpValleyArray[i] == 1)
517  {
518  clusters[cluster].maxLuminance = i-1;
519  cluster++;
520  clusters[cluster].minLuminance = i;
521  }
522  //end last cluster automatically at end
523  else if(i == 255)
524  {
525  clusters[cluster].maxLuminance = i;
526  }
527  }
528 
529  //determine cluster peaks
530  for(cluster=0; cluster<numClusters; cluster++)
531  {
532  //find max for current cluster
533  int maxIndex = clusters[cluster].minLuminance;
534  for(i=clusters[cluster].minLuminance; i<=clusters[cluster].maxLuminance; i++)
535  {
536  if(smoothLumHist[i] > smoothLumHist[maxIndex])
537  maxIndex = i;
538  }
539 
540  //mark peaks
541  int lumJND = 255/50;
542  for(i=MAX(0, maxIndex-lumJND); i<MIN(256, maxIndex+lumJND); i++)
543  {
544  clusterPeaks[i] = 1;
545  }
546  }
547 }
548 //==============================================
550 {
551  //initialize cluster stats
552  int cluster;
553  for(cluster=0; cluster<numClusters; cluster++)
554  {
555  int i;
556  for(i=0; i<256; i++)
557  {
558  clusters[cluster].edgeMagHistogram[i] = 0;
559  }
560  clusters[cluster].totalEdgeMagnitude=0.0f;
561  clusters[cluster].numPixels = 0;
562  }
563 
564  //iterate over all pixels
565  int i;
566  for(i=0; i<image->width()*image->height(); i++)
567  {
568  //skip pixels that don't belong to peaks
569  if( clusterPeaks[ lumMap[i] ] != 1)
570  continue;
571 
572  //determine cluster pixel belongs to
573  int cluster;
574  for(cluster=0; cluster<numClusters; cluster++)
575  {
576  if( lumMap[i] >= clusters[cluster].minLuminance &&
577  lumMap[i] <= clusters[cluster].maxLuminance )
578  {
579  clusters[cluster].totalEdgeMagnitude+= edgeMagMap[i];
580  clusters[cluster].numPixels++;
581  clusters[cluster].edgeMagHistogram[ MIN( MAX( (int)edgeMagMap[i], 0), 255) ]++;
582  break;
583  }
584  } //cluster
585  } //pixel i
586 
587  //iterate over clusters to determine min and max peak cluster sizes
590  for(cluster=1; cluster<numClusters; cluster++)
591  {
592  if(clusters[cluster].numPixels < minClusterSize)
593  minClusterSize = clusters[cluster].numPixels;
594 
595  if(clusters[cluster].numPixels > maxClusterSize)
596  maxClusterSize = clusters[cluster].numPixels;
597  }
598 
599  //iterate over clusters one final time to deduce normalized inputs to fuzzy logic process
600  int JND = 255/50;
601  for(cluster=0; cluster<numClusters; cluster++)
602  {
603  clusters[cluster].meanMode = MIN( clusters[cluster].totalEdgeMagnitude / clusters[cluster].numPixels,
604  3*JND );
605 
606  int i;
607  int mode = 0;
608  for(i=1; i<256; i++)
609  {
610  if( clusters[cluster].edgeMagHistogram[i] > clusters[cluster].edgeMagHistogram[ mode ] )
611  mode = i;
612  }
613  clusters[cluster].mode = MIN( mode, 2*JND );
614 
615  clusters[cluster].pixelCount = ((float)(clusters[cluster].numPixels - minClusterSize)) /
617  }
618 }
619 //==============================================
620 //compute edge thresholds for each cluster using 18-rule fuzzy logic approach
622 {
623  //iterate over each cluster
624  int cluster;
625  float S1,M1,L1;
626  float S2,M2,L2;
627  float S3,L3;
628  float outS, outM, outL;
629 
630  int JND = 255/50;
631 
632  for(cluster=0; cluster<numClusters; cluster++)
633  {
634  //----
635  //compute S,M, and L values for each input
636  //----
637  S1 = MAX( 1.0f - ((clusters[cluster].meanMode/JND) / 1.5f), 0 );
638 
639  if( (clusters[cluster].meanMode/JND) <= 1.5f )
640  M1 = MAX( (clusters[cluster].meanMode/JND) - 0.5f, 0 );
641  else
642  M1 = MAX( 2.5f - (clusters[cluster].meanMode/JND), 0 );
643 
644  L1 = MAX( ((clusters[cluster].meanMode/JND) - 1.5f) / 1.5f, 0 );
645  //----
646  S2 = MAX( 1.0f - (clusters[cluster].mode/JND), 0 );
647 
648  if( (clusters[cluster].mode/JND) <= 1.0f )
649  M2 = MAX( -1.0f + 2*(clusters[cluster].mode/JND), 0 );
650  else
651  M2 = MAX( 3.0f - 2*(clusters[cluster].mode/JND), 0 );
652 
653  L2 = MAX( (clusters[cluster].mode/JND) - 1.0, 0 );
654  //----
655  S3 = MAX( 1.0f - 2*clusters[cluster].pixelCount, 0 );
656  L3 = MAX( -1.0f + 2*clusters[cluster].pixelCount, 0 );
657  //----
658 
659  //Compute M,L for outputs using set of 18 rules.
660  //outS is inherantly S given the ruleset provided
661  outS = 0.0f;
662  outM = 0.0f;
663  outL = 0.0f;
664  //Out 1
665  if( numClusters > 2 )
666  {
667  outM += S1*S2*S3; //rule 1
668 
669  //rule 2
670  if( clusters[cluster].meanMode < clusters[cluster].mode )
671  outS += S1*S2*L3;
672  else
673  outM += S1*S2*L3;
674 
675  outM += S1*M2*S3; //rule 3
676  outM += S1*M2*L3; //rule 4
677  outM += S1*L2*S3; //rule 5
678  outM += S1*L2*L3; //rule 6
679  outM += M1*S2*S3; //rule 7
680  outM += M1*S2*L3; //rule 8
681  outM += M1*M2*S3; //rule 9
682  outL += M1*M2*L3; //rule 10
683  outM += M1*L2*S3; //rule 11
684  outL += M1*L2*L3; //rule 12
685  outM += L1*S2*S3; //rule 13
686  outL += L1*S2*L3; //rule 14
687  outM += L1*M2*S3; //rule 15
688  outL += L1*M2*L3; //rule 16
689  outL += L1*L2*S3; //rule 17
690  outL += L1*L2*L3; //rule 18
691  }
692  //Out 2
693  else
694  {
695  outL += S1*S2*S3; //rule 1
696  outL += S1*S2*L3; //rule 2
697  outM += S1*M2*S3; //rule 3
698  outL += S1*M2*L3; //rule 4
699  outM += S1*L2*S3; //rule 5
700  outM += S1*L2*L3; //rule 6
701  outL += M1*S2*S3; //rule 7
702  outL += M1*S2*L3; //rule 8
703  outL += M1*M2*S3; //rule 9
704  outL += M1*M2*L3; //rule 10
705  outL += M1*L2*S3; //rule 11
706  outL += M1*L2*L3; //rule 12
707  outL += L1*S2*S3; //rule 13
708  outL += L1*S2*L3; //rule 14
709  outL += L1*M2*S3; //rule 15
710  outL += L1*M2*L3; //rule 16
711  outL += L1*L2*S3; //rule 17
712  outL += L1*L2*L3; //rule 18
713  }
714 
715  //find centroid - Beta[k]
716  float A = outM + 0.5f;
717  float B = 2.5f - outM;
718  float C = 1.5f * (outL + 1);
719  float D = 1.5f * (outM + 1);
720  float E = 2.5f - outL;
721 
722  //---------------------------------------------------------------
723  //Case 1: Both outM and outL are above intersection point of diagonals
724  if( outM > 0.5f && outL > 0.5f )
725  {
726  //find area of 7 subregions
727  float area1 = ((A-0.5f)*outM)/2;
728  float area2 = outM * (B-A);
729  float area3 = ((2.1f-B) * (outM - 0.5)) / 2;
730  float area4 = (2.1 - B) * 0.5f;
731  float area5 = ((C - 2.1f) * (outL - 0.5)) / 2;
732  float area6 = (C - 2.1f) * 0.5f;
733  float area7 = (3.0f - C) * outL;
734 
735  //find half of total area
736  float halfArea = (area1 + area2 + area3 + area4 + area5 + area6 + area7) / 2;
737 
738  //determine which region split will be within and resulting horizontal midpoint
739 
740  //Within area 1
741  if( area1 > halfArea )
742  {
743  clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
744  }
745  //Within area 2
746  else if( area1 + area2 > halfArea )
747  {
748  clusters[cluster].beta = ((halfArea - area1) / outM) + A;
749  }
750  //Within area 3-4
751  else if( area1 + area2 + area3 + area4 > halfArea )
752  {
753  float a = -0.5f;
754  float b = 2.8f;
755  float c = area1 + area2 + area3 - halfArea - B/2 - 2.625f;
756  clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
757  }
758  //Within area 5-6
759  else if( area1 + area2 + area3 + area4 + area5 + area6 > halfArea )
760  {
761  float a = 1.0f/3;
762  float b = -0.7f;
763  float c = area1 + area2 + area3 + area4 - halfArea;
764  clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
765  }
766  //Within area 7
767  else
768  {
769  clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4 + area5 + area6) ) / outL) + C;
770  }
771  } //end case 1
772  //---------------------------------------------------------------
773  //Case 2
774  else if ( outM < 0.5f && outL > outM )
775  {
776  //find area of 5 subregions
777  float area1 = (outM*(A-0.5f)) / 2;
778  float area2 = (D-A) * outM;
779  float area3 = ((C-D) * (outL - outM)) / 2;
780  float area4 = (C-D) * outM;
781  float area5 = (3.0f - C) * outL;
782 
783  //find half of total area
784  float halfArea = (area1 + area2 + area3 + area4 + area5) / 2;
785 
786  //determine which region split will be within and resulting horizontal midpoint
787 
788  //Within area 1
789  if( area1 > halfArea )
790  {
791  clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
792  }
793  //Within area 2
794  else if( area1 + area2 > halfArea )
795  {
796  clusters[cluster].beta = ((halfArea - area1) / outM) + A;
797  }
798  //Within area 3-4
799  else if( area1 + area2 + area3 + area4 > halfArea )
800  {
801  float a = 1.0f/3.0f;
802  float b = outM - 0.5f - D/3;
803  float c = area1 + area2 - D*outM + D/2 - halfArea;
804  clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
805  }
806  //Within area5
807  else
808  {
809  clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4) ) / outL) + C;
810  }
811  } //end case 2
812  //---------------------------------------------------------------
813  //Case 3
814  else
815  {
816  //find area of 5 subregions
817  float area1 = (outM*(A-0.5f)) / 2;
818  float area2 = (B-A) * outM;
819  float area3 = ((E-B) * (outM - outL)) / 2;
820  float area4 = (E-B) * outL;
821  float area5 = (3.0f - E) * outL;
822 
823  //find half of total area
824  float halfArea = (area1 + area2 + area3 + area4 + area5) / 2;
825 
826  //determine which region split will be within and resulting horizontal midpoint
827 
828  //Within area 1
829  if( area1 > halfArea )
830  {
831  clusters[cluster].beta = 0.5f + (float)sqrt(2*halfArea);
832  }
833  //Within area 2
834  else if( area1 + area2 > halfArea )
835  {
836  clusters[cluster].beta = ((halfArea - area1) / outM) + A;
837  }
838  //Within area 3-4
839  else if( area1 + area2 + area3 + area4 > halfArea )
840  {
841  float a = -0.5f;
842  float b = E/2 + 2.5f/2;
843  float c = area3 - 2.5f*E/2;
844  clusters[cluster].beta = (-b + (float)sqrt( b*b - 4*a*c )) / (2*a);
845  }
846  //Within area5
847  else
848  {
849  clusters[cluster].beta = ((halfArea - (area1 + area2 + area3 + area4) ) / outL) + E;
850  }
851  } //end case 3
852  //---------------------------------------------------------------
853 
854  //Compute edge threshold
855  int lumJND = 255/50;
856  clusters[cluster].edgeThreshold = clusters[cluster].mode + clusters[cluster].beta*lumJND;
857 
858  } //end for cluster
859 
860 }
861 //==============================================
863 {
864  int x, y;
865  QRgb* rgb;
866 
867  uchar* scanLine;
868  for( y=0; y<image->height(); y++)
869  {
870  scanLine = image->scanLine(y);
871  for( x=0; x<image->width(); x++)
872  {
873  //initialize pixel to black
874  rgb = ((QRgb*)scanLine+x);
875  *rgb = qRgb( 0, 0, 0 );
876 
877  //lookup ESF for this pixel
878  float ESF = LUT[ GSLCmap[x + y*image->width()] ].ESF;
879 
880  //If ESF value for this pixel is 0 skip
881  if( ESF == 0.0f ) continue;
882 
883  //lookup edge magnitude threshold
884  float lum = lumMap[x + y*image->width()];
885  float edgeMagThresh = -1.0f;
886  int cluster;
887  for(cluster=0; cluster<numClusters; cluster++)
888  {
889  if(lum >= clusters[cluster].minLuminance &&
890  lum <= clusters[cluster].maxLuminance)
891  {
892  edgeMagThresh = clusters[cluster].edgeThreshold;
893  break;
894  }
895  }
896 
897  //if cluster not found bail
898  if( cluster >= numClusters )
899  {
900 // cout << "Error! Could not find cluster pixel belonged to!\n";
901  continue;
902  }
903 
904  //if edge mag below thresh then skip
905  if( edgeMagMap[x + y*image->width()] < edgeMagThresh ) continue;
906 
907  //ok, last checks implement NMS (non-maximum supression)
908  int direction = LUT[ GSLCmap[x + y*image->width()] ].direction;
909  int neighborIndex1 = -1;
910  int neighborIndex2 = -1;
911 
912  if( direction == 0)
913  {
914  if( x > 0)
915  neighborIndex1 = x-1 + y*image->width();
916  if( x < image->width() - 1 )
917  neighborIndex2 = x+1 + y*image->width();
918  }
919  else if(direction == 1)
920  {
921  if( x > 0 && y < image->height() - 1 )
922  neighborIndex1 = x-1 + (y+1)*image->width();
923  if( x < image->width() - 1 && y > 0 )
924  neighborIndex2 = x+1 + (y-1)*image->width();
925  }
926  else if(direction == 2)
927  {
928  if( y < image->height() - 1 )
929  neighborIndex1 = x + (y+1)*image->width();
930  if( y > 0)
931  neighborIndex2 = x + (y-1)*image->width();
932  }
933  else if(direction == 3)
934  {
935  if( x < image->width() - 1 && y < image->height() - 1 )
936  neighborIndex1 = x+1 + (y+1)*image->width();
937  if( x > 0 && y > 0 )
938  neighborIndex2 = x-1 + (y-1)*image->width();
939  }
940 
941  //neighbor 1 has higher confidence, skip!
942  if( neighborIndex1 != -1 &&
943  LUT[ GSLCmap[neighborIndex1] ].ESF * edgeMagMap[neighborIndex1] >
944  ESF * edgeMagMap[x + y*image->width()] )
945  continue;
946 
947  //neighbor 2 has higher confidence, skip!
948  if( neighborIndex2 != -1 &&
949  LUT[ GSLCmap[neighborIndex2] ].ESF * edgeMagMap[neighborIndex2] >
950  ESF * edgeMagMap[x + y*image->width()] )
951  continue;
952 
953  //All tests passed! Mark edge!
954  *rgb = qRgb( 255, 255, 255 );
955  } //x
956  } //y
957 
958  //blur image - all of it
959  blurImage( *image, 2.0f );
960 
961  //normalize image
963 
964 }
965 //==============================================
967 {
968  delete[] lumMap; lumMap = NULL;
969  delete[] edgeMagMap; edgeMagMap = NULL;
970  delete[] GSLCmap; GSLCmap = NULL;
971  delete[] clusters; clusters = NULL;
972 }
973 //==============================================
975 {
976  //----------------------
977  //First fill entire table with 0 ESF's and invalid directions
978  int i;
979  for(i=0; i<256; i++)
980  {
981  LUT[i].ESF = 0.0f;
982  LUT[i].direction = -1;
983  }
984  //----------------------
985  //Next code in all pattern that are highly
986  //likely to be on edges as described in the paper
987  //----------------------
988  //Pattern (a)
989 
990  // ###
991  // ##.
992  // ...
993  LUT[15].ESF = 0.179f;
994  LUT[15].direction = 3;
995 
996  // ...
997  // .##
998  // ###
999  LUT[240].ESF = 0.179f;
1000  LUT[240].direction = 3;
1001 
1002  // ###
1003  // .##
1004  // ...
1005  LUT[23].ESF = 0.179f;
1006  LUT[23].direction = 1;
1007 
1008  // ...
1009  // ##.
1010  // ###
1011  LUT[232].ESF = 0.179f;
1012  LUT[232].direction = 1;
1013 
1014  // ##.
1015  // ##.
1016  // #..
1017  LUT[43].ESF = 0.179f;
1018  LUT[43].direction = 3;
1019 
1020  // ..#
1021  // .##
1022  // .##
1023  LUT[212].ESF = 0.179f;
1024  LUT[212].direction = 3;
1025 
1026  // #..
1027  // ##.
1028  // ##.
1029  LUT[105].ESF = 0.179f;
1030  LUT[105].direction = 1;
1031 
1032  // .##
1033  // .##
1034  // ..#
1035  LUT[150].ESF = 0.179f;
1036  LUT[150].direction = 1;
1037  //----------------------
1038  //Pattern (b)
1039 
1040  // ###
1041  // ###
1042  // ...
1043  LUT[31].ESF = 0.137f;
1044  LUT[31].direction = 2;
1045 
1046  // ...
1047  // ###
1048  // ###
1049  LUT[248].ESF = 0.137f;
1050  LUT[248].direction = 2;
1051 
1052  // ##.
1053  // ##.
1054  // ##.
1055  LUT[107].ESF = 0.137f;
1056  LUT[107].direction = 0;
1057 
1058  // .##
1059  // .##
1060  // .##
1061  LUT[214].ESF = 0.137f;
1062  LUT[214].direction = 0;
1063  //----------------------
1064  //Pattern (c)
1065 
1066  // ###
1067  // .#.
1068  // ...
1069  LUT[7].ESF = 0.126f;
1070  LUT[7].direction = 2;
1071 
1072  // ...
1073  // .#.
1074  // ###
1075  LUT[224].ESF = 0.126f;
1076  LUT[224].direction = 2;
1077 
1078  // #..
1079  // ##.
1080  // #..
1081  LUT[41].ESF = 0.126f;
1082  LUT[41].direction = 0;
1083 
1084  // ..#
1085  // .##
1086  // ..#
1087  LUT[148].ESF = 0.126f;
1088  LUT[148].direction = 0;
1089  //----------------------
1090  //Pattern (d)
1091 
1092  // ###
1093  // ##.
1094  // #..
1095  LUT[47].ESF = 0.10f;
1096  LUT[47].direction = 3;
1097 
1098  // ..#
1099  // .##
1100  // ###
1101  LUT[244].ESF = 0.10f;
1102  LUT[244].direction = 3;
1103 
1104  // ###
1105  // .##
1106  // ..#
1107  LUT[151].ESF = 0.10f;
1108  LUT[151].direction = 1;
1109 
1110  // #..
1111  // ##.
1112  // ###
1113  LUT[233].ESF = 0.10f;
1114  LUT[233].direction = 1;
1115  //----------------------
1116  //Pattern (e)
1117 
1118  // ##.
1119  // ##.
1120  // ...
1121  LUT[11].ESF = 0.10f;
1122  LUT[11].direction = 3;
1123 
1124  // ...
1125  // .##
1126  // .##
1127  LUT[208].ESF = 0.10f;
1128  LUT[208].direction = 3;
1129 
1130  // .##
1131  // .##
1132  // ...
1133  LUT[22].ESF = 0.10f;
1134  LUT[22].direction = 1;
1135 
1136  // ...
1137  // ##.
1138  // ##.
1139  LUT[104].ESF = 0.10f;
1140  LUT[104].direction = 1;
1141  //----------------------
1142 }
1143 //==============================================
int minClusterSize
Definition: edgeDetect.h:116
int * getClusterMap()
Definition: edgeDetect.cpp:244
#define MAX(x, y)
Definition: edgeDetect.cpp:22
int smoothLumHist[256]
Definition: edgeDetect.h:97
float pixelCount
Definition: edgeDetect.h:38
void blurImage(QImage &image, float sigma)
Definition: blur.cpp:94
#define FILTER_SIZE
int direction
Definition: edgeDetect.h:23
int maxLuminance
Definition: edgeDetect.h:29
int * getPeaks()
Definition: edgeDetect.cpp:233
float B
Definition: blur.cpp:78
int lumHist[256]
luminosity and smooth luminosity histograms
Definition: edgeDetect.h:96
PixelCluster * clusters
Definition: edgeDetect.h:113
long b
Definition: jpegInternal.h:125
int clusterPeaks[256]
Definition: edgeDetect.h:100
QImage * enhanceImageContrast(QString filename, StatusWidget *status)
Definition: contrast.cpp:88
void findPixelClusters()
Definition: edgeDetect.cpp:432
float ESF
Definition: edgeDetect.h:22
float * edgeMagMap
Definition: edgeDetect.h:106
void deallocateObjects()
Definition: edgeDetect.cpp:966
EdgeDetect(QImage *image)
Definition: edgeDetect.cpp:195
int * lumMap
Definition: edgeDetect.h:103
void fillLumMapAndLumHistogram()
Definition: edgeDetect.cpp:294
QImage * image
Definition: edgeDetect.h:93
int width
Definition: blur.cpp:79
int edgeMagHistogram[256]
Definition: edgeDetect.h:31
float beta
Definition: edgeDetect.h:40
PixelCluster * getClusters()
Definition: edgeDetect.cpp:230
int minLuminance
Definition: edgeDetect.h:28
void constructGSLClut()
Definition: edgeDetect.cpp:974
void allocateAndInitObjects()
Definition: edgeDetect.cpp:267
void computeClusterStatistics()
Definition: edgeDetect.cpp:549
int pixelLum(int x, int y)
Definition: edgeDetect.cpp:425
float mode
Definition: edgeDetect.h:37
QImage * getEdgeImage()
Definition: edgeDetect.cpp:239
LUTentry LUT[256]
Definition: edgeDetect.h:90
int * getSmoothHist()
Definition: edgeDetect.cpp:236
void smoothLumHistogram()
Definition: edgeDetect.cpp:318
int numClusters
Definition: edgeDetect.h:112
int maxClusterSize
Definition: edgeDetect.h:116
void computeClusterThresholds()
Definition: edgeDetect.cpp:621
int getNumClusters()
Definition: edgeDetect.cpp:227
float edgeThreshold
Definition: edgeDetect.h:41
float totalEdgeMagnitude
Definition: edgeDetect.h:32
float meanMode
Definition: edgeDetect.h:36
void constructEdgeImage()
Definition: edgeDetect.cpp:862
void computeEdgeMagAndGSLCmaps()
Definition: edgeDetect.cpp:344
#define MIN(x, y)
Definition: edgeDetect.cpp:21
int height
Definition: blur.cpp:79
int * GSLCmap
Definition: edgeDetect.h:109