AlbumShaper
1.0a3
|
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 //==============================================