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

| vigra/edgedetection.hxx |
--- | --- | ---

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.5.0 (7 Dec 2006)