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