[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* The VIGRA Website is */ 00008 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00009 /* Please direct questions, bug reports, and contributions to */ 00010 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00011 /* vigra@informatik.uni-hamburg.de */ 00012 /* */ 00013 /* Permission is hereby granted, free of charge, to any person */ 00014 /* obtaining a copy of this software and associated documentation */ 00015 /* files (the "Software"), to deal in the Software without */ 00016 /* restriction, including without limitation the rights to use, */ 00017 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00018 /* sell copies of the Software, and to permit persons to whom the */ 00019 /* Software is furnished to do so, subject to the following */ 00020 /* conditions: */ 00021 /* */ 00022 /* The above copyright notice and this permission notice shall be */ 00023 /* included in all copies or substantial portions of the */ 00024 /* Software. */ 00025 /* */ 00026 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00027 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00028 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00029 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00030 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00031 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00032 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00033 /* OTHER DEALINGS IN THE SOFTWARE. */ 00034 /* */ 00035 /************************************************************************/ 00036 00037 00038 #ifndef VIGRA_EDGEDETECTION_HXX 00039 #define VIGRA_EDGEDETECTION_HXX 00040 00041 #include <vector> 00042 #include <queue> 00043 #include <cmath> // sqrt(), abs() 00044 #include "utilities.hxx" 00045 #include "numerictraits.hxx" 00046 #include "stdimage.hxx" 00047 #include "stdimagefunctions.hxx" 00048 #include "recursiveconvolution.hxx" 00049 #include "separableconvolution.hxx" 00050 #include "labelimage.hxx" 00051 #include "mathutil.hxx" 00052 #include "pixelneighborhood.hxx" 00053 #include "linear_solve.hxx" 00054 00055 00056 namespace vigra { 00057 00058 /** \addtogroup EdgeDetection Edge Detection 00059 Edge detectors based on first and second derivatives, 00060 and related post-processing. 00061 */ 00062 //@{ 00063 00064 /********************************************************/ 00065 /* */ 00066 /* differenceOfExponentialEdgeImage */ 00067 /* */ 00068 /********************************************************/ 00069 00070 /** \brief Detect and mark edges in an edge image using the Shen/Castan zero-crossing detector. 00071 00072 This operator applies an exponential filter to the source image 00073 at the given <TT>scale</TT> and subtracts the result from the original image. 00074 Zero crossings are detected in the resulting difference image. Whenever the 00075 gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>, 00076 an edge point is marked (using <TT>edge_marker</TT>) in the destination image on 00077 the darker side of the zero crossing (note that zero crossings occur 00078 <i>between</i> pixels). For example: 00079 00080 \code 00081 sign of difference image resulting edge points (*) 00082 00083 + - - * * . 00084 + + - => . * * 00085 + + + . . . 00086 \endcode 00087 00088 Non-edge pixels (<TT>.</TT>) remain untouched in the destination image. 00089 The result can be improved by the post-processing operation \ref removeShortEdges(). 00090 A more accurate edge placement can be achieved with the function 00091 \ref differenceOfExponentialCrackEdgeImage(). 00092 00093 The source value type 00094 (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition, 00095 subtraction and multiplication of the type with itself, and multiplication 00096 with double and 00097 \ref NumericTraits "NumericTraits" must 00098 be defined. In addition, this type must be less-comparable. 00099 00100 <b> Declarations:</b> 00101 00102 pass arguments explicitly: 00103 \code 00104 namespace vigra { 00105 template <class SrcIterator, class SrcAccessor, 00106 class DestIterator, class DestAccessor, 00107 class GradValue, 00108 class DestValue = DestAccessor::value_type> 00109 void differenceOfExponentialEdgeImage( 00110 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00111 DestIterator dul, DestAccessor da, 00112 double scale, GradValue gradient_threshold, 00113 DestValue edge_marker = NumericTraits<DestValue>::one()) 00114 } 00115 \endcode 00116 00117 use argument objects in conjunction with \ref ArgumentObjectFactories : 00118 \code 00119 namespace vigra { 00120 template <class SrcIterator, class SrcAccessor, 00121 class DestIterator, class DestAccessor, 00122 class GradValue, 00123 class DestValue = DestAccessor::value_type> 00124 void differenceOfExponentialEdgeImage( 00125 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00126 pair<DestIterator, DestAccessor> dest, 00127 double scale, GradValue gradient_threshold, 00128 DestValue edge_marker = NumericTraits<DestValue>::one()) 00129 } 00130 \endcode 00131 00132 <b> Usage:</b> 00133 00134 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 00135 Namespace: vigra 00136 00137 \code 00138 vigra::BImage src(w,h), edges(w,h); 00139 00140 // empty edge image 00141 edges = 0; 00142 ... 00143 00144 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 00145 vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 00146 0.8, 4.0, 1); 00147 \endcode 00148 00149 <b> Required Interface:</b> 00150 00151 \code 00152 SrcImageIterator src_upperleft, src_lowerright; 00153 DestImageIterator dest_upperleft; 00154 00155 SrcAccessor src_accessor; 00156 DestAccessor dest_accessor; 00157 00158 SrcAccessor::value_type u = src_accessor(src_upperleft); 00159 double d; 00160 GradValue gradient_threshold; 00161 00162 u = u + u 00163 u = u - u 00164 u = u * u 00165 u = d * u 00166 u < gradient_threshold 00167 00168 DestValue edge_marker; 00169 dest_accessor.set(edge_marker, dest_upperleft); 00170 \endcode 00171 00172 <b> Preconditions:</b> 00173 00174 \code 00175 scale > 0 00176 gradient_threshold > 0 00177 \endcode 00178 */ 00179 doxygen_overloaded_function(template <...> void differenceOfExponentialEdgeImage) 00180 00181 template <class SrcIterator, class SrcAccessor, 00182 class DestIterator, class DestAccessor, 00183 class GradValue, class DestValue> 00184 void differenceOfExponentialEdgeImage( 00185 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00186 DestIterator dul, DestAccessor da, 00187 double scale, GradValue gradient_threshold, DestValue edge_marker) 00188 { 00189 vigra_precondition(scale > 0, 00190 "differenceOfExponentialEdgeImage(): scale > 0 required."); 00191 00192 vigra_precondition(gradient_threshold > 0, 00193 "differenceOfExponentialEdgeImage(): " 00194 "gradient_threshold > 0 required."); 00195 00196 int w = slr.x - sul.x; 00197 int h = slr.y - sul.y; 00198 int x,y; 00199 00200 typedef typename 00201 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00202 TMPTYPE; 00203 typedef BasicImage<TMPTYPE> TMPIMG; 00204 00205 TMPIMG tmp(w,h); 00206 TMPIMG smooth(w,h); 00207 00208 recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0); 00209 recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0); 00210 00211 recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale); 00212 recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale); 00213 00214 typename TMPIMG::Iterator iy = smooth.upperLeft(); 00215 typename TMPIMG::Iterator ty = tmp.upperLeft(); 00216 DestIterator dy = dul; 00217 00218 static const Diff2D right(1, 0); 00219 static const Diff2D bottom(0, 1); 00220 00221 00222 TMPTYPE thresh = (gradient_threshold * gradient_threshold) * 00223 NumericTraits<TMPTYPE>::one(); 00224 TMPTYPE zero = NumericTraits<TMPTYPE>::zero(); 00225 00226 for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, ++dy.y) 00227 { 00228 typename TMPIMG::Iterator ix = iy; 00229 typename TMPIMG::Iterator tx = ty; 00230 DestIterator dx = dy; 00231 00232 for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x) 00233 { 00234 TMPTYPE diff = *tx - *ix; 00235 TMPTYPE gx = tx[right] - *tx; 00236 TMPTYPE gy = tx[bottom] - *tx; 00237 00238 if((gx * gx > thresh) && 00239 (diff * (tx[right] - ix[right]) < zero)) 00240 { 00241 if(gx < zero) 00242 { 00243 da.set(edge_marker, dx, right); 00244 } 00245 else 00246 { 00247 da.set(edge_marker, dx); 00248 } 00249 } 00250 if(((gy * gy > thresh) && 00251 (diff * (tx[bottom] - ix[bottom]) < zero))) 00252 { 00253 if(gy < zero) 00254 { 00255 da.set(edge_marker, dx, bottom); 00256 } 00257 else 00258 { 00259 da.set(edge_marker, dx); 00260 } 00261 } 00262 } 00263 } 00264 00265 typename TMPIMG::Iterator ix = iy; 00266 typename TMPIMG::Iterator tx = ty; 00267 DestIterator dx = dy; 00268 00269 for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x) 00270 { 00271 TMPTYPE diff = *tx - *ix; 00272 TMPTYPE gx = tx[right] - *tx; 00273 00274 if((gx * gx > thresh) && 00275 (diff * (tx[right] - ix[right]) < zero)) 00276 { 00277 if(gx < zero) 00278 { 00279 da.set(edge_marker, dx, right); 00280 } 00281 else 00282 { 00283 da.set(edge_marker, dx); 00284 } 00285 } 00286 } 00287 } 00288 00289 template <class SrcIterator, class SrcAccessor, 00290 class DestIterator, class DestAccessor, 00291 class GradValue> 00292 inline 00293 void differenceOfExponentialEdgeImage( 00294 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00295 DestIterator dul, DestAccessor da, 00296 double scale, GradValue gradient_threshold) 00297 { 00298 differenceOfExponentialEdgeImage(sul, slr, sa, dul, da, 00299 scale, gradient_threshold, 1); 00300 } 00301 00302 template <class SrcIterator, class SrcAccessor, 00303 class DestIterator, class DestAccessor, 00304 class GradValue, class DestValue> 00305 inline 00306 void differenceOfExponentialEdgeImage( 00307 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00308 pair<DestIterator, DestAccessor> dest, 00309 double scale, GradValue gradient_threshold, 00310 DestValue edge_marker) 00311 { 00312 differenceOfExponentialEdgeImage(src.first, src.second, src.third, 00313 dest.first, dest.second, 00314 scale, gradient_threshold, 00315 edge_marker); 00316 } 00317 00318 template <class SrcIterator, class SrcAccessor, 00319 class DestIterator, class DestAccessor, 00320 class GradValue> 00321 inline 00322 void differenceOfExponentialEdgeImage( 00323 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00324 pair<DestIterator, DestAccessor> dest, 00325 double scale, GradValue gradient_threshold) 00326 { 00327 differenceOfExponentialEdgeImage(src.first, src.second, src.third, 00328 dest.first, dest.second, 00329 scale, gradient_threshold, 1); 00330 } 00331 00332 /********************************************************/ 00333 /* */ 00334 /* differenceOfExponentialCrackEdgeImage */ 00335 /* */ 00336 /********************************************************/ 00337 00338 /** \brief Detect and mark edges in a crack edge image using the Shen/Castan zero-crossing detector. 00339 00340 This operator applies an exponential filter to the source image 00341 at the given <TT>scale</TT> and subtracts the result from the original image. 00342 Zero crossings are detected in the resulting difference image. Whenever the 00343 gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>, 00344 an edge point is marked (using <TT>edge_marker</TT>) in the destination image 00345 <i>between</i> the corresponding original pixels. Topologically, this means we 00346 must insert additional pixels between the original ones to represent the 00347 boundaries between the pixels (the so called zero- and one-cells, with the original 00348 pixels being two-cells). Within VIGRA, such an image is called \ref CrackEdgeImage. 00349 To allow insertion of the zero- and one-cells, the destination image must have twice the 00350 size of the original (precisely, <TT>(2*w-1)</TT> by <TT>(2*h-1)</TT> pixels). Then the algorithm 00351 proceeds as follows: 00352 00353 \code 00354 sign of difference image insert zero- and one-cells resulting edge points (*) 00355 00356 + . - . - . * . . . 00357 + - - . . . . . . * * * . 00358 + + - => + . + . - => . . . * . 00359 + + + . . . . . . . . * * 00360 + . + . + . . . . . 00361 \endcode 00362 00363 Thus the edge points are marked where they actually are - in between the pixels. 00364 An important property of the resulting edge image is that it conforms to the notion 00365 of well-composedness as defined by Latecki et al., i.e. connected regions and edges 00366 obtained by a subsequent \ref Labeling do not depend on 00367 whether 4- or 8-connectivity is used. 00368 The non-edge pixels (<TT>.</TT>) in the destination image remain unchanged. 00369 The result conformes to the requirements of a \ref CrackEdgeImage. It can be further 00370 improved by the post-processing operations \ref removeShortEdges() and 00371 \ref closeGapsInCrackEdgeImage(). 00372 00373 The source value type (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition, 00374 subtraction and multiplication of the type with itself, and multiplication 00375 with double and 00376 \ref NumericTraits "NumericTraits" must 00377 be defined. In addition, this type must be less-comparable. 00378 00379 <b> Declarations:</b> 00380 00381 pass arguments explicitly: 00382 \code 00383 namespace vigra { 00384 template <class SrcIterator, class SrcAccessor, 00385 class DestIterator, class DestAccessor, 00386 class GradValue, 00387 class DestValue = DestAccessor::value_type> 00388 void differenceOfExponentialCrackEdgeImage( 00389 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00390 DestIterator dul, DestAccessor da, 00391 double scale, GradValue gradient_threshold, 00392 DestValue edge_marker = NumericTraits<DestValue>::one()) 00393 } 00394 \endcode 00395 00396 use argument objects in conjunction with \ref ArgumentObjectFactories : 00397 \code 00398 namespace vigra { 00399 template <class SrcIterator, class SrcAccessor, 00400 class DestIterator, class DestAccessor, 00401 class GradValue, 00402 class DestValue = DestAccessor::value_type> 00403 void differenceOfExponentialCrackEdgeImage( 00404 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00405 pair<DestIterator, DestAccessor> dest, 00406 double scale, GradValue gradient_threshold, 00407 DestValue edge_marker = NumericTraits<DestValue>::one()) 00408 } 00409 \endcode 00410 00411 <b> Usage:</b> 00412 00413 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 00414 Namespace: vigra 00415 00416 \code 00417 vigra::BImage src(w,h), edges(2*w-1,2*h-1); 00418 00419 // empty edge image 00420 edges = 0; 00421 ... 00422 00423 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 00424 vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 00425 0.8, 4.0, 1); 00426 \endcode 00427 00428 <b> Required Interface:</b> 00429 00430 \code 00431 SrcImageIterator src_upperleft, src_lowerright; 00432 DestImageIterator dest_upperleft; 00433 00434 SrcAccessor src_accessor; 00435 DestAccessor dest_accessor; 00436 00437 SrcAccessor::value_type u = src_accessor(src_upperleft); 00438 double d; 00439 GradValue gradient_threshold; 00440 00441 u = u + u 00442 u = u - u 00443 u = u * u 00444 u = d * u 00445 u < gradient_threshold 00446 00447 DestValue edge_marker; 00448 dest_accessor.set(edge_marker, dest_upperleft); 00449 \endcode 00450 00451 <b> Preconditions:</b> 00452 00453 \code 00454 scale > 0 00455 gradient_threshold > 0 00456 \endcode 00457 00458 The destination image must have twice the size of the source: 00459 \code 00460 w_dest = 2 * w_src - 1 00461 h_dest = 2 * h_src - 1 00462 \endcode 00463 */ 00464 doxygen_overloaded_function(template <...> void differenceOfExponentialCrackEdgeImage) 00465 00466 template <class SrcIterator, class SrcAccessor, 00467 class DestIterator, class DestAccessor, 00468 class GradValue, class DestValue> 00469 void differenceOfExponentialCrackEdgeImage( 00470 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00471 DestIterator dul, DestAccessor da, 00472 double scale, GradValue gradient_threshold, 00473 DestValue edge_marker) 00474 { 00475 vigra_precondition(scale > 0, 00476 "differenceOfExponentialCrackEdgeImage(): scale > 0 required."); 00477 00478 vigra_precondition(gradient_threshold > 0, 00479 "differenceOfExponentialCrackEdgeImage(): " 00480 "gradient_threshold > 0 required."); 00481 00482 int w = slr.x - sul.x; 00483 int h = slr.y - sul.y; 00484 int x, y; 00485 00486 typedef typename 00487 NumericTraits<typename SrcAccessor::value_type>::RealPromote 00488 TMPTYPE; 00489 typedef BasicImage<TMPTYPE> TMPIMG; 00490 00491 TMPIMG tmp(w,h); 00492 TMPIMG smooth(w,h); 00493 00494 TMPTYPE zero = NumericTraits<TMPTYPE>::zero(); 00495 00496 static const Diff2D right(1,0); 00497 static const Diff2D bottom(0,1); 00498 static const Diff2D left(-1,0); 00499 static const Diff2D top(0,-1); 00500 00501 recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0); 00502 recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0); 00503 00504 recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale); 00505 recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale); 00506 00507 typename TMPIMG::Iterator iy = smooth.upperLeft(); 00508 typename TMPIMG::Iterator ty = tmp.upperLeft(); 00509 DestIterator dy = dul; 00510 00511 TMPTYPE thresh = (gradient_threshold * gradient_threshold) * 00512 NumericTraits<TMPTYPE>::one(); 00513 00514 // find zero crossings above threshold 00515 for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, dy.y+=2) 00516 { 00517 typename TMPIMG::Iterator ix = iy; 00518 typename TMPIMG::Iterator tx = ty; 00519 DestIterator dx = dy; 00520 00521 for(int x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2) 00522 { 00523 TMPTYPE diff = *tx - *ix; 00524 TMPTYPE gx = tx[right] - *tx; 00525 TMPTYPE gy = tx[bottom] - *tx; 00526 00527 if((gx * gx > thresh) && 00528 (diff * (tx[right] - ix[right]) < zero)) 00529 { 00530 da.set(edge_marker, dx, right); 00531 } 00532 if((gy * gy > thresh) && 00533 (diff * (tx[bottom] - ix[bottom]) < zero)) 00534 { 00535 da.set(edge_marker, dx, bottom); 00536 } 00537 } 00538 00539 TMPTYPE diff = *tx - *ix; 00540 TMPTYPE gy = tx[bottom] - *tx; 00541 00542 if((gy * gy > thresh) && 00543 (diff * (tx[bottom] - ix[bottom]) < zero)) 00544 { 00545 da.set(edge_marker, dx, bottom); 00546 } 00547 } 00548 00549 typename TMPIMG::Iterator ix = iy; 00550 typename TMPIMG::Iterator tx = ty; 00551 DestIterator dx = dy; 00552 00553 for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2) 00554 { 00555 TMPTYPE diff = *tx - *ix; 00556 TMPTYPE gx = tx[right] - *tx; 00557 00558 if((gx * gx > thresh) && 00559 (diff * (tx[right] - ix[right]) < zero)) 00560 { 00561 da.set(edge_marker, dx, right); 00562 } 00563 } 00564 00565 iy = smooth.upperLeft() + Diff2D(0,1); 00566 ty = tmp.upperLeft() + Diff2D(0,1); 00567 dy = dul + Diff2D(1,2); 00568 00569 static const Diff2D topleft(-1,-1); 00570 static const Diff2D topright(1,-1); 00571 static const Diff2D bottomleft(-1,1); 00572 static const Diff2D bottomright(1,1); 00573 00574 // find missing 1-cells below threshold (x-direction) 00575 for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2) 00576 { 00577 typename TMPIMG::Iterator ix = iy; 00578 typename TMPIMG::Iterator tx = ty; 00579 DestIterator dx = dy; 00580 00581 for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2) 00582 { 00583 if(da(dx) == edge_marker) continue; 00584 00585 TMPTYPE diff = *tx - *ix; 00586 00587 if((diff * (tx[right] - ix[right]) < zero) && 00588 (((da(dx, bottomright) == edge_marker) && 00589 (da(dx, topleft) == edge_marker)) || 00590 ((da(dx, bottomleft) == edge_marker) && 00591 (da(dx, topright) == edge_marker)))) 00592 00593 { 00594 da.set(edge_marker, dx); 00595 } 00596 } 00597 } 00598 00599 iy = smooth.upperLeft() + Diff2D(1,0); 00600 ty = tmp.upperLeft() + Diff2D(1,0); 00601 dy = dul + Diff2D(2,1); 00602 00603 // find missing 1-cells below threshold (y-direction) 00604 for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2) 00605 { 00606 typename TMPIMG::Iterator ix = iy; 00607 typename TMPIMG::Iterator tx = ty; 00608 DestIterator dx = dy; 00609 00610 for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2) 00611 { 00612 if(da(dx) == edge_marker) continue; 00613 00614 TMPTYPE diff = *tx - *ix; 00615 00616 if((diff * (tx[bottom] - ix[bottom]) < zero) && 00617 (((da(dx, bottomright) == edge_marker) && 00618 (da(dx, topleft) == edge_marker)) || 00619 ((da(dx, bottomleft) == edge_marker) && 00620 (da(dx, topright) == edge_marker)))) 00621 00622 { 00623 da.set(edge_marker, dx); 00624 } 00625 } 00626 } 00627 00628 dy = dul + Diff2D(1,1); 00629 00630 // find missing 0-cells 00631 for(y=0; y<h-1; ++y, dy.y+=2) 00632 { 00633 DestIterator dx = dy; 00634 00635 for(int x=0; x<w-1; ++x, dx.x+=2) 00636 { 00637 static const Diff2D dist[] = {right, top, left, bottom }; 00638 00639 int i; 00640 for(i=0; i<4; ++i) 00641 { 00642 if(da(dx, dist[i]) == edge_marker) break; 00643 } 00644 00645 if(i < 4) da.set(edge_marker, dx); 00646 } 00647 } 00648 } 00649 00650 template <class SrcIterator, class SrcAccessor, 00651 class DestIterator, class DestAccessor, 00652 class GradValue, class DestValue> 00653 inline 00654 void differenceOfExponentialCrackEdgeImage( 00655 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00656 pair<DestIterator, DestAccessor> dest, 00657 double scale, GradValue gradient_threshold, 00658 DestValue edge_marker) 00659 { 00660 differenceOfExponentialCrackEdgeImage(src.first, src.second, src.third, 00661 dest.first, dest.second, 00662 scale, gradient_threshold, 00663 edge_marker); 00664 } 00665 00666 /********************************************************/ 00667 /* */ 00668 /* removeShortEdges */ 00669 /* */ 00670 /********************************************************/ 00671 00672 /** \brief Remove short edges from an edge image. 00673 00674 This algorithm can be applied as a post-processing operation of 00675 \ref differenceOfExponentialEdgeImage() and \ref differenceOfExponentialCrackEdgeImage(). 00676 It removes all edges that are shorter than <TT>min_edge_length</TT>. The corresponding 00677 pixels are set to the <TT>non_edge_marker</TT>. The idea behind this algorithms is 00678 that very short edges are probably caused by noise and don't represent interesting 00679 image structure. Technically, the algorithms executes a connected components labeling, 00680 so the image's value type must be equality comparable. 00681 00682 If the source image fulfills the requirements of a \ref CrackEdgeImage, 00683 it will still do so after application of this algorithm. 00684 00685 Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 00686 i.e. on only one image. Also, the algorithm assumes that all non-edges pixels are already 00687 marked with the given <TT>non_edge_marker</TT> value. 00688 00689 <b> Declarations:</b> 00690 00691 pass arguments explicitly: 00692 \code 00693 namespace vigra { 00694 template <class Iterator, class Accessor, class SrcValue> 00695 void removeShortEdges( 00696 Iterator sul, Iterator slr, Accessor sa, 00697 int min_edge_length, SrcValue non_edge_marker) 00698 } 00699 \endcode 00700 00701 use argument objects in conjunction with \ref ArgumentObjectFactories : 00702 \code 00703 namespace vigra { 00704 template <class Iterator, class Accessor, class SrcValue> 00705 void removeShortEdges( 00706 triple<Iterator, Iterator, Accessor> src, 00707 int min_edge_length, SrcValue non_edge_marker) 00708 } 00709 \endcode 00710 00711 <b> Usage:</b> 00712 00713 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 00714 Namespace: vigra 00715 00716 \code 00717 vigra::BImage src(w,h), edges(w,h); 00718 00719 // empty edge image 00720 edges = 0; 00721 ... 00722 00723 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 00724 vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 00725 0.8, 4.0, 1); 00726 00727 // zero edges shorter than 10 pixels 00728 vigra::removeShortEdges(srcImageRange(edges), 10, 0); 00729 \endcode 00730 00731 <b> Required Interface:</b> 00732 00733 \code 00734 SrcImageIterator src_upperleft, src_lowerright; 00735 DestImageIterator dest_upperleft; 00736 00737 SrcAccessor src_accessor; 00738 DestAccessor dest_accessor; 00739 00740 SrcAccessor::value_type u = src_accessor(src_upperleft); 00741 00742 u == u 00743 00744 SrcValue non_edge_marker; 00745 src_accessor.set(non_edge_marker, src_upperleft); 00746 \endcode 00747 */ 00748 doxygen_overloaded_function(template <...> void removeShortEdges) 00749 00750 template <class Iterator, class Accessor, class Value> 00751 void removeShortEdges( 00752 Iterator sul, Iterator slr, Accessor sa, 00753 unsigned int min_edge_length, Value non_edge_marker) 00754 { 00755 int w = slr.x - sul.x; 00756 int h = slr.y - sul.y; 00757 int x,y; 00758 00759 IImage labels(w, h); 00760 labels = 0; 00761 00762 int number_of_regions = 00763 labelImageWithBackground(srcIterRange(sul,slr,sa), 00764 destImage(labels), true, non_edge_marker); 00765 00766 ArrayOfRegionStatistics<FindROISize<int> > 00767 region_stats(number_of_regions); 00768 00769 inspectTwoImages(srcImageRange(labels), srcImage(labels), region_stats); 00770 00771 IImage::Iterator ly = labels.upperLeft(); 00772 Iterator oy = sul; 00773 00774 for(y=0; y<h; ++y, ++oy.y, ++ly.y) 00775 { 00776 Iterator ox(oy); 00777 IImage::Iterator lx(ly); 00778 00779 for(x=0; x<w; ++x, ++ox.x, ++lx.x) 00780 { 00781 if(sa(ox) == non_edge_marker) continue; 00782 if((region_stats[*lx].count) < min_edge_length) 00783 { 00784 sa.set(non_edge_marker, ox); 00785 } 00786 } 00787 } 00788 } 00789 00790 template <class Iterator, class Accessor, class Value> 00791 inline 00792 void removeShortEdges( 00793 triple<Iterator, Iterator, Accessor> src, 00794 unsigned int min_edge_length, Value non_edge_marker) 00795 { 00796 removeShortEdges(src.first, src.second, src.third, 00797 min_edge_length, non_edge_marker); 00798 } 00799 00800 /********************************************************/ 00801 /* */ 00802 /* closeGapsInCrackEdgeImage */ 00803 /* */ 00804 /********************************************************/ 00805 00806 /** \brief Close one-pixel wide gaps in a cell grid edge image. 00807 00808 This algorithm is typically applied as a post-processing operation of 00809 \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill 00810 the requirements of a \ref CrackEdgeImage, and will still do so after 00811 application of this algorithm. 00812 00813 It closes one pixel wide gaps in the edges resulting from this algorithm. 00814 Since these gaps are usually caused by zero crossing slightly below the gradient 00815 threshold used in edge detection, this algorithms acts like a weak hysteresis 00816 thresholding. The newly found edge pixels are marked with the given <TT>edge_marker</TT>. 00817 The image's value type must be equality comparable. 00818 00819 Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 00820 i.e. on only one image. 00821 00822 <b> Declarations:</b> 00823 00824 pass arguments explicitly: 00825 \code 00826 namespace vigra { 00827 template <class SrcIterator, class SrcAccessor, class SrcValue> 00828 void closeGapsInCrackEdgeImage( 00829 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00830 SrcValue edge_marker) 00831 } 00832 \endcode 00833 00834 use argument objects in conjunction with \ref ArgumentObjectFactories : 00835 \code 00836 namespace vigra { 00837 template <class SrcIterator, class SrcAccessor, class SrcValue> 00838 void closeGapsInCrackEdgeImage( 00839 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00840 SrcValue edge_marker) 00841 } 00842 \endcode 00843 00844 <b> Usage:</b> 00845 00846 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 00847 Namespace: vigra 00848 00849 \code 00850 vigra::BImage src(w,h), edges(2*w-1, 2*h-1); 00851 00852 // empty edge image 00853 edges = 0; 00854 ... 00855 00856 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 00857 vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 00858 0.8, 4.0, 1); 00859 00860 // close gaps, mark with 1 00861 vigra::closeGapsInCrackEdgeImage(srcImageRange(edges), 1); 00862 00863 // zero edges shorter than 20 pixels 00864 vigra::removeShortEdges(srcImageRange(edges), 10, 0); 00865 \endcode 00866 00867 <b> Required Interface:</b> 00868 00869 \code 00870 SrcImageIterator src_upperleft, src_lowerright; 00871 00872 SrcAccessor src_accessor; 00873 DestAccessor dest_accessor; 00874 00875 SrcAccessor::value_type u = src_accessor(src_upperleft); 00876 00877 u == u 00878 u != u 00879 00880 SrcValue edge_marker; 00881 src_accessor.set(edge_marker, src_upperleft); 00882 \endcode 00883 */ 00884 doxygen_overloaded_function(template <...> void closeGapsInCrackEdgeImage) 00885 00886 template <class SrcIterator, class SrcAccessor, class SrcValue> 00887 void closeGapsInCrackEdgeImage( 00888 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 00889 SrcValue edge_marker) 00890 { 00891 int w = (slr.x - sul.x) / 2; 00892 int h = (slr.y - sul.y) / 2; 00893 int x, y; 00894 00895 int count1, count2, count3; 00896 00897 static const Diff2D right(1,0); 00898 static const Diff2D bottom(0,1); 00899 static const Diff2D left(-1,0); 00900 static const Diff2D top(0,-1); 00901 00902 static const Diff2D leftdist[] = { 00903 Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)}; 00904 static const Diff2D rightdist[] = { 00905 Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)}; 00906 static const Diff2D topdist[] = { 00907 Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)}; 00908 static const Diff2D bottomdist[] = { 00909 Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)}; 00910 00911 int i; 00912 00913 SrcIterator sy = sul + Diff2D(0,1); 00914 SrcIterator sx; 00915 00916 // close 1-pixel wide gaps (x-direction) 00917 for(y=0; y<h; ++y, sy.y+=2) 00918 { 00919 sx = sy + Diff2D(2,0); 00920 00921 for(x=2; x<w; ++x, sx.x+=2) 00922 { 00923 if(sa(sx) == edge_marker) continue; 00924 00925 if(sa(sx, left) != edge_marker) continue; 00926 if(sa(sx, right) != edge_marker) continue; 00927 00928 count1 = 0; 00929 count2 = 0; 00930 count3 = 0; 00931 00932 for(i=0; i<4; ++i) 00933 { 00934 if(sa(sx, leftdist[i]) == edge_marker) 00935 { 00936 ++count1; 00937 count3 ^= 1 << i; 00938 } 00939 if(sa(sx, rightdist[i]) == edge_marker) 00940 { 00941 ++count2; 00942 count3 ^= 1 << i; 00943 } 00944 } 00945 00946 if(count1 <= 1 || count2 <= 1 || count3 == 15) 00947 { 00948 sa.set(edge_marker, sx); 00949 } 00950 } 00951 } 00952 00953 sy = sul + Diff2D(1,2); 00954 00955 // close 1-pixel wide gaps (y-direction) 00956 for(y=2; y<h; ++y, sy.y+=2) 00957 { 00958 sx = sy; 00959 00960 for(x=0; x<w; ++x, sx.x+=2) 00961 { 00962 if(sa(sx) == edge_marker) continue; 00963 00964 if(sa(sx, top) != edge_marker) continue; 00965 if(sa(sx, bottom) != edge_marker) continue; 00966 00967 count1 = 0; 00968 count2 = 0; 00969 count3 = 0; 00970 00971 for(i=0; i<4; ++i) 00972 { 00973 if(sa(sx, topdist[i]) == edge_marker) 00974 { 00975 ++count1; 00976 count3 ^= 1 << i; 00977 } 00978 if(sa(sx, bottomdist[i]) == edge_marker) 00979 { 00980 ++count2; 00981 count3 ^= 1 << i; 00982 } 00983 } 00984 00985 if(count1 <= 1 || count2 <= 1 || count3 == 15) 00986 { 00987 sa.set(edge_marker, sx); 00988 } 00989 } 00990 } 00991 } 00992 00993 template <class SrcIterator, class SrcAccessor, class SrcValue> 00994 inline 00995 void closeGapsInCrackEdgeImage( 00996 triple<SrcIterator, SrcIterator, SrcAccessor> src, 00997 SrcValue edge_marker) 00998 { 00999 closeGapsInCrackEdgeImage(src.first, src.second, src.third, 01000 edge_marker); 01001 } 01002 01003 /********************************************************/ 01004 /* */ 01005 /* beautifyCrackEdgeImage */ 01006 /* */ 01007 /********************************************************/ 01008 01009 /** \brief Beautify crack edge image for visualization. 01010 01011 This algorithm is applied as a post-processing operation of 01012 \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill 01013 the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after 01014 application of this algorithm. In particular, the algorithm removes zero-cells 01015 marked as edges to avoid staircase effects on diagonal lines like this: 01016 01017 \code 01018 original edge points (*) resulting edge points 01019 01020 . * . . . . * . . . 01021 . * * * . . . * . . 01022 . . . * . => . . . * . 01023 . . . * * . . . . * 01024 . . . . . . . . . . 01025 \endcode 01026 01027 Therfore, this algorithm should only be applied as a vizualization aid, i.e. 01028 for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>, 01029 and background pixels with <TT>background_marker</TT>. The image's value type must be 01030 equality comparable. 01031 01032 Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 01033 i.e. on only one image. 01034 01035 <b> Declarations:</b> 01036 01037 pass arguments explicitly: 01038 \code 01039 namespace vigra { 01040 template <class SrcIterator, class SrcAccessor, class SrcValue> 01041 void beautifyCrackEdgeImage( 01042 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01043 SrcValue edge_marker, SrcValue background_marker) 01044 } 01045 \endcode 01046 01047 use argument objects in conjunction with \ref ArgumentObjectFactories : 01048 \code 01049 namespace vigra { 01050 template <class SrcIterator, class SrcAccessor, class SrcValue> 01051 void beautifyCrackEdgeImage( 01052 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01053 SrcValue edge_marker, SrcValue background_marker) 01054 } 01055 \endcode 01056 01057 <b> Usage:</b> 01058 01059 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 01060 Namespace: vigra 01061 01062 \code 01063 vigra::BImage src(w,h), edges(2*w-1, 2*h-1); 01064 01065 // empty edge image 01066 edges = 0; 01067 ... 01068 01069 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 01070 vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 01071 0.8, 4.0, 1); 01072 01073 // beautify edge image for visualization 01074 vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0); 01075 01076 // show to the user 01077 window.open(edges); 01078 \endcode 01079 01080 <b> Required Interface:</b> 01081 01082 \code 01083 SrcImageIterator src_upperleft, src_lowerright; 01084 01085 SrcAccessor src_accessor; 01086 DestAccessor dest_accessor; 01087 01088 SrcAccessor::value_type u = src_accessor(src_upperleft); 01089 01090 u == u 01091 u != u 01092 01093 SrcValue background_marker; 01094 src_accessor.set(background_marker, src_upperleft); 01095 \endcode 01096 */ 01097 doxygen_overloaded_function(template <...> void beautifyCrackEdgeImage) 01098 01099 template <class SrcIterator, class SrcAccessor, class SrcValue> 01100 void beautifyCrackEdgeImage( 01101 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01102 SrcValue edge_marker, SrcValue background_marker) 01103 { 01104 int w = (slr.x - sul.x) / 2; 01105 int h = (slr.y - sul.y) / 2; 01106 int x, y; 01107 01108 SrcIterator sy = sul + Diff2D(1,1); 01109 SrcIterator sx; 01110 01111 static const Diff2D right(1,0); 01112 static const Diff2D bottom(0,1); 01113 static const Diff2D left(-1,0); 01114 static const Diff2D top(0,-1); 01115 01116 // delete 0-cells at corners 01117 for(y=0; y<h; ++y, sy.y+=2) 01118 { 01119 sx = sy; 01120 01121 for(x=0; x<w; ++x, sx.x+=2) 01122 { 01123 if(sa(sx) != edge_marker) continue; 01124 01125 if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue; 01126 if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue; 01127 01128 sa.set(background_marker, sx); 01129 } 01130 } 01131 } 01132 01133 template <class SrcIterator, class SrcAccessor, class SrcValue> 01134 inline 01135 void beautifyCrackEdgeImage( 01136 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01137 SrcValue edge_marker, SrcValue background_marker) 01138 { 01139 beautifyCrackEdgeImage(src.first, src.second, src.third, 01140 edge_marker, background_marker); 01141 } 01142 01143 01144 /** Helper class that stores edgel attributes. 01145 */ 01146 class Edgel 01147 { 01148 public: 01149 /** The edgel's sub-pixel x coordinate. 01150 */ 01151 float x; 01152 01153 /** The edgel's sub-pixel y coordinate. 01154 */ 01155 float y; 01156 01157 /** The edgel's strength (magnitude of the gradient vector). 01158 */ 01159 float strength; 01160 01161 /** 01162 The edgel's orientation. This is the angle 01163 between the x-axis and the edge, so that the bright side of the 01164 edge is on the right. The angle is measured 01165 counter-clockwise in radians like this: 01166 01167 01168 \code 01169 01170 edgel axis 01171 \ (bright side) 01172 (dark \ 01173 side) \ /__ 01174 \\ \ orientation angle 01175 \ | 01176 +------------> x-axis 01177 | 01178 | 01179 | 01180 | 01181 y-axis V 01182 \endcode 01183 01184 So, for example a vertical edge with its dark side on the left 01185 has orientation PI/2, and a horizontal edge with dark side on top 01186 has orientation 0. Obviously, the edge's orientation changes 01187 by PI if the contrast is reversed. 01188 01189 */ 01190 float orientation; 01191 01192 Edgel() 01193 : x(0.0f), y(0.0f), strength(0.0f), orientation(0.0f) 01194 {} 01195 01196 Edgel(float ix, float iy, float is, float io) 01197 : x(ix), y(iy), strength(is), orientation(io) 01198 {} 01199 }; 01200 01201 template <class Image1, class Image2, class BackInsertable> 01202 void internalCannyFindEdgels(Image1 const & gx, 01203 Image1 const & gy, 01204 Image2 const & magnitude, 01205 BackInsertable & edgels) 01206 { 01207 typedef typename Image1::value_type PixelType; 01208 double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0); 01209 01210 for(int y=1; y<gx.height()-1; ++y) 01211 { 01212 for(int x=1; x<gx.width()-1; ++x) 01213 { 01214 PixelType gradx = gx(x,y); 01215 PixelType grady = gy(x,y); 01216 double mag = magnitude(x, y); 01217 if(mag == 0.0) 01218 continue; 01219 01220 int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5); 01221 int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5); 01222 01223 int x1 = x - dx, 01224 x2 = x + dx, 01225 y1 = y - dy, 01226 y2 = y + dy; 01227 01228 PixelType m1 = magnitude(x1, y1); 01229 PixelType m3 = magnitude(x2, y2); 01230 01231 if(m1 < mag && m3 <= mag) 01232 { 01233 Edgel edgel; 01234 01235 // local maximum => quadratic interpolation of sub-pixel location 01236 PixelType del = (m1 - m3) / 2.0 / (m1 + m3 - 2.0*mag); 01237 edgel.x = x + dx*del; 01238 edgel.y = y + dy*del; 01239 edgel.strength = mag; 01240 double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5; 01241 if(orientation < 0.0) 01242 orientation += 2.0*M_PI; 01243 edgel.orientation = orientation; 01244 edgels.push_back(edgel); 01245 } 01246 } 01247 } 01248 } 01249 01250 /********************************************************/ 01251 /* */ 01252 /* cannyEdgelList */ 01253 /* */ 01254 /********************************************************/ 01255 01256 /** \brief Simple implementation of Canny's edge detector. 01257 01258 This operator first calculates the gradient vector for each 01259 pixel of the image using first derivatives of a Gaussian at the 01260 given scale. Then a very simple non-maxima supression is performed: 01261 for each 3x3 neighborhood, it is determined whether the center pixel has 01262 larger gradient magnitude than its two neighbors in gradient direction 01263 (where the direction is rounded into octands). If this is the case, 01264 a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel 01265 edgel position is determined by fitting a parabola 01266 to the three gradient magnitude values 01267 mentioned above. The sub-pixel location of the parabola's tip 01268 and the gradient magnitude and direction (from the pixel center) 01269 are written in the newly created edgel. 01270 01271 <b> Declarations:</b> 01272 01273 pass arguments explicitly: 01274 \code 01275 namespace vigra { 01276 template <class SrcIterator, class SrcAccessor, class BackInsertable> 01277 void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src, 01278 BackInsertable & edgels, double scale); 01279 } 01280 \endcode 01281 01282 use argument objects in conjunction with \ref ArgumentObjectFactories : 01283 \code 01284 namespace vigra { 01285 template <class SrcIterator, class SrcAccessor, class BackInsertable> 01286 void 01287 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01288 BackInsertable & edgels, double scale); 01289 } 01290 \endcode 01291 01292 <b> Usage:</b> 01293 01294 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 01295 Namespace: vigra 01296 01297 \code 01298 vigra::BImage src(w,h); 01299 01300 // empty edgel list 01301 std::vector<vigra::Edgel> edgels; 01302 ... 01303 01304 // find edgels at scale 0.8 01305 vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8); 01306 \endcode 01307 01308 <b> Required Interface:</b> 01309 01310 \code 01311 SrcImageIterator src_upperleft; 01312 SrcAccessor src_accessor; 01313 01314 src_accessor(src_upperleft); 01315 01316 BackInsertable edgels; 01317 edgels.push_back(Edgel()); 01318 \endcode 01319 01320 SrcAccessor::value_type must be a type convertible to float 01321 01322 <b> Preconditions:</b> 01323 01324 \code 01325 scale > 0 01326 \endcode 01327 */ 01328 doxygen_overloaded_function(template <...> void cannyEdgelList) 01329 01330 template <class SrcIterator, class SrcAccessor, class BackInsertable> 01331 void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src, 01332 BackInsertable & edgels, double scale) 01333 { 01334 int w = lr.x - ul.x; 01335 int h = lr.y - ul.y; 01336 01337 // calculate image gradients 01338 typedef typename 01339 NumericTraits<typename SrcAccessor::value_type>::RealPromote 01340 TmpType; 01341 01342 BasicImage<TmpType> tmp(w,h), dx(w,h), dy(w,h); 01343 01344 Kernel1D<double> smooth, grad; 01345 01346 smooth.initGaussian(scale); 01347 grad.initGaussianDerivative(scale, 1); 01348 01349 separableConvolveX(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad)); 01350 separableConvolveY(srcImageRange(tmp), destImage(dx), kernel1d(smooth)); 01351 01352 separableConvolveY(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad)); 01353 separableConvolveX(srcImageRange(tmp), destImage(dy), kernel1d(smooth)); 01354 01355 combineTwoImages(srcImageRange(dx), srcImage(dy), destImage(tmp), 01356 MagnitudeFunctor<TmpType>()); 01357 01358 01359 // find edgels 01360 internalCannyFindEdgels(dx, dy, tmp, edgels); 01361 } 01362 01363 template <class SrcIterator, class SrcAccessor, class BackInsertable> 01364 inline void 01365 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src, 01366 BackInsertable & edgels, double scale) 01367 { 01368 cannyEdgelList(src.first, src.second, src.third, edgels, scale); 01369 } 01370 01371 /********************************************************/ 01372 /* */ 01373 /* cannyEdgeImage */ 01374 /* */ 01375 /********************************************************/ 01376 01377 /** \brief Detect and mark edges in an edge image using Canny's algorithm. 01378 01379 This operator first calls \ref cannyEdgelList() to generate an 01380 edgel list for the given image. Then it scans this list and selects edgels 01381 whose strength is above the given <TT>gradient_threshold</TT>. For each of these 01382 edgels, the edgel's location is rounded to the nearest pixel, and that 01383 pixel marked with the given <TT>edge_marker</TT>. 01384 01385 <b> Declarations:</b> 01386 01387 pass arguments explicitly: 01388 \code 01389 namespace vigra { 01390 template <class SrcIterator, class SrcAccessor, 01391 class DestIterator, class DestAccessor, 01392 class GradValue, class DestValue> 01393 void cannyEdgeImage( 01394 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01395 DestIterator dul, DestAccessor da, 01396 double scale, GradValue gradient_threshold, DestValue edge_marker); 01397 } 01398 \endcode 01399 01400 use argument objects in conjunction with \ref ArgumentObjectFactories : 01401 \code 01402 namespace vigra { 01403 template <class SrcIterator, class SrcAccessor, 01404 class DestIterator, class DestAccessor, 01405 class GradValue, class DestValue> 01406 void cannyEdgeImage( 01407 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01408 pair<DestIterator, DestAccessor> dest, 01409 double scale, GradValue gradient_threshold, DestValue edge_marker); 01410 } 01411 \endcode 01412 01413 <b> Usage:</b> 01414 01415 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 01416 Namespace: vigra 01417 01418 \code 01419 vigra::BImage src(w,h), edges(w,h); 01420 01421 // empty edge image 01422 edges = 0; 01423 ... 01424 01425 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 01426 vigra::cannyEdgeImage(srcImageRange(src), destImage(edges), 01427 0.8, 4.0, 1); 01428 \endcode 01429 01430 <b> Required Interface:</b> 01431 01432 see also: \ref cannyEdgelList(). 01433 01434 \code 01435 DestImageIterator dest_upperleft; 01436 DestAccessor dest_accessor; 01437 DestValue edge_marker; 01438 01439 dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1)); 01440 \endcode 01441 01442 <b> Preconditions:</b> 01443 01444 \code 01445 scale > 0 01446 gradient_threshold > 0 01447 \endcode 01448 */ 01449 doxygen_overloaded_function(template <...> void cannyEdgeImage) 01450 01451 template <class SrcIterator, class SrcAccessor, 01452 class DestIterator, class DestAccessor, 01453 class GradValue, class DestValue> 01454 void cannyEdgeImage( 01455 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01456 DestIterator dul, DestAccessor da, 01457 double scale, GradValue gradient_threshold, DestValue edge_marker) 01458 { 01459 std::vector<Edgel> edgels; 01460 01461 cannyEdgelList(sul, slr, sa, edgels, scale); 01462 01463 for(unsigned int i=0; i<edgels.size(); ++i) 01464 { 01465 if(gradient_threshold < edgels[i].strength) 01466 { 01467 Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5)); 01468 01469 da.set(edge_marker, dul, pix); 01470 } 01471 } 01472 } 01473 01474 template <class SrcIterator, class SrcAccessor, 01475 class DestIterator, class DestAccessor, 01476 class GradValue, class DestValue> 01477 inline void cannyEdgeImage( 01478 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01479 pair<DestIterator, DestAccessor> dest, 01480 double scale, GradValue gradient_threshold, DestValue edge_marker) 01481 { 01482 cannyEdgeImage(src.first, src.second, src.third, 01483 dest.first, dest.second, 01484 scale, gradient_threshold, edge_marker); 01485 } 01486 01487 /********************************************************/ 01488 01489 namespace detail { 01490 01491 template <class DestIterator> 01492 int neighborhoodConfiguration(DestIterator dul) 01493 { 01494 int v = 0; 01495 NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast); 01496 for(int i=0; i<8; ++i, --c) 01497 { 01498 v = (v << 1) | ((*c != 0) ? 1 : 0); 01499 } 01500 01501 return v; 01502 } 01503 01504 template <class GradValue> 01505 struct SimplePoint 01506 { 01507 Diff2D point; 01508 GradValue grad; 01509 01510 SimplePoint(Diff2D const & p, GradValue g) 01511 : point(p), grad(g) 01512 {} 01513 01514 bool operator<(SimplePoint const & o) const 01515 { 01516 return grad < o.grad; 01517 } 01518 01519 bool operator>(SimplePoint const & o) const 01520 { 01521 return grad > o.grad; 01522 } 01523 }; 01524 01525 template <class SrcIterator, class SrcAccessor, 01526 class DestIterator, class DestAccessor, 01527 class GradValue, class DestValue> 01528 void cannyEdgeImageFromGrad( 01529 SrcIterator sul, SrcIterator slr, SrcAccessor grad, 01530 DestIterator dul, DestAccessor da, 01531 GradValue gradient_threshold, DestValue edge_marker) 01532 { 01533 typedef typename SrcAccessor::value_type PixelType; 01534 typedef typename NormTraits<PixelType>::SquaredNormType NormType; 01535 01536 NormType zero = NumericTraits<NormType>::zero(); 01537 double tan22_5 = M_SQRT2 - 1.0; 01538 typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold); 01539 01540 int w = slr.x - sul.x; 01541 int h = slr.y - sul.y; 01542 01543 sul += Diff2D(1,1); 01544 dul += Diff2D(1,1); 01545 Diff2D p(0,0); 01546 01547 for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y) 01548 { 01549 SrcIterator sx = sul; 01550 DestIterator dx = dul; 01551 for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x) 01552 { 01553 PixelType g = grad(sx); 01554 NormType g2n = squaredNorm(g); 01555 if(g2n < g2thresh) 01556 continue; 01557 01558 NormType g2n1, g2n3; 01559 // find out quadrant 01560 if(abs(g[1]) < tan22_5*abs(g[0])) 01561 { 01562 // north-south edge 01563 g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0))); 01564 g2n3 = squaredNorm(grad(sx, Diff2D(1, 0))); 01565 } 01566 else if(abs(g[0]) < tan22_5*abs(g[1])) 01567 { 01568 // west-east edge 01569 g2n1 = squaredNorm(grad(sx, Diff2D(0, -1))); 01570 g2n3 = squaredNorm(grad(sx, Diff2D(0, 1))); 01571 } 01572 else if(g[0]*g[1] < zero) 01573 { 01574 // north-west-south-east edge 01575 g2n1 = squaredNorm(grad(sx, Diff2D(1, -1))); 01576 g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1))); 01577 } 01578 else 01579 { 01580 // north-east-south-west edge 01581 g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1))); 01582 g2n3 = squaredNorm(grad(sx, Diff2D(1, 1))); 01583 } 01584 01585 if(g2n1 < g2n && g2n3 <= g2n) 01586 { 01587 da.set(edge_marker, dx); 01588 } 01589 } 01590 } 01591 } 01592 01593 } // namespace detail 01594 01595 /********************************************************/ 01596 /* */ 01597 /* cannyEdgeImageWithThinning */ 01598 /* */ 01599 /********************************************************/ 01600 01601 /** \brief Detect and mark edges in an edge image using Canny's algorithm. 01602 01603 The input pixels of this algorithms must be vectors of length 2 (see Required Interface below). 01604 It first searches for all pixels whose gradient magnitude is larger 01605 than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors 01606 in gradient direction (where these neighbors are determined by nearest neighbor 01607 interpolation, i.e. according to the octant where the gradient points into). 01608 The resulting edge pixel candidates are then subjected to topological thinning 01609 so that the remaining edge pixels can be linked into edgel chains with a provable, 01610 non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient 01611 magnitude survive. Optionally, the outermost pixels are marked as edge pixels 01612 as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination 01613 image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched). 01614 01615 <b> Declarations:</b> 01616 01617 pass arguments explicitly: 01618 \code 01619 namespace vigra { 01620 template <class SrcIterator, class SrcAccessor, 01621 class DestIterator, class DestAccessor, 01622 class GradValue, class DestValue> 01623 void cannyEdgeImageFromGradWithThinning( 01624 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01625 DestIterator dul, DestAccessor da, 01626 GradValue gradient_threshold, 01627 DestValue edge_marker, bool addBorder = true); 01628 } 01629 \endcode 01630 01631 use argument objects in conjunction with \ref ArgumentObjectFactories : 01632 \code 01633 namespace vigra { 01634 template <class SrcIterator, class SrcAccessor, 01635 class DestIterator, class DestAccessor, 01636 class GradValue, class DestValue> 01637 void cannyEdgeImageFromGradWithThinning( 01638 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01639 pair<DestIterator, DestAccessor> dest, 01640 GradValue gradient_threshold, 01641 DestValue edge_marker, bool addBorder = true); 01642 } 01643 \endcode 01644 01645 <b> Usage:</b> 01646 01647 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 01648 Namespace: vigra 01649 01650 \code 01651 vigra::BImage src(w,h), edges(w,h); 01652 01653 vigra::FVector2Image grad(w,h); 01654 // compute the image gradient at scale 0.8 01655 vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8); 01656 01657 // empty edge image 01658 edges = 0; 01659 // find edges gradient larger than 4.0, mark with 1, and add border 01660 vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges), 01661 4.0, 1, true); 01662 \endcode 01663 01664 <b> Required Interface:</b> 01665 01666 \code 01667 // the input pixel type must be a vector with two elements 01668 SrcImageIterator src_upperleft; 01669 SrcAccessor src_accessor; 01670 typedef SrcAccessor::value_type SrcPixel; 01671 typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType; 01672 01673 SrcPixel g = src_accessor(src_upperleft); 01674 SrcPixel::value_type g0 = g[0]; 01675 SrcSquaredNormType gn = squaredNorm(g); 01676 01677 DestImageIterator dest_upperleft; 01678 DestAccessor dest_accessor; 01679 DestValue edge_marker; 01680 01681 dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1)); 01682 \endcode 01683 01684 <b> Preconditions:</b> 01685 01686 \code 01687 gradient_threshold > 0 01688 \endcode 01689 */ 01690 doxygen_overloaded_function(template <...> void cannyEdgeImageFromGradWithThinning) 01691 01692 template <class SrcIterator, class SrcAccessor, 01693 class DestIterator, class DestAccessor, 01694 class GradValue, class DestValue> 01695 void cannyEdgeImageFromGradWithThinning( 01696 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01697 DestIterator dul, DestAccessor da, 01698 GradValue gradient_threshold, 01699 DestValue edge_marker, bool addBorder) 01700 { 01701 int w = slr.x - sul.x; 01702 int h = slr.y - sul.y; 01703 01704 BImage edgeImage(w, h, BImage::value_type(0)); 01705 BImage::traverser eul = edgeImage.upperLeft(); 01706 BImage::Accessor ea = edgeImage.accessor(); 01707 if(addBorder) 01708 initImageBorder(destImageRange(edgeImage), 1, 1); 01709 detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1); 01710 01711 static bool isSimplePoint[256] = { 01712 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 01713 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 01714 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 01715 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 01716 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 01717 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 01718 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 01719 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 01720 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 01721 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 01722 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 01723 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 01724 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 01725 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 01726 1, 0, 1, 0 }; 01727 01728 eul += Diff2D(1,1); 01729 sul += Diff2D(1,1); 01730 int w2 = w-2; 01731 int h2 = h-2; 01732 01733 typedef detail::SimplePoint<GradValue> SP; 01734 // use std::greater becaus we need the smallest gradients at the top of the queue 01735 std::priority_queue<SP, std::vector<SP>, std::greater<SP> > pqueue; 01736 01737 Diff2D p(0,0); 01738 for(; p.y < h2; ++p.y) 01739 { 01740 for(p.x = 0; p.x < w2; ++p.x) 01741 { 01742 BImage::traverser e = eul + p; 01743 if(*e == 0) 01744 continue; 01745 int v = detail::neighborhoodConfiguration(e); 01746 if(isSimplePoint[v]) 01747 { 01748 pqueue.push(SP(p, norm(sa(sul+p)))); 01749 *e = 2; // remember that it is already in queue 01750 } 01751 } 01752 } 01753 01754 static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1), 01755 Diff2D(1,0), Diff2D(0,1) }; 01756 01757 while(pqueue.size()) 01758 { 01759 p = pqueue.top().point; 01760 pqueue.pop(); 01761 01762 BImage::traverser e = eul + p; 01763 int v = detail::neighborhoodConfiguration(e); 01764 if(!isSimplePoint[v]) 01765 continue; // point may no longer be simple because its neighbors changed 01766 01767 *e = 0; // delete simple point 01768 01769 for(int i=0; i<4; ++i) 01770 { 01771 Diff2D pneu = p + dist[i]; 01772 if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2) 01773 continue; // do not remove points at the border 01774 01775 BImage::traverser eneu = eul + pneu; 01776 if(*eneu == 1) // point is boundary and not yet in the queue 01777 { 01778 int v = detail::neighborhoodConfiguration(eneu); 01779 if(isSimplePoint[v]) 01780 { 01781 pqueue.push(SP(pneu, norm(sa(sul+pneu)))); 01782 *eneu = 2; // remember that it is already in queue 01783 } 01784 } 01785 } 01786 } 01787 01788 initImageIf(destIterRange(dul, dul+Diff2D(w,h), da), 01789 maskImage(edgeImage), edge_marker); 01790 } 01791 01792 template <class SrcIterator, class SrcAccessor, 01793 class DestIterator, class DestAccessor, 01794 class GradValue, class DestValue> 01795 inline void cannyEdgeImageFromGradWithThinning( 01796 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01797 pair<DestIterator, DestAccessor> dest, 01798 GradValue gradient_threshold, 01799 DestValue edge_marker, bool addBorder) 01800 { 01801 cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third, 01802 dest.first, dest.second, 01803 gradient_threshold, edge_marker, addBorder); 01804 } 01805 01806 template <class SrcIterator, class SrcAccessor, 01807 class DestIterator, class DestAccessor, 01808 class GradValue, class DestValue> 01809 inline void cannyEdgeImageFromGradWithThinning( 01810 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01811 DestIterator dul, DestAccessor da, 01812 GradValue gradient_threshold, DestValue edge_marker) 01813 { 01814 cannyEdgeImageFromGradWithThinning(sul, slr, sa, 01815 dul, da, 01816 gradient_threshold, edge_marker, true); 01817 } 01818 01819 template <class SrcIterator, class SrcAccessor, 01820 class DestIterator, class DestAccessor, 01821 class GradValue, class DestValue> 01822 inline void cannyEdgeImageFromGradWithThinning( 01823 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01824 pair<DestIterator, DestAccessor> dest, 01825 GradValue gradient_threshold, DestValue edge_marker) 01826 { 01827 cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third, 01828 dest.first, dest.second, 01829 gradient_threshold, edge_marker, true); 01830 } 01831 01832 /********************************************************/ 01833 /* */ 01834 /* cannyEdgeImageWithThinning */ 01835 /* */ 01836 /********************************************************/ 01837 01838 /** \brief Detect and mark edges in an edge image using Canny's algorithm. 01839 01840 This operator first calls \ref gaussianGradient() to compute the gradient of the input 01841 image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an 01842 edge image. See there for more detailed documentation. 01843 01844 <b> Declarations:</b> 01845 01846 pass arguments explicitly: 01847 \code 01848 namespace vigra { 01849 template <class SrcIterator, class SrcAccessor, 01850 class DestIterator, class DestAccessor, 01851 class GradValue, class DestValue> 01852 void cannyEdgeImageWithThinning( 01853 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01854 DestIterator dul, DestAccessor da, 01855 double scale, GradValue gradient_threshold, 01856 DestValue edge_marker, bool addBorder = true); 01857 } 01858 \endcode 01859 01860 use argument objects in conjunction with \ref ArgumentObjectFactories : 01861 \code 01862 namespace vigra { 01863 template <class SrcIterator, class SrcAccessor, 01864 class DestIterator, class DestAccessor, 01865 class GradValue, class DestValue> 01866 void cannyEdgeImageWithThinning( 01867 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01868 pair<DestIterator, DestAccessor> dest, 01869 double scale, GradValue gradient_threshold, 01870 DestValue edge_marker, bool addBorder = true); 01871 } 01872 \endcode 01873 01874 <b> Usage:</b> 01875 01876 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 01877 Namespace: vigra 01878 01879 \code 01880 vigra::BImage src(w,h), edges(w,h); 01881 01882 // empty edge image 01883 edges = 0; 01884 ... 01885 01886 // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border 01887 vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges), 01888 0.8, 4.0, 1, true); 01889 \endcode 01890 01891 <b> Required Interface:</b> 01892 01893 see also: \ref cannyEdgelList(). 01894 01895 \code 01896 DestImageIterator dest_upperleft; 01897 DestAccessor dest_accessor; 01898 DestValue edge_marker; 01899 01900 dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1)); 01901 \endcode 01902 01903 <b> Preconditions:</b> 01904 01905 \code 01906 scale > 0 01907 gradient_threshold > 0 01908 \endcode 01909 */ 01910 doxygen_overloaded_function(template <...> void cannyEdgeImageWithThinning) 01911 01912 template <class SrcIterator, class SrcAccessor, 01913 class DestIterator, class DestAccessor, 01914 class GradValue, class DestValue> 01915 void cannyEdgeImageWithThinning( 01916 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01917 DestIterator dul, DestAccessor da, 01918 double scale, GradValue gradient_threshold, 01919 DestValue edge_marker, bool addBorder) 01920 { 01921 // mark pixels that are higher than their neighbors in gradient direction 01922 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 01923 BasicImage<TinyVector<TmpType, 2> > grad(slr-sul); 01924 gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale); 01925 cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da), 01926 gradient_threshold, edge_marker, addBorder); 01927 } 01928 01929 template <class SrcIterator, class SrcAccessor, 01930 class DestIterator, class DestAccessor, 01931 class GradValue, class DestValue> 01932 inline void cannyEdgeImageWithThinning( 01933 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01934 pair<DestIterator, DestAccessor> dest, 01935 double scale, GradValue gradient_threshold, 01936 DestValue edge_marker, bool addBorder) 01937 { 01938 cannyEdgeImageWithThinning(src.first, src.second, src.third, 01939 dest.first, dest.second, 01940 scale, gradient_threshold, edge_marker, addBorder); 01941 } 01942 01943 template <class SrcIterator, class SrcAccessor, 01944 class DestIterator, class DestAccessor, 01945 class GradValue, class DestValue> 01946 inline void cannyEdgeImageWithThinning( 01947 SrcIterator sul, SrcIterator slr, SrcAccessor sa, 01948 DestIterator dul, DestAccessor da, 01949 double scale, GradValue gradient_threshold, DestValue edge_marker) 01950 { 01951 cannyEdgeImageWithThinning(sul, slr, sa, 01952 dul, da, 01953 scale, gradient_threshold, edge_marker, true); 01954 } 01955 01956 template <class SrcIterator, class SrcAccessor, 01957 class DestIterator, class DestAccessor, 01958 class GradValue, class DestValue> 01959 inline void cannyEdgeImageWithThinning( 01960 triple<SrcIterator, SrcIterator, SrcAccessor> src, 01961 pair<DestIterator, DestAccessor> dest, 01962 double scale, GradValue gradient_threshold, DestValue edge_marker) 01963 { 01964 cannyEdgeImageWithThinning(src.first, src.second, src.third, 01965 dest.first, dest.second, 01966 scale, gradient_threshold, edge_marker, true); 01967 } 01968 01969 /********************************************************/ 01970 01971 template <class Image1, class Image2, class BackInsertable> 01972 void internalCannyFindEdgels3x3(Image1 const & grad, 01973 Image2 const & mask, 01974 BackInsertable & edgels) 01975 { 01976 typedef typename Image1::value_type PixelType; 01977 typedef typename PixelType::value_type ValueType; 01978 01979 for(int y=1; y<grad.height()-1; ++y) 01980 { 01981 for(int x=1; x<grad.width()-1; ++x) 01982 { 01983 if(!mask(x,y)) 01984 continue; 01985 01986 ValueType gradx = grad(x,y)[0]; 01987 ValueType grady = grad(x,y)[1]; 01988 double mag = hypot(gradx, grady); 01989 if(mag == 0.0) 01990 continue; 01991 double c = gradx / mag, 01992 s = grady / mag; 01993 01994 Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1); 01995 l(0,0) = 1.0; 01996 01997 for(int yy = -1; yy <= 1; ++yy) 01998 { 01999 for(int xx = -1; xx <= 1; ++xx) 02000 { 02001 double u = c*xx + s*yy; 02002 double v = norm(grad(x+xx, y+yy)); 02003 l(1,0) = u; 02004 l(2,0) = u*u; 02005 ml += outer(l); 02006 mr += v*l; 02007 } 02008 } 02009 02010 linearSolve(ml, mr, r); 02011 02012 Edgel edgel; 02013 02014 // local maximum => quadratic interpolation of sub-pixel location 02015 double del = -r(1,0) / 2.0 / r(2,0); 02016 if(std::fabs(del) > 1.5) // don't move by more than about a pixel diameter 02017 del = 0.0; 02018 edgel.x = x + c*del; 02019 edgel.y = y + s*del; 02020 edgel.strength = mag; 02021 double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5; 02022 if(orientation < 0.0) 02023 orientation += 2.0*M_PI; 02024 edgel.orientation = orientation; 02025 edgels.push_back(edgel); 02026 } 02027 } 02028 } 02029 02030 02031 /********************************************************/ 02032 /* */ 02033 /* cannyEdgelList3x3 */ 02034 /* */ 02035 /********************************************************/ 02036 02037 /** \brief Improved implementation of Canny's edge detector. 02038 02039 This operator first computes pixels which are crossed by the edge using 02040 cannyEdgeImageWithThinning(). The gradient magnitude in the 3x3 neighborhood of these 02041 pixels are then projected onto the normal of the edge (as determined 02042 by the gradient direction). The edgel's subpixel location is found by fitting a 02043 parabola through the 9 gradient values and determining the parabola's tip. 02044 A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola 02045 is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher. 02046 02047 <b> Declarations:</b> 02048 02049 pass arguments explicitly: 02050 \code 02051 namespace vigra { 02052 template <class SrcIterator, class SrcAccessor, class BackInsertable> 02053 void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src, 02054 BackInsertable & edgels, double scale); 02055 } 02056 \endcode 02057 02058 use argument objects in conjunction with \ref ArgumentObjectFactories : 02059 \code 02060 namespace vigra { 02061 template <class SrcIterator, class SrcAccessor, class BackInsertable> 02062 void 02063 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src, 02064 BackInsertable & edgels, double scale); 02065 } 02066 \endcode 02067 02068 <b> Usage:</b> 02069 02070 <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br> 02071 Namespace: vigra 02072 02073 \code 02074 vigra::BImage src(w,h); 02075 02076 // empty edgel list 02077 std::vector<vigra::Edgel> edgels; 02078 ... 02079 02080 // find edgels at scale 0.8 02081 vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8); 02082 \endcode 02083 02084 <b> Required Interface:</b> 02085 02086 \code 02087 SrcImageIterator src_upperleft; 02088 SrcAccessor src_accessor; 02089 02090 src_accessor(src_upperleft); 02091 02092 BackInsertable edgels; 02093 edgels.push_back(Edgel()); 02094 \endcode 02095 02096 SrcAccessor::value_type must be a type convertible to float 02097 02098 <b> Preconditions:</b> 02099 02100 \code 02101 scale > 0 02102 \endcode 02103 */ 02104 doxygen_overloaded_function(template <...> void cannyEdgelList3x3) 02105 02106 template <class SrcIterator, class SrcAccessor, class BackInsertable> 02107 void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src, 02108 BackInsertable & edgels, double scale) 02109 { 02110 int w = lr.x - ul.x; 02111 int h = lr.y - ul.y; 02112 02113 typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType; 02114 BasicImage<TinyVector<TmpType, 2> > grad(lr-ul); 02115 gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale); 02116 02117 UInt8Image edges(lr-ul); 02118 cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges), 02119 0.0, 1, false); 02120 02121 // find edgels 02122 internalCannyFindEdgels3x3(grad, edges, edgels); 02123 } 02124 02125 template <class SrcIterator, class SrcAccessor, class BackInsertable> 02126 inline void 02127 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src, 02128 BackInsertable & edgels, double scale) 02129 { 02130 cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale); 02131 } 02132 02133 02134 02135 //@} 02136 02137 /** \page CrackEdgeImage Crack Edge Image 02138 02139 Crack edges are marked <i>between</i> the pixels of an image. 02140 A Crack Edge Image is an image that represents these edges. In order 02141 to accomodate the cracks, the Crack Edge Image must be twice as large 02142 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image 02143 can easily be derived from a binary image or from the signs of the 02144 response of a Laplacean filter. Consider the following sketch, where 02145 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and 02146 <TT>*</TT> the resulting crack edges. 02147 02148 \code 02149 sign of difference image insert cracks resulting CrackEdgeImage 02150 02151 + . - . - . * . . . 02152 + - - . . . . . . * * * . 02153 + + - => + . + . - => . . . * . 02154 + + + . . . . . . . . * * 02155 + . + . + . . . . . 02156 \endcode 02157 02158 Starting from the original binary image (left), we insert crack pixels 02159 to get to the double-sized image (center). Finally, we mark all 02160 crack pixels whose non-crack neighbors have different signs as 02161 crack edge points, while all other pixels (crack and non-crack) become 02162 region pixels. 02163 02164 <b>Requirements on a Crack Edge Image:</b> 02165 02166 <ul> 02167 <li>Crack Edge Images have odd width and height. 02168 <li>Crack pixels have at least one odd coordinate. 02169 <li>Only crack pixels may be marked as edge points. 02170 <li>Crack pixels with two odd coordinates must be marked as edge points 02171 whenever any of their neighboring crack pixels was marked. 02172 </ul> 02173 02174 The last two requirements ensure that both edges and regions are 4-connected. 02175 Thus, 4-connectivity and 8-connectivity yield identical connected 02176 components in a Crack Edge Image (so called <i>well-composedness</i>). 02177 This ensures that Crack Edge Images have nice topological properties 02178 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000). 02179 */ 02180 02181 02182 } // namespace vigra 02183 02184 #endif // VIGRA_EDGEDETECTION_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|