[ 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.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)