University Of Massachusetts at Lowell
Robotics Lab
00001 /* --------------------------------------------------------------------------- 00002 Phission : 00003 Realtime Vision Processing System 00004 00005 Copyright (C) 2003-2006 Philip D.S. Thoren (pthoren@cs.uml.edu) 00006 University of Massachusetts at Lowell, 00007 Laboratory for Artificial Intelligence and Robotics 00008 00009 This file is part of Phission. 00010 00011 Phission is free software; you can redistribute it and/or modify 00012 it under the terms of the GNU Lesser General Public License as published by 00013 the Free Software Foundation; either version 2 of the License, or 00014 (at your option) any later version. 00015 00016 Phission is distributed in the hope that it will be useful, 00017 but WITHOUT ANY WARRANTY; without even the implied warranty of 00018 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00019 GNU Lesser General Public License for more details. 00020 00021 You should have received a copy of the GNU Lesser General Public License 00022 along with Phission; if not, write to the Free Software 00023 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00024 00025 ---------------------------------------------------------------------------*/ 00026 #ifdef HAVE_CONFIG_H 00027 #include <phissionconfig.h> 00028 #endif 00029 00030 #include <phStandard.h> 00031 00032 #include <canny_Filter.h> 00033 00034 #define UNROLL_LOOP 1 00035 00036 #define TIMESTUFF() 0 00037 00038 #include <phMath.h> 00039 #include <phError.h> 00040 #include <phMemory.h> 00041 #include <phPrint.h> 00042 00043 /* ---------------------------------------------------------------------- */ 00044 canny_Filter::canny_Filter( uint32_t lowThreshold, 00045 uint32_t highThreshold ) : 00046 phFilter("canny_Filter") 00047 00048 { 00049 this->m_format = ( phImageRGB24 | 00050 phImageRGBA32 | 00051 phImageGREY8 ); 00052 00053 this->m_gauss = NULL; 00054 this->m_edge = NULL; 00055 this->m_gaussSize = 0; 00056 this->m_edgeSize = 0; 00057 00058 this->m_useGaussian = 1; 00059 this->m_useSobelApprox = 1; 00060 00061 this->set(lowThreshold, highThreshold); 00062 } 00063 00064 /* ---------------------------------------------------------------------- */ 00065 canny_Filter::~canny_Filter() 00066 { 00067 phFree(this->m_gauss); 00068 phFree(this->m_edge); 00069 } 00070 00071 /* ------------------------------------------------------------------------ */ 00072 phFilter *canny_Filter::cloneFilter() 00073 { 00074 phFUNCTION("canny_Filter::cloneFilter") 00075 int locked = 0; 00076 canny_Filter *canny = new canny_Filter(); 00077 00078 phTHIS_LOOSE_LOCK(locked); 00079 00080 canny->set(this->m_highThreshold, this->m_lowThreshold); 00081 00082 phTHIS_LOOSE_UNLOCK(locked); 00083 00084 return (phFilter *)canny; 00085 } 00086 00087 /* ---------------------------------------------------------------------- */ 00088 int canny_Filter::set( uint32_t lowThreshold, 00089 uint32_t highThreshold ) 00090 { 00091 phFUNCTION("canny_Filter::set") 00092 int locked = 0; 00093 00094 phTHIS_LOOSE_LOCK(locked); 00095 00096 this->m_highThreshold = highThreshold; 00097 this->m_lowThreshold = lowThreshold; 00098 00099 phTHIS_LOOSE_UNLOCK(locked); 00100 00101 return phSUCCESS; 00102 } 00103 00104 /* ------------------------------------------------------------------------ */ 00105 void canny_Filter::enableGaussian( int enable ) 00106 { 00107 this->m_useGaussian = (enable ? 1 : 0); 00108 } 00109 00110 /* ------------------------------------------------------------------------ */ 00111 void canny_Filter::disableGaussian( int disable ) 00112 { 00113 this->m_useGaussian = (disable ? 0 : 1); 00114 } 00115 00116 /* ------------------------------------------------------------------------ */ 00117 void canny_Filter::enableSobelApprox( int enable ) 00118 { 00119 this->m_useSobelApprox = (enable ? 1 : 0); 00120 } 00121 00122 /* ------------------------------------------------------------------------ */ 00123 void canny_Filter::disableSobelApprox( int disable ) 00124 { 00125 this->m_useSobelApprox = (disable ? 0 : 1); 00126 } 00127 00128 /* ---------------------------------------------------------------------- */ 00129 int canny_Filter::filter() 00130 { 00131 phFUNCTION("canny_Filter::filter") 00132 00133 /* Filter variables */ 00134 uint32_t r, c = 0; 00135 uint32_t row, col, offset = 0; 00136 uint32_t inrow, incol = 0; 00137 uint32_t row_length = 0; 00138 uint8_t *addr = NULL; 00139 uint8_t *toaddr = NULL; 00140 uint8_t *fromaddr = NULL; 00141 double Y = 0.0; 00142 int I, J = 0; 00143 int sumX = 0; 00144 int sumY = 0; 00145 int SUM = 0; 00146 int leftPixel = 0; 00147 int rightPixel = 0; 00148 int P1, P2, P3, P4, P5, P6, P7, P8 = 0; 00149 float ORIENT = 0.0; 00150 int edgeDirection = 0; 00151 /* 5x5 Gaussian mask. Ref: http://www.cee.hw.ac.uk/hipr/html/gsmooth.html */ 00152 int GAUSS_MASK[5][5] = { 00153 { 2, 4, 5, 4, 2}, 00154 { 4, 9,12, 9, 4}, 00155 { 5,12,15,12, 5}, 00156 { 4, 9,12, 9, 4}, 00157 { 2, 4, 5, 4, 2} 00158 }; 00159 /* 3x3 GX Sobel mask. Ref: www.cee.hw.ac.uk/hipr/html/sobel.html */ 00160 int GX[3][3] = { 00161 {-1,0,1}, 00162 {-2,0,2}, 00163 {-1,0,1} 00164 }; 00165 /* 3x3 GY Sobel mask. Ref: www.cee.hw.ac.uk/hipr/html/sobel.html */ 00166 int GY[3][3] = { 00167 { 1, 2, 1}, 00168 { 0, 0, 0}, 00169 {-1,-2,-1} 00170 }; 00171 #if TIMESTUFF() 00172 /* timing structures */ 00173 #endif 00174 /* End filter variables */ 00175 00176 /* Begin Filter */ 00177 phDALLOC_RESIZE(this->m_gauss, 00178 this->m_gaussSize, 00179 width * height * depth, 00180 uint8_t ); 00181 phDALLOC_RESIZE(this->m_edge, 00182 this->m_edgeSize, 00183 width * height * depth, 00184 uint8_t ); 00185 00186 row_length = width * depth; 00187 00188 /* 00189 * Source of code for Canny edge detector: canny.c by Bill Green 00190 * 00191 * Code altered to be faster and easier to read. 00192 * The code is also relatively generic with only 00193 * 00194 * 6 inputs: 00195 * width, 00196 * height, 00197 * depth, 00198 * *data, 00199 * high threshold, 00200 * low threshold 00201 * 00202 * Program Steps: 00203 * 1.) 5x5 Gaussian mask for noise removal 00204 * 2.) Use Sobel masks to approximate gradient 00205 * 3.) Find orientation of gradient 00206 * 4.) Implement nonmaximum suppression to assign edges 00207 * 5.) Hysteresis thresholding (2 thresholds) 00208 */ 00209 00210 /* RGB to Grey scale */ 00211 /* Common Luminance equation: Y = 0.3R + 0.59G + 0.11B */ 00212 if (this->m_workspaceImage->getFormat() & (phImageRGB24 | phImageRGBA32)) 00213 { 00214 uint32_t limit = height * width; 00215 #if TIMESTUFF() 00216 /* phPROGRESS("Luminance Conversion\n"); */ 00217 /* time start */ 00218 #endif 00219 addr = Image; 00220 for( r = 0; r < limit; r++, addr += depth ) 00221 { 00222 Y = (0.3 * (*(addr)) ) + /* R */ 00223 (0.59 * (*(addr + 1))) + /* G */ 00224 (0.11 * (*(addr + 2))); /* B */ 00225 00226 if (Y > 255) Y = 255; 00227 phMemset((void *)addr, 00228 (uint8_t)Y, 00229 sizeof(uint8_t)*depth); 00230 } 00231 #if TIMESTUFF() 00232 /* time end */ 00233 #endif 00234 } 00235 00236 if (this->m_useGaussian) 00237 { 00238 /* Gaussian */ 00239 uint32_t r_limit = height - 2; 00240 uint32_t c_limit = width - 2; 00241 #if TIMESTUFF() 00242 /* phPROGRESS("Gaussian Blur\n"); */ 00243 /* time start */ 00244 #endif 00245 00246 /* Start 2 rows down and 2 pixels in */ 00247 toaddr = this->m_gauss + (2 * row_length) + (2 * depth); 00248 00249 /* toaddr += 4 * depth: skip border 2 pixels on end of row 00250 * and front of next row */ 00251 for( r = 2; r < r_limit; r++, toaddr += 4 * depth ) 00252 { 00253 /* toaddr += depth: goto the next pixel */ 00254 for( c = 2; c < c_limit; c++, toaddr += depth ) 00255 { 00256 SUM = 0; 00257 00258 fromaddr = Image + 00259 ((r - 2) * row_length) + 00260 ((c - 2) * depth); 00261 #if UNROLL_LOOP 00262 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 0 ][ 0 ])); 00263 fromaddr += depth; 00264 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 0 ][ 1 ])); 00265 fromaddr += depth; 00266 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 0 ][ 2 ])); 00267 fromaddr += depth; 00268 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 0 ][ 3 ])); 00269 fromaddr += depth; 00270 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 0 ][ 4 ])); 00271 fromaddr += row_length - (4 * depth); 00272 00273 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 1 ][ 0 ])); 00274 fromaddr += depth; 00275 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 1 ][ 1 ])); 00276 fromaddr += depth; 00277 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 1 ][ 2 ])); 00278 fromaddr += depth; 00279 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 1 ][ 3 ])); 00280 fromaddr += depth; 00281 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 1 ][ 4 ])); 00282 fromaddr += (row_length - (4 * depth)); 00283 00284 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 2 ][ 0 ])); 00285 fromaddr += depth; 00286 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 2 ][ 1 ])); 00287 fromaddr += depth; 00288 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 2 ][ 2 ])); 00289 fromaddr += depth; 00290 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 2 ][ 3 ])); 00291 fromaddr += depth; 00292 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 2 ][ 4 ])); 00293 fromaddr += (row_length - (4 * depth)); 00294 00295 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 3 ][ 0 ])); 00296 fromaddr += depth; 00297 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 3 ][ 1 ])); 00298 fromaddr += depth; 00299 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 3 ][ 2 ])); 00300 fromaddr += depth; 00301 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 3 ][ 3 ])); 00302 fromaddr += depth; 00303 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 3 ][ 4 ])); 00304 fromaddr += (row_length - (4 * depth)); 00305 00306 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 4 ][ 0 ])); 00307 fromaddr += depth; 00308 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 4 ][ 1 ])); 00309 fromaddr += depth; 00310 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 4 ][ 2 ])); 00311 fromaddr += depth; 00312 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 4 ][ 3 ])); 00313 fromaddr += depth; 00314 SUM += (int)((*fromaddr) * (GAUSS_MASK[ 4 ][ 4 ])); 00315 #else /* UNROLL_LOOP */ 00316 for( I = 0; I < 5; I++ ) 00317 { 00318 for( J = 0; J < 5; J++ ) 00319 { 00320 SUM += (int)((*fromaddr) * (GAUSS_MASK[I][J])); 00321 fromaddr += depth; 00322 } 00323 /* one row down - width of MASKing box */ 00324 fromaddr += row_length - (5 * depth); 00325 } 00326 #endif /* UNROLL_LOOP */ 00327 00328 SUM = SUM / 115; 00329 if (SUM > 255) SUM = 255; 00330 00331 phMemset((void *)(toaddr), SUM, 00332 sizeof(uint8_t)*depth); 00333 } 00334 } 00335 #if TIMESTUFF() 00336 /* time end */ 00337 #endif 00338 } 00339 else 00340 { 00341 phMemcpy(this->m_gauss, 00342 Image, 00343 width*height*depth*sizeof(uint8_t)); 00344 } 00345 00346 00347 if (this->m_useSobelApprox) 00348 { 00349 /* Sobel gradient appox. */ 00350 #if TIMESTUFF() 00351 /* phPROGRESS("Sobel gradient approximation\n"); */ 00352 /* time start */ 00353 #endif 00354 for( r = 2; r < (height - 2); r++ ) 00355 { 00356 row = r * row_length; 00357 for(c = 2; c < (width - 2); c++ ) 00358 { 00359 col = c * depth; 00360 00361 sumX = 0; 00362 sumY = 0; 00363 00364 /* image boundaries */ 00365 if((r == 0) || (r == (height-1))) 00366 { 00367 SUM = 0; 00368 } 00369 else if ((c == 0) || (c == (width-1))) 00370 { 00371 SUM = 0; 00372 } 00373 /* Convolution starts here */ 00374 else 00375 { 00376 /* pthoren - I think 'J' and 'I' are being applied backwards */ 00377 /* gradient approximations */ 00378 for( I = -1; I <= 1; I++ ) 00379 { 00380 inrow = (r + I) * row_length; 00381 for( J = -1; J <= 1; J++ ) 00382 { 00383 incol = (c + J) * depth; 00384 offset = inrow + incol; 00385 /* X gradient approximation */ 00386 sumX += (int)( (*(this->m_gauss + offset)) * GX[I+1][J+1]); 00387 /* Y gradient approximation */ 00388 sumY += (int)( (*(this->m_gauss + offset)) * GY[I+1][J+1]); 00389 } 00390 } 00391 00392 /* GRADIENT MAGNITUDE APPROXIMATION (Myler p.218) */ 00393 SUM = abs(sumX) + abs(sumY); 00394 00395 if ( SUM > 255 ) SUM=255; 00396 /* SUm can not be negative */ 00397 /* if ( SUM < 0 ) SUM=0; */ 00398 00399 /* Magnitude orientation */ 00400 /* Cannot divide by zero */ 00401 if(sumX == 0) 00402 { 00403 if ( sumY == 0 ) 00404 { 00405 ORIENT = 0.0; 00406 } 00407 else if (sumY<0) 00408 { 00409 sumY = -sumY; 00410 ORIENT = 90.0; 00411 } 00412 else 00413 { 00414 ORIENT = 90.0; 00415 } 00416 } 00417 /* Can't take invtan of angle in 2nd Quad */ 00418 else if ((sumX < 0) && (sumY > 0)) 00419 { 00420 sumX = -sumX; 00421 ORIENT = 180 - ((atan((float)(sumY)/(float)(sumX))) * (float)(180/M_PI)); 00422 } 00423 00424 /* Can't take invtan of angle in 4th Quad */ 00425 else if((sumX > 0) && (sumY < 0)) 00426 { 00427 sumY = -sumY; 00428 ORIENT = 180 - ((atan((float)(sumY)/(float)(sumX))) * (float)(180/M_PI)); 00429 } 00430 /* else angle is in 1st or 3rd Quad */ 00431 else 00432 { 00433 ORIENT = (atan((float)(sumY)/(float)(sumX))) * (float)(180/M_PI); 00434 } 00435 00436 /*----------------------------------------------------- 00437 * Find edgeDirection by assigning ORIENT a value of 00438 * either 0, 45, 90 or 135 degrees, depending on which 00439 * value ORIENT is closest to 00440 * ---------------------------------------------------*/ 00441 if (ORIENT < 22.5) 00442 edgeDirection = 0; 00443 00444 else if (ORIENT < 67.5) 00445 edgeDirection = 45; 00446 00447 else if (ORIENT < 112.5) 00448 edgeDirection = 90; 00449 00450 else if (ORIENT < 157.5) 00451 edgeDirection = 135; 00452 00453 else 00454 edgeDirection = 0; 00455 00456 /* Obtain values of 2 adjacent pixels in edge direction. */ 00457 addr = this->m_gauss + row + col; 00458 if(edgeDirection == 0) 00459 { 00460 /* one pixel to the left */ 00461 leftPixel = (int)(*((addr - depth))); 00462 /* one pixel to the rigby */ 00463 rightPixel = (int)(*((addr + depth))); 00464 } 00465 else if(edgeDirection == 45) 00466 { 00467 /* one row down, one pixel to the left */ 00468 leftPixel = (int)(*((addr + row_length - depth))); 00469 /* one row up, one pixel to the right */ 00470 rightPixel = (int)(*((addr - row_length + depth))); 00471 } 00472 else if(edgeDirection == 90) 00473 { 00474 /* one row up */ 00475 leftPixel = (int)(*((addr - row_length))); 00476 /* one row down */ 00477 rightPixel = (int)(*((addr + row_length))); 00478 } 00479 else 00480 { 00481 /* One row up, one pixel to the left */ 00482 leftPixel = (int)(*((addr - row_length - depth))); 00483 /* One row down, one pixel to the right */ 00484 rightPixel = (int)(*((addr + row_length + depth))); 00485 } 00486 /* ------------------------------------------------- 00487 * Compare current magnitude to both adjacent 00488 * pixel values. And if it is less than either 00489 * of the 2 adjacent values - suppress it and make 00490 * a nonedge. 00491 * ------------------------------------------------- */ 00492 00493 if ((SUM < leftPixel) || (SUM < rightPixel)) 00494 { 00495 SUM = 0; 00496 } 00497 else 00498 { 00499 /* Hysteresis */ 00500 if ( SUM >= this->m_highThreshold ) 00501 { 00502 SUM = 255; /* edge */ 00503 } 00504 else if ( SUM <= this->m_lowThreshold ) 00505 { 00506 SUM = 0; /* nonedge */ 00507 } 00508 /* SUM is between T1 & T2 */ 00509 else 00510 { 00511 /* This is the old way to do it, as you can see 00512 * there are a lot more calculations performed 00513 * in the following as opposed to the newer version 00514 * below. 00515 */ 00516 00517 #if 0 00518 addr = this->m_gauss + r * row_length; 00519 /* Determine values of neighboring pixels */ 00520 P1 = (int)(*(addr + ((c - width - 1) * depth))); 00521 P2 = (int)(*(addr + ((c - width) * depth))); 00522 P3 = (int)(*(addr + ((c - width + 1) * depth))); 00523 P4 = (int)(*(addr + ((c - 1) * depth))); 00524 P5 = (int)(*(addr + ((c + 1) * depth))); 00525 P6 = (int)(*(addr + ((c + width - 1) * depth))); 00526 P7 = (int)(*(addr + ((c + width) * depth))); 00527 P8 = (int)(*(addr + ((c + width + 1) * depth))); 00528 #else 00529 00530 /* Old way- subs:6 adds:15 muls:9 00531 * potential-loads: 33 00532 * 00533 * New way- subs:4 adds:7 muls:3 00534 * potential-loads: 28 00535 * 00536 * += requires 2 loads, one for the target before and 00537 * one for the operand. 00538 * 00539 * This could still be improved by making 'addr' 00540 * change very little by taking it's calculation from 00541 * here to outside the loop and incrementing it by 00542 * depth. TODO 00543 */ 00544 00545 addr = (this->m_gauss + r * row_length) - 00546 ((width - c - 1) * depth); 00547 00548 /* 00549 * + ((c - width - 1) * depth); 00550 * 00551 * Using the following when adding to the address 00552 * will not work as expected. The left-hand side will 00553 * evaluate to a negative number and doesn't get added 00554 * correctly to the address type. However, if you switch 00555 * to "width - c", the value will always be positive and 00556 * the addition is performed correctly. 00557 */ 00558 00559 /* [R, C ] * 00560 * [0,[0,1,2]] * 00561 * [1,[0,x,2]] * 00562 * [2,[0,1,2]] */ 00563 00564 /* [0,0] first pixel of kernel */ 00565 P1 = (int)*(addr); 00566 00567 /* [0,1] next column:pixel */ 00568 addr += depth; 00569 P2 = (int)*(addr); 00570 00571 /* [0,2] next column:pixel */ 00572 addr += depth; 00573 P3 = (int)*(addr); 00574 00575 /* [1,2] next row */ 00576 addr += row_length; 00577 P5 = (int)*(addr); 00578 00579 /* [1,0] first pixel of row */ 00580 addr -= 2 * depth; 00581 P4 = (int)*(addr); 00582 00583 /* [2,0] next row */ 00584 addr += row_length; 00585 P6 = (int)*(addr); 00586 00587 /* [2,1] next column:pixel */ 00588 addr += depth; 00589 P7 = (int)*(addr); 00590 00591 /* [2,2] next column:pixel */ 00592 addr += depth; 00593 P8 = (int)*(addr); 00594 #endif 00595 00596 /* Check to see if neighboring pixel values are edges */ 00597 if ((P1 > this->m_highThreshold) || 00598 (P2 > this->m_highThreshold) || 00599 (P3 > this->m_highThreshold) || 00600 (P4 > this->m_highThreshold) || 00601 (P5 > this->m_highThreshold) || 00602 (P6 > this->m_highThreshold) || 00603 (P7 > this->m_highThreshold) || 00604 (P8 > this->m_highThreshold)) 00605 { 00606 SUM = 0; /* make edge */ 00607 } 00608 else 00609 { 00610 SUM = 255; /* make nonedge */ 00611 } 00612 } 00613 } 00614 } /* else loop ends here (starts after b.c.) */ 00615 00616 /* Set an edge for all channels of the image */ 00617 phMemset((void *)(this->m_edge + row + col), 00618 255 - (uint8_t)(SUM), 00619 sizeof(uint8_t) * depth); 00620 } 00621 } 00622 00623 /* copy the edge image to the Image output */ 00624 phMemcpy(Image, 00625 this->m_edge, 00626 width * height * depth * sizeof(uint8_t)); 00627 00628 #if TIMESTUFF() 00629 /* time end */ 00630 #endif 00631 } 00632 else 00633 { 00634 /* If we just want to see the output of the 5x5 gaussian, just copy 00635 * the gaussian to the Image */ 00636 phMemcpy(Image, 00637 this->m_gauss, 00638 width * height * depth * sizeof(uint8_t)); 00639 } 00640 00641 /* End Filter */ 00642 00643 return phSUCCESS; 00644 00645 error: 00646 /* make sure to free up any space here, not including Image */ 00647 00648 return phFAIL; 00649 } 00650
| Copyright (C) 2002 - 2007 |
Philip D.S. Thoren ( pthoren@users.sourceforge.net ) University Of Massachusetts at Lowell Robotics Lab |
|
1.4.4