vigra - vigra/edgedetection.hxx Source File

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

vigra/edgedetection.hxx |

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)