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