[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]

details vigra/edgedetection.hxx VIGRA

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.3.2, Jan 27 2005 )                                    */
00008 /*    You may use, modify, and distribute this software according       */
00009 /*    to the terms stated in the LICENSE file included in               */
00010 /*    the VIGRA distribution.                                           */
00011 /*                                                                      */
00012 /*    The VIGRA Website is                                              */
00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00014 /*    Please direct questions, bug reports, and contributions to        */
00015 /*        koethe@informatik.uni-hamburg.de                              */
00016 /*                                                                      */
00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
00020 /*                                                                      */
00021 /************************************************************************/
00022  
00023  
00024 #ifndef VIGRA_EDGEDETECTION_HXX
00025 #define VIGRA_EDGEDETECTION_HXX
00026 
00027 #include <vector>
00028 #include <cmath>     // sqrt(), abs()
00029 #include "vigra/utilities.hxx"
00030 #include "vigra/numerictraits.hxx"
00031 #include "vigra/stdimage.hxx"
00032 #include "vigra/stdimagefunctions.hxx"
00033 #include "vigra/recursiveconvolution.hxx"
00034 #include "vigra/separableconvolution.hxx"
00035 #include "vigra/labelimage.hxx"
00036 
00037 
00038 namespace vigra {
00039 
00040 /** \addtogroup EdgeDetection Edge Detection
00041     Edge detectors based on first and second derivatives,
00042           and related post-processing.
00043 */
00044 //@{ 
00045                                     
00046 /********************************************************/
00047 /*                                                      */
00048 /*           differenceOfExponentialEdgeImage           */
00049 /*                                                      */
00050 /********************************************************/
00051 
00052 /** \brief Detect and mark edges in an edge image using the Shen/Castan zero-crossing detector.
00053 
00054     This operator applies an exponential filter to the source image 
00055     at the given <TT>scale</TT> and subtracts the result from the original image. 
00056     Zero crossings are detected in the resulting difference image. Whenever the
00057     gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
00058     an edge point is marked (using <TT>edge_marker</TT>) in the destination image on
00059     the darker side of the zero crossing (note that zero crossings occur 
00060     <i>between</i> pixels). For example:
00061     
00062     \code
00063     sign of difference image     resulting edge points (*)
00064     
00065         + - -                          * * .
00066         + + -               =>         . * *
00067         + + +                          . . .
00068     \endcode
00069     
00070     Non-edge pixels (<TT>.</TT>) remain untouched in the destination image. 
00071     The result can be improved by the post-processing operation \ref removeShortEdges().
00072     A more accurate edge placement can be achieved with the function 
00073     \ref differenceOfExponentialCrackEdgeImage(). 
00074 
00075     The source value type 
00076     (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition, 
00077     subtraction and multiplication of the type with itself, and multiplication 
00078     with double and 
00079     \ref NumericTraits "NumericTraits" must 
00080     be defined. In addition, this type must be less-comparable.
00081     
00082     <b> Declarations:</b>
00083     
00084     pass arguments explicitly:
00085     \code
00086     namespace vigra {
00087         template <class SrcIterator, class SrcAccessor,
00088               class DestIterator, class DestAccessor,
00089               class GradValue,
00090               class DestValue = DestAccessor::value_type>
00091         void differenceOfExponentialEdgeImage(
00092                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00093                DestIterator dul, DestAccessor da,
00094                double scale, GradValue gradient_threshold, 
00095                DestValue edge_marker = NumericTraits<DestValue>::one())
00096     }
00097     \endcode
00098     
00099     use argument objects in conjunction with \ref ArgumentObjectFactories:
00100     \code
00101     namespace vigra {
00102         template <class SrcIterator, class SrcAccessor,
00103               class DestIterator, class DestAccessor, 
00104               class GradValue,
00105               class DestValue = DestAccessor::value_type>
00106         inline 
00107         void differenceOfExponentialEdgeImage(
00108                triple<SrcIterator, SrcIterator, SrcAccessor> src,
00109                pair<DestIterator, DestAccessor> dest,
00110                double scale, GradValue gradient_threshold,
00111                DestValue edge_marker = NumericTraits<DestValue>::one())
00112     }
00113     \endcode
00114     
00115     <b> Usage:</b>
00116     
00117         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00118     Namespace: vigra
00119     
00120     \code
00121     vigra::BImage src(w,h), edges(w,h);
00122     
00123     // empty edge image
00124     edges = 0;
00125     ...
00126     
00127     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00128     vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 
00129                                      0.8, 4.0, 1);
00130     \endcode
00131 
00132     <b> Required Interface:</b>
00133     
00134     \code
00135     SrcImageIterator src_upperleft, src_lowerright;
00136     DestImageIterator dest_upperleft;
00137     
00138     SrcAccessor src_accessor;
00139     DestAccessor dest_accessor;
00140     
00141     SrcAccessor::value_type u = src_accessor(src_upperleft);
00142     double d;
00143     GradValue gradient_threshold;
00144     
00145     u = u + u
00146     u = u - u
00147     u = u * u
00148     u = d * u
00149     u < gradient_threshold
00150     
00151     DestValue edge_marker;
00152     dest_accessor.set(edge_marker, dest_upperleft);
00153     \endcode
00154     
00155     <b> Preconditions:</b>
00156     
00157     \code
00158     scale > 0
00159     gradient_threshold > 0
00160     \endcode
00161 */
00162 template <class SrcIterator, class SrcAccessor,
00163           class DestIterator, class DestAccessor, 
00164           class GradValue, class DestValue>
00165 void differenceOfExponentialEdgeImage(
00166                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00167            DestIterator dul, DestAccessor da,
00168            double scale, GradValue gradient_threshold, DestValue edge_marker)
00169 {
00170     vigra_precondition(scale > 0,
00171                  "differenceOfExponentialEdgeImage(): scale > 0 required.");
00172          
00173     vigra_precondition(gradient_threshold > 0,
00174                  "differenceOfExponentialEdgeImage(): "
00175          "gradient_threshold > 0 required.");
00176          
00177     int w = slr.x - sul.x;
00178     int h = slr.y - sul.y;
00179     int x,y;
00180 
00181     typedef typename 
00182         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00183     TMPTYPE;
00184     typedef BasicImage<TMPTYPE> TMPIMG;
00185 
00186     TMPIMG tmp(w,h);
00187     TMPIMG smooth(w,h);
00188     
00189     recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
00190     recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
00191 
00192     recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
00193     recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
00194 
00195     typename TMPIMG::Iterator iy = smooth.upperLeft();
00196     typename TMPIMG::Iterator ty = tmp.upperLeft();
00197     DestIterator              dy = dul;
00198     
00199     static const Diff2D right(1, 0);
00200     static const Diff2D bottom(0, 1);
00201     
00202     
00203     TMPTYPE thresh = (gradient_threshold * gradient_threshold) * 
00204                      NumericTraits<TMPTYPE>::one();
00205     TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
00206 
00207     for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, ++dy.y)
00208     {
00209         typename TMPIMG::Iterator ix = iy;
00210         typename TMPIMG::Iterator tx = ty;
00211         DestIterator              dx = dy;
00212 
00213         for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
00214         {
00215             TMPTYPE diff = *tx - *ix;
00216             TMPTYPE gx = tx[right] - *tx;
00217             TMPTYPE gy = tx[bottom] - *tx;
00218 
00219             if((gx * gx > thresh) &&
00220                 (diff * (tx[right] - ix[right]) < zero))
00221             {
00222                 if(gx < zero)
00223                 {
00224                     da.set(edge_marker, dx, right);
00225                 }
00226                 else
00227                 {
00228                     da.set(edge_marker, dx);
00229                 }
00230             }
00231             if(((gy * gy > thresh) &&
00232                 (diff * (tx[bottom] - ix[bottom]) < zero)))
00233             {
00234                 if(gy < zero)
00235                 {
00236                     da.set(edge_marker, dx, bottom);
00237                 }
00238                 else
00239                 {
00240                     da.set(edge_marker, dx);
00241                 }
00242             }
00243         }
00244     }
00245     
00246     typename TMPIMG::Iterator ix = iy;
00247     typename TMPIMG::Iterator tx = ty;
00248     DestIterator              dx = dy;
00249     
00250     for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
00251     {
00252         TMPTYPE diff = *tx - *ix;
00253         TMPTYPE gx = tx[right] - *tx;
00254 
00255         if((gx * gx > thresh) &&
00256            (diff * (tx[right] - ix[right]) < zero))
00257         {
00258             if(gx < zero)
00259             {
00260                 da.set(edge_marker, dx, right);
00261             }
00262             else
00263             {
00264                 da.set(edge_marker, dx);
00265             }
00266         }
00267     }
00268 }
00269 
00270 template <class SrcIterator, class SrcAccessor,
00271           class DestIterator, class DestAccessor,
00272           class GradValue>
00273 inline 
00274 void differenceOfExponentialEdgeImage(
00275                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00276            DestIterator dul, DestAccessor da,
00277            double scale, GradValue gradient_threshold)
00278 {
00279     differenceOfExponentialEdgeImage(sul, slr, sa, dul, da, 
00280                                         scale, gradient_threshold, 1);
00281 }
00282 
00283 template <class SrcIterator, class SrcAccessor,
00284           class DestIterator, class DestAccessor, 
00285       class GradValue, class DestValue>
00286 inline 
00287 void differenceOfExponentialEdgeImage(
00288            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00289        pair<DestIterator, DestAccessor> dest,
00290        double scale, GradValue gradient_threshold,
00291        DestValue edge_marker)
00292 {
00293     differenceOfExponentialEdgeImage(src.first, src.second, src.third,
00294                                         dest.first, dest.second,
00295                     scale, gradient_threshold,
00296                     edge_marker);
00297 }
00298 
00299 template <class SrcIterator, class SrcAccessor,
00300           class DestIterator, class DestAccessor,
00301           class GradValue>
00302 inline 
00303 void differenceOfExponentialEdgeImage(
00304            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00305        pair<DestIterator, DestAccessor> dest,
00306        double scale, GradValue gradient_threshold)
00307 {
00308     differenceOfExponentialEdgeImage(src.first, src.second, src.third,
00309                                         dest.first, dest.second,
00310                     scale, gradient_threshold, 1);
00311 }
00312 
00313 /********************************************************/
00314 /*                                                      */
00315 /*        differenceOfExponentialCrackEdgeImage         */
00316 /*                                                      */
00317 /********************************************************/
00318 
00319 /** \brief Detect and mark edges in a crack edge image using the Shen/Castan zero-crossing detector.
00320 
00321     This operator applies an exponential filter to the source image 
00322     at the given <TT>scale</TT> and subtracts the result from the original image. 
00323     Zero crossings are detected in the resulting difference image. Whenever the
00324     gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
00325     an edge point is marked (using <TT>edge_marker</TT>) in the destination image 
00326     <i>between</i> the corresponding original pixels. Topologically, this means we 
00327     must insert additional pixels between the original ones to represent the
00328     boundaries between the pixels (the so called zero- and one-cells, with the original
00329     pixels being two-cells). Within VIGRA, such an image is called \ref CrackEdgeImage.
00330     To allow insertion of the zero- and one-cells, the destination image must have twice the 
00331     size of the original (precisely, <TT>(2*w-1)</TT> by <TT>(2*h-1)</TT> pixels). Then the algorithm 
00332     proceeds as follows:
00333     
00334     \code
00335 sign of difference image     insert zero- and one-cells     resulting edge points (*)
00336 
00337                                      + . - . -                   . * . . .
00338       + - -                          . . . . .                   . * * * .
00339       + + -               =>         + . + . -           =>      . . . * .
00340       + + +                          . . . . .                   . . . * *
00341                                      + . + . +                   . . . . .
00342     \endcode
00343     
00344     Thus the edge points are marked where they actually are - in between the pixels. 
00345     An important property of the resulting edge image is that it conforms to the notion 
00346     of well-composedness as defined by Latecki et al., i.e. connected regions and edges 
00347     obtained by a subsequent \ref Labeling do not depend on 
00348     whether 4- or 8-connectivity is used.
00349     The non-edge pixels (<TT>.</TT>) in the destination image remain unchanged. 
00350     The result conformes to the requirements of a \ref CrackEdgeImage. It can be further
00351     improved by the post-processing operations \ref removeShortEdges() and
00352     \ref closeGapsInCrackEdgeImage().
00353     
00354     The source value type (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition, 
00355     subtraction and multiplication of the type with itself, and multiplication 
00356     with double and 
00357     \ref NumericTraits "NumericTraits" must 
00358     be defined. In addition, this type must be less-comparable.
00359     
00360     <b> Declarations:</b>
00361     
00362     pass arguments explicitly:
00363     \code
00364     namespace vigra {
00365         template <class SrcIterator, class SrcAccessor,
00366               class DestIterator, class DestAccessor,
00367               class GradValue,
00368               class DestValue = DestAccessor::value_type>
00369         void differenceOfExponentialCrackEdgeImage(
00370                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00371                DestIterator dul, DestAccessor da,
00372                double scale, GradValue gradient_threshold, 
00373                DestValue edge_marker = NumericTraits<DestValue>::one())
00374     }
00375     \endcode
00376     
00377     use argument objects in conjunction with \ref ArgumentObjectFactories:
00378     \code
00379     namespace vigra {
00380         template <class SrcIterator, class SrcAccessor,
00381               class DestIterator, class DestAccessor, 
00382               class GradValue,
00383               class DestValue = DestAccessor::value_type>
00384         inline 
00385         void differenceOfExponentialCrackEdgeImage(
00386                triple<SrcIterator, SrcIterator, SrcAccessor> src,
00387                pair<DestIterator, DestAccessor> dest,
00388                double scale, GradValue gradient_threshold,
00389                DestValue edge_marker = NumericTraits<DestValue>::one())
00390     }
00391     \endcode
00392     
00393     <b> Usage:</b>
00394     
00395         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00396     Namespace: vigra
00397     
00398     \code
00399     vigra::BImage src(w,h), edges(2*w-1,2*h-1);
00400     
00401     // empty edge image
00402     edges = 0;
00403     ...
00404     
00405     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00406     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 
00407                                      0.8, 4.0, 1);
00408     \endcode
00409 
00410     <b> Required Interface:</b>
00411     
00412     \code
00413     SrcImageIterator src_upperleft, src_lowerright;
00414     DestImageIterator dest_upperleft;
00415     
00416     SrcAccessor src_accessor;
00417     DestAccessor dest_accessor;
00418     
00419     SrcAccessor::value_type u = src_accessor(src_upperleft);
00420     double d;
00421     GradValue gradient_threshold;
00422     
00423     u = u + u
00424     u = u - u
00425     u = u * u
00426     u = d * u
00427     u < gradient_threshold
00428     
00429     DestValue edge_marker;
00430     dest_accessor.set(edge_marker, dest_upperleft);
00431     \endcode
00432     
00433     <b> Preconditions:</b>
00434     
00435     \code
00436     scale > 0
00437     gradient_threshold > 0
00438     \endcode
00439     
00440     The destination image must have twice the size of the source:
00441     \code
00442     w_dest = 2 * w_src - 1
00443     h_dest = 2 * h_src - 1
00444     \endcode
00445 */
00446 template <class SrcIterator, class SrcAccessor,
00447           class DestIterator, class DestAccessor,
00448           class GradValue, class DestValue>
00449 void differenceOfExponentialCrackEdgeImage(
00450                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00451            DestIterator dul, DestAccessor da,
00452            double scale, GradValue gradient_threshold, 
00453            DestValue edge_marker)
00454 {
00455     vigra_precondition(scale > 0,
00456                  "differenceOfExponentialCrackEdgeImage(): scale > 0 required.");
00457          
00458     vigra_precondition(gradient_threshold > 0,
00459                  "differenceOfExponentialCrackEdgeImage(): "
00460          "gradient_threshold > 0 required.");
00461          
00462     int w = slr.x - sul.x;
00463     int h = slr.y - sul.y;
00464     int x, y;
00465 
00466     typedef typename 
00467         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00468     TMPTYPE;
00469     typedef BasicImage<TMPTYPE> TMPIMG;
00470 
00471     TMPIMG tmp(w,h);
00472     TMPIMG smooth(w,h);
00473     
00474     TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
00475     
00476     static const Diff2D right(1,0);
00477     static const Diff2D bottom(0,1);
00478     static const Diff2D left(-1,0);
00479     static const Diff2D top(0,-1);
00480     
00481     recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
00482     recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
00483 
00484     recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
00485     recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
00486 
00487     typename TMPIMG::Iterator iy = smooth.upperLeft();
00488     typename TMPIMG::Iterator ty = tmp.upperLeft();
00489     DestIterator              dy = dul;
00490     
00491     TMPTYPE thresh = (gradient_threshold * gradient_threshold) * 
00492                      NumericTraits<TMPTYPE>::one();
00493 
00494     // find zero crossings above threshold
00495     for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, dy.y+=2)
00496     {
00497         typename TMPIMG::Iterator ix = iy;
00498         typename TMPIMG::Iterator tx = ty;
00499         DestIterator              dx = dy;
00500 
00501         for(int x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
00502         {
00503             TMPTYPE diff = *tx - *ix;
00504             TMPTYPE gx = tx[right] - *tx;
00505             TMPTYPE gy = tx[bottom] - *tx;
00506 
00507             if((gx * gx > thresh) &&
00508                (diff * (tx[right] - ix[right]) < zero))
00509             {
00510                 da.set(edge_marker, dx, right);
00511             }
00512             if((gy * gy > thresh) &&
00513                (diff * (tx[bottom] - ix[bottom]) < zero))
00514             {
00515                 da.set(edge_marker, dx, bottom);
00516             }
00517         }
00518 
00519         TMPTYPE diff = *tx - *ix;
00520         TMPTYPE gy = tx[bottom] - *tx;
00521 
00522         if((gy * gy > thresh) &&
00523            (diff * (tx[bottom] - ix[bottom]) < zero))
00524         {
00525             da.set(edge_marker, dx, bottom);
00526         }
00527     }
00528     
00529     typename TMPIMG::Iterator ix = iy;
00530     typename TMPIMG::Iterator tx = ty;
00531     DestIterator              dx = dy;
00532     
00533     for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
00534     {
00535         TMPTYPE diff = *tx - *ix;
00536         TMPTYPE gx = tx[right] - *tx;
00537 
00538         if((gx * gx > thresh) &&
00539            (diff * (tx[right] - ix[right]) < zero))
00540         {
00541             da.set(edge_marker, dx, right);
00542         }
00543     }
00544 
00545     iy = smooth.upperLeft() + Diff2D(0,1);
00546     ty = tmp.upperLeft() + Diff2D(0,1);
00547     dy = dul + Diff2D(1,2);
00548     
00549     static const Diff2D topleft(-1,-1);
00550     static const Diff2D topright(1,-1);
00551     static const Diff2D bottomleft(-1,1);
00552     static const Diff2D bottomright(1,1);
00553 
00554     // find missing 1-cells below threshold (x-direction)
00555     for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
00556     {
00557         typename TMPIMG::Iterator ix = iy;
00558         typename TMPIMG::Iterator tx = ty;
00559         DestIterator              dx = dy;
00560 
00561         for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
00562         {
00563             if(da(dx) == edge_marker) continue;
00564 
00565             TMPTYPE diff = *tx - *ix;
00566 
00567             if((diff * (tx[right] - ix[right]) < zero) &&
00568                (((da(dx, bottomright) == edge_marker) && 
00569                  (da(dx, topleft) == edge_marker)) ||
00570                 ((da(dx, bottomleft) == edge_marker) && 
00571                  (da(dx, topright) == edge_marker))))
00572 
00573             {
00574                 da.set(edge_marker, dx);
00575             }
00576         }
00577     }
00578     
00579     iy = smooth.upperLeft() + Diff2D(1,0);
00580     ty = tmp.upperLeft() + Diff2D(1,0);
00581     dy = dul + Diff2D(2,1);
00582 
00583     // find missing 1-cells below threshold (y-direction)
00584     for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
00585     {
00586         typename TMPIMG::Iterator ix = iy;
00587         typename TMPIMG::Iterator tx = ty;
00588         DestIterator              dx = dy;
00589 
00590         for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
00591         {
00592             if(da(dx) == edge_marker) continue;
00593 
00594             TMPTYPE diff = *tx - *ix;
00595 
00596             if((diff * (tx[bottom] - ix[bottom]) < zero) &&
00597                (((da(dx, bottomright) == edge_marker) && 
00598                  (da(dx, topleft) == edge_marker)) ||
00599                 ((da(dx, bottomleft) == edge_marker) && 
00600                  (da(dx, topright) == edge_marker))))
00601 
00602             {
00603                 da.set(edge_marker, dx);
00604             }
00605         }
00606     }
00607     
00608     dy = dul + Diff2D(1,1);
00609 
00610     // find missing 0-cells 
00611     for(y=0; y<h-1; ++y, dy.y+=2)
00612     {
00613         DestIterator              dx = dy;
00614 
00615         for(int x=0; x<w-1; ++x, dx.x+=2)
00616         {
00617             static const Diff2D dist[] = {right, top, left, bottom };
00618 
00619             int i;
00620             for(i=0; i<4; ++i)
00621             {
00622             if(da(dx, dist[i]) == edge_marker) break;
00623             }
00624 
00625             if(i < 4) da.set(edge_marker, dx);
00626         }
00627     }
00628 }
00629 
00630 template <class SrcIterator, class SrcAccessor,
00631           class DestIterator, class DestAccessor,
00632           class GradValue, class DestValue>
00633 inline 
00634 void differenceOfExponentialCrackEdgeImage(
00635            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00636        pair<DestIterator, DestAccessor> dest,
00637        double scale, GradValue gradient_threshold,
00638        DestValue edge_marker)
00639 {
00640     differenceOfExponentialCrackEdgeImage(src.first, src.second, src.third,
00641                                         dest.first, dest.second,
00642                     scale, gradient_threshold,
00643                     edge_marker);
00644 }
00645 
00646 /********************************************************/
00647 /*                                                      */
00648 /*                  removeShortEdges                    */
00649 /*                                                      */
00650 /********************************************************/
00651 
00652 /** \brief Remove short edges from an edge image.
00653 
00654     This algorithm can be applied as a post-processing operation of 
00655     \ref differenceOfExponentialEdgeImage() and \ref differenceOfExponentialCrackEdgeImage(). 
00656     It removes all edges that are shorter than <TT>min_edge_length</TT>. The corresponding
00657     pixels are set to the <TT>non_edge_marker</TT>. The idea behind this algorithms is
00658     that very short edges are probably caused by noise and don't represent interesting
00659     image structure. Technically, the algorithms executes a connected components labeling,
00660     so the image's value type must be equality comparable. 
00661     
00662     If the source image fulfills the requirements of a \ref CrackEdgeImage,
00663     it will still do so after application of this algorithm.
00664     
00665     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 
00666     i.e. on only one image. Also, the algorithm assumes that all non-edges pixels are already
00667     marked with the given <TT>non_edge_marker</TT> value.
00668     
00669     <b> Declarations:</b>
00670     
00671     pass arguments explicitly:
00672     \code
00673     namespace vigra {
00674         template <class Iterator, class Accessor, class SrcValue>
00675         void removeShortEdges(
00676                Iterator sul, Iterator slr, Accessor sa,
00677                int min_edge_length, SrcValue non_edge_marker)
00678     }
00679     \endcode
00680     
00681     use argument objects in conjunction with \ref ArgumentObjectFactories:
00682     \code
00683     namespace vigra {
00684         template <class Iterator, class Accessor, class SrcValue>
00685         inline 
00686         void removeShortEdges(
00687                triple<Iterator, Iterator, Accessor> src,
00688                int min_edge_length, SrcValue non_edge_marker)
00689     }
00690     \endcode
00691     
00692     <b> Usage:</b>
00693     
00694         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00695     Namespace: vigra
00696     
00697     \code
00698     vigra::BImage src(w,h), edges(w,h);
00699     
00700     // empty edge image
00701     edges = 0;
00702     ...
00703     
00704     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00705     vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 
00706                                      0.8, 4.0, 1);
00707                     
00708     // zero edges shorter than 10 pixels
00709     vigra::removeShortEdges(srcImageRange(edges), 10, 0);
00710     \endcode
00711 
00712     <b> Required Interface:</b>
00713     
00714     \code
00715     SrcImageIterator src_upperleft, src_lowerright;
00716     DestImageIterator dest_upperleft;
00717     
00718     SrcAccessor src_accessor;
00719     DestAccessor dest_accessor;
00720     
00721     SrcAccessor::value_type u = src_accessor(src_upperleft);
00722 
00723     u == u
00724     
00725     SrcValue non_edge_marker;
00726     src_accessor.set(non_edge_marker, src_upperleft);
00727     \endcode
00728 */
00729 template <class Iterator, class Accessor, class Value>
00730 void removeShortEdges(
00731                Iterator sul, Iterator slr, Accessor sa,
00732            unsigned int min_edge_length, Value non_edge_marker)
00733 {
00734     int w = slr.x - sul.x;
00735     int h = slr.y - sul.y;
00736     int x,y;
00737 
00738     IImage labels(w, h);
00739     labels = 0;
00740 
00741     int number_of_regions = 
00742                 labelImageWithBackground(srcIterRange(sul,slr,sa), 
00743                                      destImage(labels), true, non_edge_marker);
00744     
00745     ArrayOfRegionStatistics<FindROISize<int> > 
00746                                          region_stats(number_of_regions);
00747     
00748     inspectTwoImages(srcImageRange(labels), srcImage(labels), region_stats);
00749              
00750     IImage::Iterator ly = labels.upperLeft();
00751     Iterator oy = sul;
00752     
00753     for(y=0; y<h; ++y, ++oy.y, ++ly.y)
00754     {
00755         Iterator ox(oy);
00756         IImage::Iterator lx(ly);
00757 
00758         for(x=0; x<w; ++x, ++ox.x, ++lx.x)
00759         {
00760             if(sa(ox) == non_edge_marker) continue;
00761             if((region_stats[*lx].count) < min_edge_length)
00762             {
00763                  sa.set(non_edge_marker, ox);
00764             }
00765         }
00766     }
00767 }
00768 
00769 template <class Iterator, class Accessor, class Value>
00770 inline 
00771 void removeShortEdges(
00772            triple<Iterator, Iterator, Accessor> src,
00773        unsigned int min_edge_length, Value non_edge_marker)
00774 {
00775     removeShortEdges(src.first, src.second, src.third,
00776                      min_edge_length, non_edge_marker);
00777 }
00778 
00779 /********************************************************/
00780 /*                                                      */
00781 /*             closeGapsInCrackEdgeImage                */
00782 /*                                                      */
00783 /********************************************************/
00784 
00785 /** \brief Close one-pixel wide gaps in a cell grid edge image.
00786 
00787     This algorithm is typically applied as a post-processing operation of 
00788     \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
00789     the requirements of a \ref CrackEdgeImage, and will still do so after 
00790     application of this algorithm.
00791 
00792     It closes one pixel wide gaps in the edges resulting from this algorithm. 
00793     Since these gaps are usually caused by zero crossing slightly below the gradient 
00794     threshold used in edge detection, this algorithms acts like a weak hysteresis 
00795     thresholding. The newly found edge pixels are marked with the given <TT>edge_marker</TT>.
00796     The image's value type must be equality comparable. 
00797     
00798     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 
00799     i.e. on only one image.
00800     
00801     <b> Declarations:</b>
00802     
00803     pass arguments explicitly:
00804     \code
00805     namespace vigra {
00806         template <class SrcIterator, class SrcAccessor, class SrcValue>
00807         void closeGapsInCrackEdgeImage(
00808                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
00809                SrcValue edge_marker)
00810     }
00811     \endcode
00812     
00813     use argument objects in conjunction with \ref ArgumentObjectFactories:
00814     \code
00815     namespace vigra {
00816         template <class SrcIterator, class SrcAccessor, class SrcValue>
00817         inline
00818         void closeGapsInCrackEdgeImage(
00819                triple<SrcIterator, SrcIterator, SrcAccessor> src,
00820                SrcValue edge_marker)
00821     }
00822     \endcode
00823     
00824     <b> Usage:</b>
00825     
00826         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00827     Namespace: vigra
00828     
00829     \code
00830     vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
00831     
00832     // empty edge image
00833     edges = 0;
00834     ...
00835     
00836     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00837     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 
00838                                          0.8, 4.0, 1);
00839                     
00840     // close gaps, mark with 1
00841     vigra::closeGapsInCrackEdgeImage(srcImageRange(edges), 1);
00842                     
00843     // zero edges shorter than 20 pixels
00844     vigra::removeShortEdges(srcImageRange(edges), 10, 0);
00845     \endcode
00846 
00847     <b> Required Interface:</b>
00848     
00849     \code
00850     SrcImageIterator src_upperleft, src_lowerright;
00851     
00852     SrcAccessor src_accessor;
00853     DestAccessor dest_accessor;
00854     
00855     SrcAccessor::value_type u = src_accessor(src_upperleft);
00856 
00857     u == u
00858     u != u
00859     
00860     SrcValue edge_marker;
00861     src_accessor.set(edge_marker, src_upperleft);
00862     \endcode
00863 */
00864 template <class SrcIterator, class SrcAccessor, class SrcValue>
00865 void closeGapsInCrackEdgeImage(
00866                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
00867            SrcValue edge_marker)
00868 {
00869     int w = (slr.x - sul.x) / 2;
00870     int h = (slr.y - sul.y) / 2;
00871     int x, y;
00872 
00873     int count1, count2, count3;
00874 
00875     static const Diff2D right(1,0);
00876     static const Diff2D bottom(0,1);
00877     static const Diff2D left(-1,0);
00878     static const Diff2D top(0,-1);
00879     
00880     static const Diff2D leftdist[] = { 
00881         Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
00882     static const Diff2D rightdist[] = { 
00883         Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
00884     static const Diff2D topdist[] = { 
00885         Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
00886     static const Diff2D bottomdist[] = { 
00887         Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};
00888 
00889     int i;
00890 
00891     SrcIterator sy = sul + Diff2D(0,1);
00892     SrcIterator sx;
00893     
00894     // close 1-pixel wide gaps (x-direction)
00895     for(y=0; y<h; ++y, sy.y+=2)
00896     {
00897         sx = sy + Diff2D(2,0);
00898 
00899         for(x=2; x<w; ++x, sx.x+=2)
00900         {
00901             if(sa(sx) == edge_marker) continue;
00902 
00903             if(sa(sx, left) != edge_marker) continue;
00904             if(sa(sx, right) != edge_marker) continue;
00905 
00906             count1 = 0;
00907             count2 = 0;
00908             count3 = 0;
00909 
00910             for(i=0; i<4; ++i)
00911             {
00912                 if(sa(sx, leftdist[i]) == edge_marker) 
00913                 {
00914                     ++count1;
00915                     count3 ^= 1 << i;
00916                 }
00917                 if(sa(sx, rightdist[i]) == edge_marker) 
00918                 {
00919                     ++count2;
00920                     count3 ^= 1 << i;
00921                 }
00922             }
00923 
00924             if(count1 <= 1 || count2 <= 1 || count3 == 15) 
00925             {
00926                 sa.set(edge_marker, sx);
00927             }
00928         }
00929    }
00930     
00931     sy = sul + Diff2D(1,2);
00932 
00933     // close 1-pixel wide gaps (y-direction)
00934     for(y=2; y<h; ++y, sy.y+=2)
00935     {
00936         sx = sy;
00937 
00938         for(x=0; x<w; ++x, sx.x+=2)
00939         {
00940             if(sa(sx) == edge_marker) continue;
00941 
00942             if(sa(sx, top) != edge_marker) continue;
00943             if(sa(sx, bottom) != edge_marker) continue;
00944 
00945             count1 = 0;
00946             count2 = 0;
00947             count3 = 0;
00948 
00949             for(i=0; i<4; ++i)
00950             {
00951                 if(sa(sx, topdist[i]) == edge_marker)
00952                 {
00953                     ++count1;
00954                     count3 ^= 1 << i;
00955                 }
00956                 if(sa(sx, bottomdist[i]) == edge_marker)
00957                 {
00958                     ++count2;
00959                     count3 ^= 1 << i;
00960                 }
00961             }
00962 
00963             if(count1 <= 1 || count2 <= 1 || count3 == 15)
00964             {
00965                 sa.set(edge_marker, sx);
00966             }
00967         }
00968     }
00969 }
00970 
00971 template <class SrcIterator, class SrcAccessor, class SrcValue>
00972 inline
00973 void closeGapsInCrackEdgeImage(
00974            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00975        SrcValue edge_marker)
00976 {
00977     closeGapsInCrackEdgeImage(src.first, src.second, src.third,
00978                     edge_marker);
00979 }
00980 
00981 /********************************************************/
00982 /*                                                      */
00983 /*              beautifyCrackEdgeImage                  */
00984 /*                                                      */
00985 /********************************************************/
00986 
00987 /** \brief Beautify crack edge image for visualization.
00988 
00989     This algorithm is applied as a post-processing operation of 
00990     \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
00991     the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after 
00992     application of this algorithm. In particular, the algorithm removes zero-cells 
00993     marked as edges to avoid staircase effects on diagonal lines like this:
00994     
00995     \code
00996     original edge points (*)     resulting edge points
00997     
00998           . * . . .                   . * . . .
00999           . * * * .                   . . * . .
01000           . . . * .           =>      . . . * .
01001           . . . * *                   . . . . *
01002           . . . . .                   . . . . .
01003     \endcode
01004     
01005     Therfore, this algorithm should only be applied as a vizualization aid, i.e. 
01006     for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>, 
01007     and background pixels with <TT>background_marker</TT>. The image's value type must be 
01008     equality comparable. 
01009     
01010     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 
01011     i.e. on only one image.
01012     
01013     <b> Declarations:</b>
01014     
01015     pass arguments explicitly:
01016     \code
01017     namespace vigra {
01018         template <class SrcIterator, class SrcAccessor, class SrcValue>
01019         void beautifyCrackEdgeImage(
01020                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
01021                SrcValue edge_marker, SrcValue background_marker)
01022     }
01023     \endcode
01024     
01025     use argument objects in conjunction with \ref ArgumentObjectFactories:
01026     \code
01027     namespace vigra {
01028         template <class SrcIterator, class SrcAccessor, class SrcValue>
01029         inline
01030         void beautifyCrackEdgeImage(
01031                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01032                SrcValue edge_marker, SrcValue background_marker)
01033     }
01034     \endcode
01035     
01036     <b> Usage:</b>
01037     
01038         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01039     Namespace: vigra
01040     
01041     \code
01042     vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
01043     
01044     // empty edge image
01045     edges = 0;
01046     ...
01047     
01048     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
01049     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 
01050                                          0.8, 4.0, 1);
01051                     
01052     // beautify edge image for visualization
01053     vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
01054     
01055     // show to the user
01056     window.open(edges);
01057     \endcode
01058 
01059     <b> Required Interface:</b>
01060     
01061     \code
01062     SrcImageIterator src_upperleft, src_lowerright;
01063     
01064     SrcAccessor src_accessor;
01065     DestAccessor dest_accessor;
01066     
01067     SrcAccessor::value_type u = src_accessor(src_upperleft);
01068 
01069     u == u
01070     u != u
01071     
01072     SrcValue background_marker;
01073     src_accessor.set(background_marker, src_upperleft);
01074     \endcode
01075 */
01076 template <class SrcIterator, class SrcAccessor, class SrcValue>
01077 void beautifyCrackEdgeImage(
01078                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
01079            SrcValue edge_marker, SrcValue background_marker)
01080 {
01081     int w = (slr.x - sul.x) / 2;
01082     int h = (slr.y - sul.y) / 2;
01083     int x, y;
01084 
01085     SrcIterator sy = sul + Diff2D(1,1);
01086     SrcIterator sx;
01087     
01088     static const Diff2D right(1,0);
01089     static const Diff2D bottom(0,1);
01090     static const Diff2D left(-1,0);
01091     static const Diff2D top(0,-1);
01092     
01093     //  delete 0-cells at corners
01094     for(y=0; y<h; ++y, sy.y+=2)
01095     {
01096         sx = sy;
01097 
01098         for(x=0; x<w; ++x, sx.x+=2)
01099         {
01100             if(sa(sx) != edge_marker) continue;
01101 
01102             if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
01103             if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;
01104 
01105             sa.set(background_marker, sx);
01106         }
01107     }
01108 }
01109 
01110 template <class SrcIterator, class SrcAccessor, class SrcValue>
01111 inline
01112 void beautifyCrackEdgeImage(
01113            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01114            SrcValue edge_marker, SrcValue background_marker)
01115 {
01116     beautifyCrackEdgeImage(src.first, src.second, src.third,
01117                     edge_marker, background_marker);
01118 }
01119 
01120 
01121 /** Helper class that stores edgel attributes.
01122 */
01123 class Edgel
01124 {
01125   public:
01126         /** The edgel's sub-pixel x coordinate.
01127         */    
01128     float x; 
01129 
01130         /** The edgel's sub-pixel y coordinate.
01131         */    
01132     float y;
01133 
01134         /** The edgel's strength (magnitude of the gradient vector).
01135         */    
01136     float strength;
01137     
01138         /**
01139         The edgel's orientation. This is the angle 
01140         between the x-axis and the edge, so that the bright side of the 
01141         edge is on the right. The angle is measured
01142         counter-clockwise in radians like this:
01143 
01144 
01145         \code
01146 
01147   edgel axis
01148       \  (bright side)
01149  (dark \
01150  side)  \ /__ 
01151          \\  \ orientation angle
01152           \  |
01153            +------------> x-axis
01154            |
01155            |
01156            |
01157            |
01158     y-axis V
01159         \endcode
01160 
01161         So, for example a vertical edge with its dark side on the left
01162         has orientation PI/2, and a horizontal edge with dark side on top
01163         has orientation 0. Obviously, the edge's orientation changes
01164         by PI if the contrast is reversed.
01165 
01166         */
01167     float orientation;
01168     
01169     Edgel()
01170     : x(0.0f), y(0.0f), strength(0.0f), orientation(0.0f)
01171     {}
01172     
01173     Edgel(float ix, float iy, float is, float io)
01174     : x(ix), y(iy), strength(is), orientation(io)
01175     {}
01176 };
01177 
01178 template <class Image1, class Image2>
01179 void internalCannyFindEdgels(Image1 const & dx,
01180                              Image1 const & dy,
01181                              Image2 const & magnitude,
01182                              std::vector<Edgel> & edgels)
01183 {
01184     typedef typename Image1::value_type PixelType;
01185     
01186     PixelType zero = NumericTraits<PixelType>::zero();
01187     double tan22_5 = M_SQRT2 - 1.0;
01188     
01189     for(int y=1; y<dx.height()-1; ++y)
01190     {
01191         for(int x=1; x<dx.width()-1; ++x)
01192         {
01193             bool maximum_found = false;
01194             Edgel edgel;
01195             
01196             PixelType gradx = dx(x,y);
01197             PixelType grady = dy(x,y);
01198             
01199             // find out quadrant
01200             if(abs(grady) < tan22_5*abs(gradx))
01201             {
01202                 // north-south edge
01203                 PixelType m1 = magnitude(x-1, y);
01204                 PixelType m2 = magnitude(x, y);
01205                 PixelType m3 = magnitude(x+1, y);
01206                 
01207                 if(m1 < m2 && m3 <= m2)
01208                 {
01209                     edgel.y = y;
01210                 
01211                     // local maximum => quadratic interpolation of sub-pixel location
01212                     PixelType del = (m1 - m3) / 2.0;
01213                     del /= (m1 + m3 - 2.0*m2);
01214                     edgel.x = x + del;
01215                     edgel.strength = m2;
01216                     
01217                     maximum_found = true;                    
01218                 }
01219             }
01220             else if(abs(gradx) < tan22_5*abs(grady))
01221             {
01222                 // west-east edge
01223                 PixelType m1 = magnitude(x, y-1);
01224                 PixelType m2 = magnitude(x, y);
01225                 PixelType m3 = magnitude(x, y+1);
01226                 
01227                 if(m1 < m2 && m3 <= m2)
01228                 {
01229                     edgel.x = x;
01230                 
01231                     // local maximum => quadratic interpolation of sub-pixel location
01232                     PixelType del = (m1 - m3) / 2.0;
01233                     del /= (m1 + m3 - 2.0*m2);
01234                     edgel.y = y + del;
01235                     edgel.strength = m2;
01236                     
01237                     maximum_found = true;                    
01238                 }
01239             }
01240             else if(gradx*grady < zero)
01241             {
01242                 // north-west-south-east edge
01243                 PixelType m1 = magnitude(x+1, y-1);
01244                 PixelType m2 = magnitude(x, y);
01245                 PixelType m3 = magnitude(x-1, y+1);
01246                 
01247                 if(m1 < m2 && m3 <= m2)
01248                 {
01249                     // local maximum => quadratic interpolation of sub-pixel location
01250                     PixelType del = (m1 - m3) / 2.0;
01251                     del /= (m1 + m3 - 2.0*m2);
01252                     edgel.x = x - del;
01253                     edgel.y = y + del;
01254                     edgel.strength = m2;
01255                     
01256                     maximum_found = true;                    
01257                 }
01258             }
01259             else
01260             {
01261                 // north-east-south-west edge
01262                 PixelType m1 = magnitude(x-1, y-1);
01263                 PixelType m2 = magnitude(x, y);
01264                 PixelType m3 = magnitude(x+1, y+1);
01265                 
01266                 if(m1 < m2 && m3 <= m2)
01267                 {
01268                     // local maximum => quadratic interpolation of sub-pixel location
01269                     PixelType del = (m1 - m3) / 2.0;
01270                     del /= (m1 + m3 - 2.0*m2);
01271                     edgel.x = x + del;
01272                     edgel.y = y + del;
01273                     edgel.strength = m2;
01274                     
01275                     maximum_found = true;                    
01276                 }
01277             }
01278             
01279             if(maximum_found)
01280             {
01281                 double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
01282                 if(orientation < 0.0)
01283                     orientation += 2.0*M_PI;
01284                 edgel.orientation = orientation;
01285                 edgels.push_back(edgel);
01286             }
01287         }
01288     }
01289 }
01290 
01291 /********************************************************/
01292 /*                                                      */
01293 /*                      cannyEdgelList                  */
01294 /*                                                      */
01295 /********************************************************/
01296 
01297 /** \brief Simple implementation of Canny's edge detector.
01298 
01299     This operator first calculates the gradient vector for each
01300     pixel of the image using first derivatives of a Gaussian at the 
01301     given scale. Then a very simple non-maxima supression is performed: 
01302     for each 3x3 neighborhood, it is determined whether the center pixel has 
01303     larger gradient magnitude than its two neighbors in gradient direction
01304     (where the direction is rounded into octands). If this is the case,
01305     a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
01306     edgel position is determined by fitting a parabola 
01307     to the three gradient magnitude values 
01308     mentioned above. The sub-pixel location of the parabola's tip 
01309     and the gradient magnitude and direction are written in the newly created edgel.
01310     
01311     <b> Declarations:</b>
01312     
01313     pass arguments explicitly:
01314     \code
01315     namespace vigra {
01316         template <class SrcIterator, class SrcAccessor>
01317         void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01318                                 std::vector<Edgel> & edgels, double scale);
01319     }
01320     \endcode
01321     
01322     use argument objects in conjunction with \ref ArgumentObjectFactories:
01323     \code
01324     namespace vigra {
01325         template <class SrcIterator, class SrcAccessor>
01326         void 
01327         cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01328                        std::vector<Edgel> & edgels, double scale);
01329     }
01330     \endcode
01331     
01332     <b> Usage:</b>
01333     
01334     <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01335     Namespace: vigra
01336     
01337     \code
01338     vigra::BImage src(w,h);
01339     
01340     // empty edgel list
01341     std::vector<vigra::Edgel> edgels;
01342     ...
01343     
01344     // find edgels at scale 0.8  
01345     vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
01346     \endcode
01347 
01348     <b> Required Interface:</b>
01349     
01350     \code
01351     SrcImageIterator src_upperleft;
01352     SrcAccessor src_accessor;
01353     
01354     src_accessor(src_upperleft);
01355     \endcode
01356     
01357     SrcAccessor::value_type must be a type convertible to float
01358     
01359     <b> Preconditions:</b>
01360     
01361     \code
01362     scale > 0
01363     \endcode
01364 */
01365 template <class SrcIterator, class SrcAccessor>
01366 void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01367                         std::vector<Edgel> & edgels, double scale)
01368 {
01369     int w = lr.x - ul.x;
01370     int h = lr.y - ul.y;
01371     
01372     // calculate image gradients
01373     typedef typename 
01374         NumericTraits<typename SrcAccessor::value_type>::RealPromote
01375         TmpType;
01376         
01377     BasicImage<TmpType> tmp(w,h), dx(w,h), dy(w,h);
01378     
01379     Kernel1D<double> smooth, grad;
01380     
01381     smooth.initGaussian(scale);
01382     grad.initGaussianDerivative(scale, 1);
01383     
01384     separableConvolveX(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
01385     separableConvolveY(srcImageRange(tmp), destImage(dx), kernel1d(smooth));
01386     
01387     separableConvolveY(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
01388     separableConvolveX(srcImageRange(tmp), destImage(dy), kernel1d(smooth));
01389     
01390     combineTwoImages(srcImageRange(dx), srcImage(dy), destImage(tmp),
01391                      MagnitudeFunctor<TmpType>());
01392     
01393     
01394     // find edgels
01395     internalCannyFindEdgels(dx, dy, tmp, edgels);
01396 }
01397 
01398 template <class SrcIterator, class SrcAccessor>
01399 inline void 
01400 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01401                std::vector<Edgel> & edgels, double scale)
01402 {
01403     cannyEdgelList(src.first, src.second, src.third, edgels, scale);
01404 }
01405 
01406 /********************************************************/
01407 /*                                                      */
01408 /*                       cannyEdgeImage                 */
01409 /*                                                      */
01410 /********************************************************/
01411 
01412 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
01413 
01414     This operator first calls \ref cannyEdgelList() to generate an 
01415     edgel list for the given image. Than it scans this list and selects edgels
01416     whose strength is above the given <TT>gradient_threshold</TT>. For each of these 
01417     edgels, the edgel's location is rounded to the nearest pixel, and that
01418     pixel marked with the given <TT>edge_marker</TT>.
01419     
01420     <b> Declarations:</b>
01421     
01422     pass arguments explicitly:
01423     \code
01424     namespace vigra {
01425         template <class SrcIterator, class SrcAccessor,
01426                   class DestIterator, class DestAccessor, 
01427                   class GradValue, class DestValue>
01428         void cannyEdgeImage(
01429                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01430                    DestIterator dul, DestAccessor da,
01431                    double scale, GradValue gradient_threshold, DestValue edge_marker);
01432     }
01433     \endcode
01434     
01435     use argument objects in conjunction with \ref ArgumentObjectFactories:
01436     \code
01437     namespace vigra {
01438         template <class SrcIterator, class SrcAccessor,
01439                   class DestIterator, class DestAccessor, 
01440                   class GradValue, class DestValue>
01441         inline void cannyEdgeImage(
01442                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01443                    pair<DestIterator, DestAccessor> dest,
01444                    double scale, GradValue gradient_threshold, DestValue edge_marker);
01445     }
01446     \endcode
01447     
01448     <b> Usage:</b>
01449     
01450     <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01451     Namespace: vigra
01452     
01453     \code
01454     vigra::BImage src(w,h), edges(w,h);
01455     
01456     // empty edge image
01457     edges = 0;
01458     ...
01459     
01460     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
01461     vigra::cannyEdgeImage(srcImageRange(src), destImage(edges), 
01462                                      0.8, 4.0, 1);
01463     \endcode
01464 
01465     <b> Required Interface:</b>
01466     
01467     see also: \ref cannyEdgelList().
01468     
01469     \code
01470     DestImageIterator dest_upperleft;
01471     DestAccessor dest_accessor;
01472     DestValue edge_marker;
01473     
01474     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
01475     \endcode
01476     
01477     <b> Preconditions:</b>
01478     
01479     \code
01480     scale > 0
01481     gradient_threshold > 0
01482     \endcode
01483 */
01484 template <class SrcIterator, class SrcAccessor,
01485           class DestIterator, class DestAccessor, 
01486           class GradValue, class DestValue>
01487 void cannyEdgeImage(
01488            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01489            DestIterator dul, DestAccessor da,
01490            double scale, GradValue gradient_threshold, DestValue edge_marker)
01491 {
01492     std::vector<Edgel> edgels;
01493     
01494     cannyEdgelList(sul, slr, sa, edgels, scale);
01495     
01496     for(unsigned int i=0; i<edgels.size(); ++i)
01497     {
01498         if(gradient_threshold < edgels[i].strength)
01499         {
01500             Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
01501             
01502             da.set(edge_marker, dul, pix);
01503         }
01504     }
01505 }
01506 
01507 template <class SrcIterator, class SrcAccessor,
01508           class DestIterator, class DestAccessor, 
01509           class GradValue, class DestValue>
01510 inline void cannyEdgeImage(
01511            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01512            pair<DestIterator, DestAccessor> dest,
01513            double scale, GradValue gradient_threshold, DestValue edge_marker)
01514 {
01515     cannyEdgeImage(src.first, src.second, src.third,
01516                    dest.first, dest.second,
01517                    scale, gradient_threshold, edge_marker);
01518 }
01519 
01520 //@}
01521 
01522 /** \page CrackEdgeImage Crack Edge Image
01523 
01524 Crack edges are marked <i>between</i> the pixels of an image. 
01525 A Crack Edge Image is an image that represents these edges. In order
01526 to accomodate the cracks, the Crack Edge Image must be twice as large
01527 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
01528 can easily be derived from a binary image or from the signs of the 
01529 response of a Laplacean filter. Consider the following sketch, where
01530 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
01531 <TT>*</TT> the resulting crack edges.
01532 
01533     \code
01534 sign of difference image         insert cracks         resulting CrackEdgeImage
01535 
01536                                    + . - . -              . * . . .
01537       + - -                        . . . . .              . * * * .
01538       + + -               =>       + . + . -      =>      . . . * .
01539       + + +                        . . . . .              . . . * *
01540                                    + . + . +              . . . . .
01541     \endcode
01542 
01543 Starting from the original binary image (left), we insert crack pixels
01544 to get to the double-sized image (center). Finally, we mark all 
01545 crack pixels whose non-crack neighbors have different signs as 
01546 crack edge points, while all other pixels (crack and non-crack) become 
01547 region pixels.
01548 
01549 <b>Requirements on a Crack Edge Image:</b>
01550 
01551 <ul>
01552     <li>Crack Edge Images have odd width and height.
01553     <li>Crack pixels have at least one odd coordinate.
01554     <li>Only crack pixels may be marked as edge points.
01555     <li>Crack pixels with two odd coordinates must be marked as edge points
01556         whenever any of their neighboring crack pixels was marked.  
01557 </ul>
01558 
01559 The last two requirements ensure that both edges and regions are 4-connected. 
01560 Thus, 4-connectivity and 8-connectivity yield identical connected 
01561 components in a Crack Edge Image (so called <i>well-composedness</i>).
01562 This ensures that Crack Edge Images have nice topological properties
01563 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000). 
01564 */
01565 
01566 
01567 } // namespace vigra
01568 
01569 #endif // VIGRA_EDGEDETECTION_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.2 (27 Jan 2005)