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

# vigra/edgedetection.hxx

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
    00008 /*    The VIGRA Website is                                              */
    00009 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00010 /*    Please direct questions, bug reports, and contributions to        */
    00011 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00012 /*        vigra@informatik.uni-hamburg.de                               */
    00013 /*                                                                      */
    00014 /*    Permission is hereby granted, free of charge, to any person       */
    00015 /*    obtaining a copy of this software and associated documentation    */
    00016 /*    files (the "Software"), to deal in the Software without           */
    00017 /*    restriction, including without limitation the rights to use,      */
    00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00019 /*    sell copies of the Software, and to permit persons to whom the    */
    00020 /*    Software is furnished to do so, subject to the following          */
    00021 /*    conditions:                                                       */
    00022 /*                                                                      */
    00023 /*    The above copyright notice and this permission notice shall be    */
    00024 /*    included in all copies or substantial portions of the             */
    00025 /*    Software.                                                         */
    00026 /*                                                                      */
    00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00035 /*                                                                      */
    00036 /************************************************************************/
    00037
    00038
    00039 #ifndef VIGRA_EDGEDETECTION_HXX
    00040 #define VIGRA_EDGEDETECTION_HXX
    00041
    00042 #include <vector>
    00043 #include <queue>
    00044 #include <cmath>     // sqrt(), abs()
    00045 #include "utilities.hxx"
    00046 #include "numerictraits.hxx"
    00047 #include "stdimage.hxx"
    00048 #include "stdimagefunctions.hxx"
    00049 #include "recursiveconvolution.hxx"
    00050 #include "separableconvolution.hxx"
    00051 #include "labelimage.hxx"
    00052 #include "mathutil.hxx"
    00053 #include "pixelneighborhood.hxx"
    00054 #include "linear_solve.hxx"
    00055
    00056
    00057 namespace vigra {
    00058
    00059 /** \addtogroup EdgeDetection Edge Detection
    00060     Edge detectors based on first and second derivatives,
    00061           and related post-processing.
    00062 */
    00063 //@{
    00064
    00065 /********************************************************/
    00066 /*                                                      */
    00067 /*           differenceOfExponentialEdgeImage           */
    00068 /*                                                      */
    00069 /********************************************************/
    00070
    00071 /** \brief Detect and mark edges in an edge image using the Shen/Castan zero-crossing detector.
    00072
    00073     This operator applies an exponential filter to the source image
    00074     at the given <TT>scale</TT> and subtracts the result from the original image.
    00075     Zero crossings are detected in the resulting difference image. Whenever the
    00076     gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
    00077     an edge point is marked (using <TT>edge_marker</TT>) in the destination image on
    00078     the darker side of the zero crossing (note that zero crossings occur
    00079     <i>between</i> pixels). For example:
    00080
    00081     \code
    00082     sign of difference image     resulting edge points (*)
    00083
    00084         + - -                          * * .
    00085         + + -               =>         . * *
    00086         + + +                          . . .
    00087     \endcode
    00088
    00089     Non-edge pixels (<TT>.</TT>) remain untouched in the destination image.
    00090     The result can be improved by the post-processing operation \ref removeShortEdges().
    00091     A more accurate edge placement can be achieved with the function
    00092     \ref differenceOfExponentialCrackEdgeImage().
    00093
    00094     The source value type
    00095     (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition,
    00096     subtraction and multiplication of the type with itself, and multiplication
    00097     with double and
    00098     \ref NumericTraits "NumericTraits" must
    00099     be defined. In addition, this type must be less-comparable.
    00100
    00101     <b> Declarations:</b>
    00102
    00103     pass arguments explicitly:
    00104     \code
    00105     namespace vigra {
    00106         template <class SrcIterator, class SrcAccessor,
    00107               class DestIterator, class DestAccessor,
    00108               class GradValue,
    00109               class DestValue = DestAccessor::value_type>
    00110         void differenceOfExponentialEdgeImage(
    00111                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    00112                DestIterator dul, DestAccessor da,
    00113                double scale, GradValue gradient_threshold,
    00114                DestValue edge_marker = NumericTraits<DestValue>::one())
    00115     }
    00116     \endcode
    00117
    00118     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00119     \code
    00120     namespace vigra {
    00121         template <class SrcIterator, class SrcAccessor,
    00122               class DestIterator, class DestAccessor,
    00123               class GradValue,
    00124               class DestValue = DestAccessor::value_type>
    00125         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> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><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 = (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 conformes 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> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    00415     Namespace: vigra
    00416
    00417     \code
    00418     vigra::BImage src(w,h), edges(2*w-1,2*h-1);
    00419
    00420     // empty edge image
    00421     edges = 0;
    00422     ...
    00423
    00424     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
    00425     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
    00426                                      0.8, 4.0, 1);
    00427     \endcode
    00428
    00429     <b> Required Interface:</b>
    00430
    00431     \code
    00432     SrcImageIterator src_upperleft, src_lowerright;
    00433     DestImageIterator dest_upperleft;
    00434
    00435     SrcAccessor src_accessor;
    00436     DestAccessor dest_accessor;
    00437
    00438     SrcAccessor::value_type u = src_accessor(src_upperleft);
    00439     double d;
    00440     GradValue gradient_threshold;
    00441
    00442     u = u + u
    00443     u = u - u
    00444     u = u * u
    00445     u = d * u
    00446     u < gradient_threshold
    00447
    00448     DestValue edge_marker;
    00449     dest_accessor.set(edge_marker, dest_upperleft);
    00450     \endcode
    00451
    00452     <b> Preconditions:</b>
    00453
    00454     \code
    00455     scale > 0
    00456     gradient_threshold > 0
    00457     \endcode
    00458
    00459     The destination image must have twice the size of the source:
    00460     \code
    00461     w_dest = 2 * w_src - 1
    00462     h_dest = 2 * h_src - 1
    00463     \endcode
    00464 */
    00465 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 = (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> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><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> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><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) / 2;
    00893     int h = (slr.y - sul.y) / 2;
    00894     int x, y;
    00895
    00896     int count1, count2, count3;
    00897
    00898     static const Diff2D right(1,0);
    00899     static const Diff2D bottom(0,1);
    00900     static const Diff2D left(-1,0);
    00901     static const Diff2D top(0,-1);
    00902
    00903     static const Diff2D leftdist[] = {
    00904         Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
    00905     static const Diff2D rightdist[] = {
    00906         Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
    00907     static const Diff2D topdist[] = {
    00908         Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
    00909     static const Diff2D bottomdist[] = {
    00910         Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};
    00911
    00912     int i;
    00913
    00914     SrcIterator sy = sul + Diff2D(0,1);
    00915     SrcIterator sx;
    00916
    00917     // close 1-pixel wide gaps (x-direction)
    00918     for(y=0; y<h; ++y, sy.y+=2)
    00919     {
    00920         sx = sy + Diff2D(2,0);
    00921
    00922         for(x=2; x<w; ++x, sx.x+=2)
    00923         {
    00924             if(sa(sx) == edge_marker) continue;
    00925
    00926             if(sa(sx, left) != edge_marker) continue;
    00927             if(sa(sx, right) != edge_marker) continue;
    00928
    00929             count1 = 0;
    00930             count2 = 0;
    00931             count3 = 0;
    00932
    00933             for(i=0; i<4; ++i)
    00934             {
    00935                 if(sa(sx, leftdist[i]) == edge_marker)
    00936                 {
    00937                     ++count1;
    00938                     count3 ^= 1 << i;
    00939                 }
    00940                 if(sa(sx, rightdist[i]) == edge_marker)
    00941                 {
    00942                     ++count2;
    00943                     count3 ^= 1 << i;
    00944                 }
    00945             }
    00946
    00947             if(count1 <= 1 || count2 <= 1 || count3 == 15)
    00948             {
    00949                 sa.set(edge_marker, sx);
    00950             }
    00951         }
    00952    }
    00953
    00954     sy = sul + Diff2D(1,2);
    00955
    00956     // close 1-pixel wide gaps (y-direction)
    00957     for(y=2; y<h; ++y, sy.y+=2)
    00958     {
    00959         sx = sy;
    00960
    00961         for(x=0; x<w; ++x, sx.x+=2)
    00962         {
    00963             if(sa(sx) == edge_marker) continue;
    00964
    00965             if(sa(sx, top) != edge_marker) continue;
    00966             if(sa(sx, bottom) != edge_marker) continue;
    00967
    00968             count1 = 0;
    00969             count2 = 0;
    00970             count3 = 0;
    00971
    00972             for(i=0; i<4; ++i)
    00973             {
    00974                 if(sa(sx, topdist[i]) == edge_marker)
    00975                 {
    00976                     ++count1;
    00977                     count3 ^= 1 << i;
    00978                 }
    00979                 if(sa(sx, bottomdist[i]) == edge_marker)
    00980                 {
    00981                     ++count2;
    00982                     count3 ^= 1 << i;
    00983                 }
    00984             }
    00985
    00986             if(count1 <= 1 || count2 <= 1 || count3 == 15)
    00987             {
    00988                 sa.set(edge_marker, sx);
    00989             }
    00990         }
    00991     }
    00992 }
    00993
    00994 template <class SrcIterator, class SrcAccessor, class SrcValue>
    00995 inline
    00996 void closeGapsInCrackEdgeImage(
    00997            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00998        SrcValue edge_marker)
    00999 {
    01000     closeGapsInCrackEdgeImage(src.first, src.second, src.third,
    01001                     edge_marker);
    01002 }
    01003
    01004 /********************************************************/
    01005 /*                                                      */
    01006 /*              beautifyCrackEdgeImage                  */
    01007 /*                                                      */
    01008 /********************************************************/
    01009
    01010 /** \brief Beautify crack edge image for visualization.
    01011
    01012     This algorithm is applied as a post-processing operation of
    01013     \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
    01014     the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after
    01015     application of this algorithm. In particular, the algorithm removes zero-cells
    01016     marked as edges to avoid staircase effects on diagonal lines like this:
    01017
    01018     \code
    01019     original edge points (*)     resulting edge points
    01020
    01021           . * . . .                   . * . . .
    01022           . * * * .                   . . * . .
    01023           . . . * .           =>      . . . * .
    01024           . . . * *                   . . . . *
    01025           . . . . .                   . . . . .
    01026     \endcode
    01027
    01028     Therfore, this algorithm should only be applied as a vizualization aid, i.e.
    01029     for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>,
    01030     and background pixels with <TT>background_marker</TT>. The image's value type must be
    01031     equality comparable.
    01032
    01033     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
    01034     i.e. on only one image.
    01035
    01036     <b> Declarations:</b>
    01037
    01038     pass arguments explicitly:
    01039     \code
    01040     namespace vigra {
    01041         template <class SrcIterator, class SrcAccessor, class SrcValue>
    01042         void beautifyCrackEdgeImage(
    01043                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01044                SrcValue edge_marker, SrcValue background_marker)
    01045     }
    01046     \endcode
    01047
    01048     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01049     \code
    01050     namespace vigra {
    01051         template <class SrcIterator, class SrcAccessor, class SrcValue>
    01052         void beautifyCrackEdgeImage(
    01053                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01054                SrcValue edge_marker, SrcValue background_marker)
    01055     }
    01056     \endcode
    01057
    01058     <b> Usage:</b>
    01059
    01060         <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    01061     Namespace: vigra
    01062
    01063     \code
    01064     vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
    01065
    01066     // empty edge image
    01067     edges = 0;
    01068     ...
    01069
    01070     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
    01071     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
    01072                                          0.8, 4.0, 1);
    01073
    01074     // beautify edge image for visualization
    01075     vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
    01076
    01077     // show to the user
    01078     window.open(edges);
    01079     \endcode
    01080
    01081     <b> Required Interface:</b>
    01082
    01083     \code
    01084     SrcImageIterator src_upperleft, src_lowerright;
    01085
    01086     SrcAccessor src_accessor;
    01087     DestAccessor dest_accessor;
    01088
    01089     SrcAccessor::value_type u = src_accessor(src_upperleft);
    01090
    01091     u == u
    01092     u != u
    01093
    01094     SrcValue background_marker;
    01095     src_accessor.set(background_marker, src_upperleft);
    01096     \endcode
    01097 */
    01098 doxygen_overloaded_function(template <...> void beautifyCrackEdgeImage)
    01099
    01100 template <class SrcIterator, class SrcAccessor, class SrcValue>
    01101 void beautifyCrackEdgeImage(
    01102                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01103            SrcValue edge_marker, SrcValue background_marker)
    01104 {
    01105     int w = (slr.x - sul.x) / 2;
    01106     int h = (slr.y - sul.y) / 2;
    01107     int x, y;
    01108
    01109     SrcIterator sy = sul + Diff2D(1,1);
    01110     SrcIterator sx;
    01111
    01112     static const Diff2D right(1,0);
    01113     static const Diff2D bottom(0,1);
    01114     static const Diff2D left(-1,0);
    01115     static const Diff2D top(0,-1);
    01116
    01117     //  delete 0-cells at corners
    01118     for(y=0; y<h; ++y, sy.y+=2)
    01119     {
    01120         sx = sy;
    01121
    01122         for(x=0; x<w; ++x, sx.x+=2)
    01123         {
    01124             if(sa(sx) != edge_marker) continue;
    01125
    01126             if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
    01127             if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;
    01128
    01129             sa.set(background_marker, sx);
    01130         }
    01131     }
    01132 }
    01133
    01134 template <class SrcIterator, class SrcAccessor, class SrcValue>
    01135 inline
    01136 void beautifyCrackEdgeImage(
    01137            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01138            SrcValue edge_marker, SrcValue background_marker)
    01139 {
    01140     beautifyCrackEdgeImage(src.first, src.second, src.third,
    01141                     edge_marker, background_marker);
    01142 }
    01143
    01144
    01145 /** Helper class that stores edgel attributes.
    01146 */
    01147 class Edgel
    01148 {
    01149   public:
    01150         /** The edgel's sub-pixel x coordinate.
    01151         */
    01152     float x;
    01153
    01154         /** The edgel's sub-pixel y coordinate.
    01155         */
    01156     float y;
    01157
    01158         /** The edgel's strength (magnitude of the gradient vector).
    01159         */
    01160     float strength;
    01161
    01162         /**
    01163         The edgel's orientation. This is the angle
    01164         between the x-axis and the edge, so that the bright side of the
    01165         edge is on the right. The angle is measured
    01166         counter-clockwise in radians like this:
    01167
    01168
    01169         \code
    01170
    01171   edgel axis
    01172       \  (bright side)
    01173  (dark \
    01174  side)  \ /__
    01175          \\  \ orientation angle
    01176           \  |
    01177            +------------> x-axis
    01178            |
    01179            |
    01180            |
    01181            |
    01182     y-axis V
    01183         \endcode
    01184
    01185         So, for example a vertical edge with its dark side on the left
    01186         has orientation PI/2, and a horizontal edge with dark side on top
    01187         has orientation 0. Obviously, the edge's orientation changes
    01188         by PI if the contrast is reversed.
    01189
    01190         */
    01191     float orientation;
    01192
    01193     Edgel()
    01194     : x(0.0f), y(0.0f), strength(0.0f), orientation(0.0f)
    01195     {}
    01196
    01197     Edgel(float ix, float iy, float is, float io)
    01198     : x(ix), y(iy), strength(is), orientation(io)
    01199     {}
    01200 };
    01201
    01202 template <class Image1, class Image2, class BackInsertable>
    01203 void internalCannyFindEdgels(Image1 const & gx,
    01204                              Image1 const & gy,
    01205                              Image2 const & magnitude,
    01206                              BackInsertable & edgels)
    01207 {
    01208     typedef typename Image1::value_type PixelType;
    01209     double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);
    01210
    01211     for(int y=1; y<gx.height()-1; ++y)
    01212     {
    01213         for(int x=1; x<gx.width()-1; ++x)
    01214         {
    01215             PixelType gradx = gx(x,y);
    01216             PixelType grady = gy(x,y);
    01217             double mag = magnitude(x, y);
    01218             if(mag == 0.0)
    01219                    continue;
    01220
    01221             int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
    01222             int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);
    01223
    01224             int x1 = x - dx,
    01225                 x2 = x + dx,
    01226                 y1 = y - dy,
    01227                 y2 = y + dy;
    01228
    01229             PixelType m1 = magnitude(x1, y1);
    01230             PixelType m3 = magnitude(x2, y2);
    01231
    01232             if(m1 < mag && m3 <= mag)
    01233             {
    01234                 Edgel edgel;
    01235
    01236                 // local maximum => quadratic interpolation of sub-pixel location
    01237                 PixelType del = (m1 - m3) / 2.0 / (m1 + m3 - 2.0*mag);
    01238                 edgel.x = x + dx*del;
    01239                 edgel.y = y + dy*del;
    01240                 edgel.strength = mag;
    01241                 double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
    01242                 if(orientation < 0.0)
    01243                     orientation += 2.0*M_PI;
    01244                 edgel.orientation = orientation;
    01245                 edgels.push_back(edgel);
    01246             }
    01247         }
    01248     }
    01249 }
    01250
    01251 /********************************************************/
    01252 /*                                                      */
    01253 /*                      cannyEdgelList                  */
    01254 /*                                                      */
    01255 /********************************************************/
    01256
    01257 /** \brief Simple implementation of Canny's edge detector.
    01258
    01259     This operator first calculates the gradient vector for each
    01260     pixel of the image using first derivatives of a Gaussian at the
    01261     given scale. Then a very simple non-maxima supression is performed:
    01262     for each 3x3 neighborhood, it is determined whether the center pixel has
    01263     larger gradient magnitude than its two neighbors in gradient direction
    01264     (where the direction is rounded into octands). If this is the case,
    01265     a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
    01266     edgel position is determined by fitting a parabola
    01267     to the three gradient magnitude values
    01268     mentioned above. The sub-pixel location of the parabola's tip
    01269     and the gradient magnitude and direction (from the pixel center)
    01270     are written in the newly created edgel.
    01271
    01272     <b> Declarations:</b>
    01273
    01274     pass arguments explicitly:
    01275     \code
    01276     namespace vigra {
    01277         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01278         void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
    01279                             BackInsertable & edgels, double scale);
    01280     }
    01281     \endcode
    01282
    01283     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01284     \code
    01285     namespace vigra {
    01286         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01287         void
    01288         cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01289                        BackInsertable & edgels, double scale);
    01290     }
    01291     \endcode
    01292
    01293     <b> Usage:</b>
    01294
    01295     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    01296     Namespace: vigra
    01297
    01298     \code
    01299     vigra::BImage src(w,h);
    01300
    01301     // empty edgel list
    01302     std::vector<vigra::Edgel> edgels;
    01303     ...
    01304
    01305     // find edgels at scale 0.8
    01306     vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
    01307     \endcode
    01308
    01309     <b> Required Interface:</b>
    01310
    01311     \code
    01312     SrcImageIterator src_upperleft;
    01313     SrcAccessor src_accessor;
    01314
    01315     src_accessor(src_upperleft);
    01316
    01317     BackInsertable edgels;
    01318     edgels.push_back(Edgel());
    01319     \endcode
    01320
    01321     SrcAccessor::value_type must be a type convertible to float
    01322
    01323     <b> Preconditions:</b>
    01324
    01325     \code
    01326     scale > 0
    01327     \endcode
    01328 */
    01329 doxygen_overloaded_function(template <...> void cannyEdgelList)
    01330
    01331 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01332 void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
    01333                         BackInsertable & edgels, double scale)
    01334 {
    01335     int w = lr.x - ul.x;
    01336     int h = lr.y - ul.y;
    01337
    01338     // calculate image gradients
    01339     typedef typename
    01340         NumericTraits<typename SrcAccessor::value_type>::RealPromote
    01341         TmpType;
    01342
    01343     BasicImage<TmpType> tmp(w,h), dx(w,h), dy(w,h);
    01344
    01345     Kernel1D<double> smooth, grad;
    01346
    01347     smooth.initGaussian(scale);
    01348     grad.initGaussianDerivative(scale, 1);
    01349
    01350     separableConvolveX(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
    01351     separableConvolveY(srcImageRange(tmp), destImage(dx), kernel1d(smooth));
    01352
    01353     separableConvolveY(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
    01354     separableConvolveX(srcImageRange(tmp), destImage(dy), kernel1d(smooth));
    01355
    01356     combineTwoImages(srcImageRange(dx), srcImage(dy), destImage(tmp),
    01357                      MagnitudeFunctor<TmpType>());
    01358
    01359
    01360     // find edgels
    01361     internalCannyFindEdgels(dx, dy, tmp, edgels);
    01362 }
    01363
    01364 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01365 inline void
    01366 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01367                BackInsertable & edgels, double scale)
    01368 {
    01369     cannyEdgelList(src.first, src.second, src.third, edgels, scale);
    01370 }
    01371
    01372 /********************************************************/
    01373 /*                                                      */
    01374 /*                       cannyEdgeImage                 */
    01375 /*                                                      */
    01376 /********************************************************/
    01377
    01378 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
    01379
    01380     This operator first calls \ref cannyEdgelList() to generate an
    01381     edgel list for the given image. Then it scans this list and selects edgels
    01382     whose strength is above the given <TT>gradient_threshold</TT>. For each of these
    01383     edgels, the edgel's location is rounded to the nearest pixel, and that
    01384     pixel marked with the given <TT>edge_marker</TT>.
    01385
    01386     <b> Declarations:</b>
    01387
    01388     pass arguments explicitly:
    01389     \code
    01390     namespace vigra {
    01391         template <class SrcIterator, class SrcAccessor,
    01392                   class DestIterator, class DestAccessor,
    01393                   class GradValue, class DestValue>
    01394         void cannyEdgeImage(
    01395                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01396                    DestIterator dul, DestAccessor da,
    01397                    double scale, GradValue gradient_threshold, DestValue edge_marker);
    01398     }
    01399     \endcode
    01400
    01401     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01402     \code
    01403     namespace vigra {
    01404         template <class SrcIterator, class SrcAccessor,
    01405                   class DestIterator, class DestAccessor,
    01406                   class GradValue, class DestValue>
    01407         void cannyEdgeImage(
    01408                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01409                    pair<DestIterator, DestAccessor> dest,
    01410                    double scale, GradValue gradient_threshold, DestValue edge_marker);
    01411     }
    01412     \endcode
    01413
    01414     <b> Usage:</b>
    01415
    01416     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    01417     Namespace: vigra
    01418
    01419     \code
    01420     vigra::BImage src(w,h), edges(w,h);
    01421
    01422     // empty edge image
    01423     edges = 0;
    01424     ...
    01425
    01426     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
    01427     vigra::cannyEdgeImage(srcImageRange(src), destImage(edges),
    01428                                      0.8, 4.0, 1);
    01429     \endcode
    01430
    01431     <b> Required Interface:</b>
    01432
    01433     see also: \ref cannyEdgelList().
    01434
    01435     \code
    01436     DestImageIterator dest_upperleft;
    01437     DestAccessor dest_accessor;
    01438     DestValue edge_marker;
    01439
    01440     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
    01441     \endcode
    01442
    01443     <b> Preconditions:</b>
    01444
    01445     \code
    01446     scale > 0
    01447     gradient_threshold > 0
    01448     \endcode
    01449 */
    01450 doxygen_overloaded_function(template <...> void cannyEdgeImage)
    01451
    01452 template <class SrcIterator, class SrcAccessor,
    01453           class DestIterator, class DestAccessor,
    01454           class GradValue, class DestValue>
    01455 void cannyEdgeImage(
    01456            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01457            DestIterator dul, DestAccessor da,
    01458            double scale, GradValue gradient_threshold, DestValue edge_marker)
    01459 {
    01460     std::vector<Edgel> edgels;
    01461
    01462     cannyEdgelList(sul, slr, sa, edgels, scale);
    01463
    01464     for(unsigned int i=0; i<edgels.size(); ++i)
    01465     {
    01466         if(gradient_threshold < edgels[i].strength)
    01467         {
    01468             Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
    01469
    01470             da.set(edge_marker, dul, pix);
    01471         }
    01472     }
    01473 }
    01474
    01475 template <class SrcIterator, class SrcAccessor,
    01476           class DestIterator, class DestAccessor,
    01477           class GradValue, class DestValue>
    01478 inline void cannyEdgeImage(
    01479            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01480            pair<DestIterator, DestAccessor> dest,
    01481            double scale, GradValue gradient_threshold, DestValue edge_marker)
    01482 {
    01483     cannyEdgeImage(src.first, src.second, src.third,
    01484                    dest.first, dest.second,
    01485                    scale, gradient_threshold, edge_marker);
    01486 }
    01487
    01488 /********************************************************/
    01489
    01490 namespace detail {
    01491
    01492 template <class DestIterator>
    01493 int neighborhoodConfiguration(DestIterator dul)
    01494 {
    01495     int v = 0;
    01496     NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
    01497     for(int i=0; i<8; ++i, --c)
    01498     {
    01499         v = (v << 1) | ((*c != 0) ? 1 : 0);
    01500     }
    01501
    01502     return v;
    01503 }
    01504
    01505 template <class GradValue>
    01506 struct SimplePoint
    01507 {
    01508     Diff2D point;
    01509     GradValue grad;
    01510
    01511     SimplePoint(Diff2D const & p, GradValue g)
    01512     : point(p), grad(g)
    01513     {}
    01514
    01515     bool operator<(SimplePoint const & o) const
    01516     {
    01517         return grad < o.grad;
    01518     }
    01519
    01520     bool operator>(SimplePoint const & o) const
    01521     {
    01522         return grad > o.grad;
    01523     }
    01524 };
    01525
    01526 template <class SrcIterator, class SrcAccessor,
    01527           class DestIterator, class DestAccessor,
    01528           class GradValue, class DestValue>
    01529 void cannyEdgeImageFromGrad(
    01530            SrcIterator sul, SrcIterator slr, SrcAccessor grad,
    01531            DestIterator dul, DestAccessor da,
    01532            GradValue gradient_threshold, DestValue edge_marker)
    01533 {
    01534     typedef typename SrcAccessor::value_type PixelType;
    01535     typedef typename NormTraits<PixelType>::SquaredNormType NormType;
    01536
    01537     NormType zero = NumericTraits<NormType>::zero();
    01538     double tan22_5 = M_SQRT2 - 1.0;
    01539     typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold);
    01540
    01541     int w = slr.x - sul.x;
    01542     int h = slr.y - sul.y;
    01543
    01544     sul += Diff2D(1,1);
    01545     dul += Diff2D(1,1);
    01546     Diff2D p(0,0);
    01547
    01548     for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y)
    01549     {
    01550         SrcIterator sx = sul;
    01551         DestIterator dx = dul;
    01552         for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x)
    01553         {
    01554             PixelType g = grad(sx);
    01555             NormType g2n = squaredNorm(g);
    01556             if(g2n < g2thresh)
    01557                 continue;
    01558
    01559             NormType g2n1, g2n3;
    01560             // find out quadrant
    01561             if(abs(g[1]) < tan22_5*abs(g[0]))
    01562             {
    01563                 // north-south edge
    01564                 g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0)));
    01565                 g2n3 = squaredNorm(grad(sx, Diff2D(1, 0)));
    01566             }
    01567             else if(abs(g[0]) < tan22_5*abs(g[1]))
    01568             {
    01569                 // west-east edge
    01570                 g2n1 = squaredNorm(grad(sx, Diff2D(0, -1)));
    01571                 g2n3 = squaredNorm(grad(sx, Diff2D(0, 1)));
    01572             }
    01573             else if(g[0]*g[1] < zero)
    01574             {
    01575                 // north-west-south-east edge
    01576                 g2n1 = squaredNorm(grad(sx, Diff2D(1, -1)));
    01577                 g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1)));
    01578             }
    01579             else
    01580             {
    01581                 // north-east-south-west edge
    01582                 g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1)));
    01583                 g2n3 = squaredNorm(grad(sx, Diff2D(1, 1)));
    01584             }
    01585
    01586             if(g2n1 < g2n && g2n3 <= g2n)
    01587             {
    01588                 da.set(edge_marker, dx);
    01589             }
    01590         }
    01591     }
    01592 }
    01593
    01594 } // namespace detail
    01595
    01596 /********************************************************/
    01597 /*                                                      */
    01598 /*              cannyEdgeImageWithThinning              */
    01599 /*                                                      */
    01600 /********************************************************/
    01601
    01602 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
    01603
    01604     The input pixels of this algorithms must be vectors of length 2 (see Required Interface below).
    01605     It first searches for all pixels whose gradient magnitude is larger
    01606     than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors
    01607     in gradient direction (where these neighbors are determined by nearest neighbor
    01608     interpolation, i.e. according to the octant where the gradient points into).
    01609     The resulting edge pixel candidates are then subjected to topological thinning
    01610     so that the remaining edge pixels can be linked into edgel chains with a provable,
    01611     non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient
    01612     magnitude survive. Optionally, the outermost pixels are marked as edge pixels
    01613     as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination
    01614     image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched).
    01615
    01616     <b> Declarations:</b>
    01617
    01618     pass arguments explicitly:
    01619     \code
    01620     namespace vigra {
    01621         template <class SrcIterator, class SrcAccessor,
    01622                   class DestIterator, class DestAccessor,
    01623                   class GradValue, class DestValue>
    01624         void cannyEdgeImageFromGradWithThinning(
    01625                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01626                    DestIterator dul, DestAccessor da,
    01627                    GradValue gradient_threshold,
    01628                    DestValue edge_marker, bool addBorder = true);
    01629     }
    01630     \endcode
    01631
    01632     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01633     \code
    01634     namespace vigra {
    01635         template <class SrcIterator, class SrcAccessor,
    01636                   class DestIterator, class DestAccessor,
    01637                   class GradValue, class DestValue>
    01638         void cannyEdgeImageFromGradWithThinning(
    01639                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01640                    pair<DestIterator, DestAccessor> dest,
    01641                    GradValue gradient_threshold,
    01642                    DestValue edge_marker, bool addBorder = true);
    01643     }
    01644     \endcode
    01645
    01646     <b> Usage:</b>
    01647
    01648     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    01649     Namespace: vigra
    01650
    01651     \code
    01652     vigra::BImage src(w,h), edges(w,h);
    01653
    01654     vigra::FVector2Image grad(w,h);
    01655     // compute the image gradient at scale 0.8
    01656     vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8);
    01657
    01658     // empty edge image
    01659     edges = 0;
    01660     // find edges gradient larger than 4.0, mark with 1, and add border
    01661     vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
    01662                                               4.0, 1, true);
    01663     \endcode
    01664
    01665     <b> Required Interface:</b>
    01666
    01667     \code
    01668     // the input pixel type must be a vector with two elements
    01669     SrcImageIterator src_upperleft;
    01670     SrcAccessor src_accessor;
    01671     typedef SrcAccessor::value_type SrcPixel;
    01672     typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType;
    01673
    01674     SrcPixel g = src_accessor(src_upperleft);
    01675     SrcPixel::value_type g0 = g[0];
    01676     SrcSquaredNormType gn = squaredNorm(g);
    01677
    01678     DestImageIterator dest_upperleft;
    01679     DestAccessor dest_accessor;
    01680     DestValue edge_marker;
    01681
    01682     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
    01683     \endcode
    01684
    01685     <b> Preconditions:</b>
    01686
    01687     \code
    01688     gradient_threshold > 0
    01689     \endcode
    01690 */
    01691 doxygen_overloaded_function(template <...> void cannyEdgeImageFromGradWithThinning)
    01692
    01693 template <class SrcIterator, class SrcAccessor,
    01694           class DestIterator, class DestAccessor,
    01695           class GradValue, class DestValue>
    01696 void cannyEdgeImageFromGradWithThinning(
    01697            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01698            DestIterator dul, DestAccessor da,
    01699            GradValue gradient_threshold,
    01700            DestValue edge_marker, bool addBorder)
    01701 {
    01702     int w = slr.x - sul.x;
    01703     int h = slr.y - sul.y;
    01704
    01705     BImage edgeImage(w, h, BImage::value_type(0));
    01706     BImage::traverser eul = edgeImage.upperLeft();
    01707     BImage::Accessor ea = edgeImage.accessor();
    01708     if(addBorder)
    01709         initImageBorder(destImageRange(edgeImage), 1, 1);
    01710     detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1);
    01711
    01712     static bool isSimplePoint[256] = {
    01713         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
    01714         0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
    01715         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
    01716         1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
    01717         0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
    01718         0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0,
    01719         0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,
    01720         1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
    01721         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
    01722         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    01723         0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
    01724         0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0,
    01725         1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
    01726         0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
    01727         1, 0, 1, 0 };
    01728
    01729     eul += Diff2D(1,1);
    01730     sul += Diff2D(1,1);
    01731     int w2 = w-2;
    01732     int h2 = h-2;
    01733
    01734     typedef detail::SimplePoint<GradValue> SP;
    01735     // use std::greater becaus we need the smallest gradients at the top of the queue
    01736     std::priority_queue<SP, std::vector<SP>, std::greater<SP> >  pqueue;
    01737
    01738     Diff2D p(0,0);
    01739     for(; p.y < h2; ++p.y)
    01740     {
    01741         for(p.x = 0; p.x < w2; ++p.x)
    01742         {
    01743             BImage::traverser e = eul + p;
    01744             if(*e == 0)
    01745                 continue;
    01746             int v = detail::neighborhoodConfiguration(e);
    01747             if(isSimplePoint[v])
    01748             {
    01749                 pqueue.push(SP(p, norm(sa(sul+p))));
    01750                 *e = 2; // remember that it is already in queue
    01751             }
    01752         }
    01753     }
    01754
    01755     static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
    01756                                    Diff2D(1,0),  Diff2D(0,1) };
    01757
    01758     while(pqueue.size())
    01759     {
    01760         p = pqueue.top().point;
    01761         pqueue.pop();
    01762
    01763         BImage::traverser e = eul + p;
    01764         int v = detail::neighborhoodConfiguration(e);
    01765         if(!isSimplePoint[v])
    01766             continue; // point may no longer be simple because its neighbors changed
    01767
    01768         *e = 0; // delete simple point
    01769
    01770         for(int i=0; i<4; ++i)
    01771         {
    01772             Diff2D pneu = p + dist[i];
    01773             if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
    01774                 continue; // do not remove points at the border
    01775
    01776             BImage::traverser eneu = eul + pneu;
    01777             if(*eneu == 1) // point is boundary and not yet in the queue
    01778             {
    01779                 int v = detail::neighborhoodConfiguration(eneu);
    01780                 if(isSimplePoint[v])
    01781                 {
    01782                     pqueue.push(SP(pneu, norm(sa(sul+pneu))));
    01783                     *eneu = 2; // remember that it is already in queue
    01784                 }
    01785             }
    01786         }
    01787     }
    01788
    01789     initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
    01790                 maskImage(edgeImage), edge_marker);
    01791 }
    01792
    01793 template <class SrcIterator, class SrcAccessor,
    01794           class DestIterator, class DestAccessor,
    01795           class GradValue, class DestValue>
    01796 inline void cannyEdgeImageFromGradWithThinning(
    01797            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01798            pair<DestIterator, DestAccessor> dest,
    01799            GradValue gradient_threshold,
    01800            DestValue edge_marker, bool addBorder)
    01801 {
    01802     cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
    01803                                dest.first, dest.second,
    01804                                gradient_threshold, edge_marker, addBorder);
    01805 }
    01806
    01807 template <class SrcIterator, class SrcAccessor,
    01808           class DestIterator, class DestAccessor,
    01809           class GradValue, class DestValue>
    01810 inline void cannyEdgeImageFromGradWithThinning(
    01811            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01812            DestIterator dul, DestAccessor da,
    01813            GradValue gradient_threshold, DestValue edge_marker)
    01814 {
    01815     cannyEdgeImageFromGradWithThinning(sul, slr, sa,
    01816                                dul, da,
    01817                                gradient_threshold, edge_marker, true);
    01818 }
    01819
    01820 template <class SrcIterator, class SrcAccessor,
    01821           class DestIterator, class DestAccessor,
    01822           class GradValue, class DestValue>
    01823 inline void cannyEdgeImageFromGradWithThinning(
    01824            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01825            pair<DestIterator, DestAccessor> dest,
    01826            GradValue gradient_threshold, DestValue edge_marker)
    01827 {
    01828     cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
    01829                                dest.first, dest.second,
    01830                                gradient_threshold, edge_marker, true);
    01831 }
    01832
    01833 /********************************************************/
    01834 /*                                                      */
    01835 /*              cannyEdgeImageWithThinning              */
    01836 /*                                                      */
    01837 /********************************************************/
    01838
    01839 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
    01840
    01841     This operator first calls \ref gaussianGradient() to compute the gradient of the input
    01842     image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an
    01843     edge image. See there for more detailed documentation.
    01844
    01845     <b> Declarations:</b>
    01846
    01847     pass arguments explicitly:
    01848     \code
    01849     namespace vigra {
    01850         template <class SrcIterator, class SrcAccessor,
    01851                   class DestIterator, class DestAccessor,
    01852                   class GradValue, class DestValue>
    01853         void cannyEdgeImageWithThinning(
    01854                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01855                    DestIterator dul, DestAccessor da,
    01856                    double scale, GradValue gradient_threshold,
    01857                    DestValue edge_marker, bool addBorder = true);
    01858     }
    01859     \endcode
    01860
    01861     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01862     \code
    01863     namespace vigra {
    01864         template <class SrcIterator, class SrcAccessor,
    01865                   class DestIterator, class DestAccessor,
    01866                   class GradValue, class DestValue>
    01867         void cannyEdgeImageWithThinning(
    01868                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01869                    pair<DestIterator, DestAccessor> dest,
    01870                    double scale, GradValue gradient_threshold,
    01871                    DestValue edge_marker, bool addBorder = true);
    01872     }
    01873     \endcode
    01874
    01875     <b> Usage:</b>
    01876
    01877     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    01878     Namespace: vigra
    01879
    01880     \code
    01881     vigra::BImage src(w,h), edges(w,h);
    01882
    01883     // empty edge image
    01884     edges = 0;
    01885     ...
    01886
    01887     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
    01888     vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges),
    01889                                      0.8, 4.0, 1, true);
    01890     \endcode
    01891
    01892     <b> Required Interface:</b>
    01893
    01894     see also: \ref cannyEdgelList().
    01895
    01896     \code
    01897     DestImageIterator dest_upperleft;
    01898     DestAccessor dest_accessor;
    01899     DestValue edge_marker;
    01900
    01901     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
    01902     \endcode
    01903
    01904     <b> Preconditions:</b>
    01905
    01906     \code
    01907     scale > 0
    01908     gradient_threshold > 0
    01909     \endcode
    01910 */
    01911 doxygen_overloaded_function(template <...> void cannyEdgeImageWithThinning)
    01912
    01913 template <class SrcIterator, class SrcAccessor,
    01914           class DestIterator, class DestAccessor,
    01915           class GradValue, class DestValue>
    01916 void cannyEdgeImageWithThinning(
    01917            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01918            DestIterator dul, DestAccessor da,
    01919            double scale, GradValue gradient_threshold,
    01920            DestValue edge_marker, bool addBorder)
    01921 {
    01922     // mark pixels that are higher than their neighbors in gradient direction
    01923     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    01924     BasicImage<TinyVector<TmpType, 2> > grad(slr-sul);
    01925     gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale);
    01926     cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da),
    01927                                gradient_threshold, edge_marker, addBorder);
    01928 }
    01929
    01930 template <class SrcIterator, class SrcAccessor,
    01931           class DestIterator, class DestAccessor,
    01932           class GradValue, class DestValue>
    01933 inline void cannyEdgeImageWithThinning(
    01934            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01935            pair<DestIterator, DestAccessor> dest,
    01936            double scale, GradValue gradient_threshold,
    01937            DestValue edge_marker, bool addBorder)
    01938 {
    01939     cannyEdgeImageWithThinning(src.first, src.second, src.third,
    01940                                dest.first, dest.second,
    01941                                scale, gradient_threshold, edge_marker, addBorder);
    01942 }
    01943
    01944 template <class SrcIterator, class SrcAccessor,
    01945           class DestIterator, class DestAccessor,
    01946           class GradValue, class DestValue>
    01947 inline void cannyEdgeImageWithThinning(
    01948            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
    01949            DestIterator dul, DestAccessor da,
    01950            double scale, GradValue gradient_threshold, DestValue edge_marker)
    01951 {
    01952     cannyEdgeImageWithThinning(sul, slr, sa,
    01953                                dul, da,
    01954                                scale, gradient_threshold, edge_marker, true);
    01955 }
    01956
    01957 template <class SrcIterator, class SrcAccessor,
    01958           class DestIterator, class DestAccessor,
    01959           class GradValue, class DestValue>
    01960 inline void cannyEdgeImageWithThinning(
    01961            triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01962            pair<DestIterator, DestAccessor> dest,
    01963            double scale, GradValue gradient_threshold, DestValue edge_marker)
    01964 {
    01965     cannyEdgeImageWithThinning(src.first, src.second, src.third,
    01966                                dest.first, dest.second,
    01967                                scale, gradient_threshold, edge_marker, true);
    01968 }
    01969
    01970 /********************************************************/
    01971
    01972 template <class Image1, class Image2, class BackInsertable>
    01973 void internalCannyFindEdgels3x3(Image1 const & grad,
    01974                                 Image2 const & mask,
    01975                                 BackInsertable & edgels)
    01976 {
    01977     typedef typename Image1::value_type PixelType;
    01978     typedef typename PixelType::value_type ValueType;
    01979
    01980     for(int y=1; y<grad.height()-1; ++y)
    01981     {
    01982         for(int x=1; x<grad.width()-1; ++x)
    01983         {
    01984             if(!mask(x,y))
    01985                 continue;
    01986
    01987             ValueType gradx = grad(x,y)[0];
    01988             ValueType grady = grad(x,y)[1];
    01989             double mag = hypot(gradx, grady);
    01990             if(mag == 0.0)
    01991                    continue;
    01992             double c = gradx / mag,
    01993                    s = grady / mag;
    01994
    01995             Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1);
    01996             l(0,0) = 1.0;
    01997
    01998             for(int yy = -1; yy <= 1; ++yy)
    01999             {
    02000                 for(int xx = -1; xx <= 1; ++xx)
    02001                 {
    02002                     double u = c*xx + s*yy;
    02003                     double v = norm(grad(x+xx, y+yy));
    02004                     l(1,0) = u;
    02005                     l(2,0) = u*u;
    02006                     ml += outer(l);
    02007                     mr += v*l;
    02008                 }
    02009             }
    02010
    02011             linearSolve(ml, mr, r);
    02012
    02013             Edgel edgel;
    02014
    02015             // local maximum => quadratic interpolation of sub-pixel location
    02016             double del = -r(1,0) / 2.0 / r(2,0);
    02017             if(std::fabs(del) > 1.5)  // don't move by more than about a pixel diameter
    02018                 del = 0.0;
    02019             edgel.x = x + c*del;
    02020             edgel.y = y + s*del;
    02021             edgel.strength = mag;
    02022             double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
    02023             if(orientation < 0.0)
    02024                 orientation += 2.0*M_PI;
    02025             edgel.orientation = orientation;
    02026             edgels.push_back(edgel);
    02027         }
    02028     }
    02029 }
    02030
    02031
    02032 /********************************************************/
    02033 /*                                                      */
    02034 /*                   cannyEdgelList3x3                  */
    02035 /*                                                      */
    02036 /********************************************************/
    02037
    02038 /** \brief Improved implementation of Canny's edge detector.
    02039
    02040     This operator first computes pixels which are crossed by the edge using
    02041     cannyEdgeImageWithThinning(). The gradient magnitude in the 3x3 neighborhood of these
    02042     pixels are then projected onto the normal of the edge (as determined
    02043     by the gradient direction). The edgel's subpixel location is found by fitting a
    02044     parabola through the 9 gradient values and determining the parabola's tip.
    02045     A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola
    02046     is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher.
    02047
    02048     <b> Declarations:</b>
    02049
    02050     pass arguments explicitly:
    02051     \code
    02052     namespace vigra {
    02053         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    02054         void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
    02055                                BackInsertable & edgels, double scale);
    02056     }
    02057     \endcode
    02058
    02059     use argument objects in conjunction with \ref ArgumentObjectFactories :
    02060     \code
    02061     namespace vigra {
    02062         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    02063         void
    02064         cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    02065                           BackInsertable & edgels, double scale);
    02066     }
    02067     \endcode
    02068
    02069     <b> Usage:</b>
    02070
    02071     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
    02072     Namespace: vigra
    02073
    02074     \code
    02075     vigra::BImage src(w,h);
    02076
    02077     // empty edgel list
    02078     std::vector<vigra::Edgel> edgels;
    02079     ...
    02080
    02081     // find edgels at scale 0.8
    02082     vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
    02083     \endcode
    02084
    02085     <b> Required Interface:</b>
    02086
    02087     \code
    02088     SrcImageIterator src_upperleft;
    02089     SrcAccessor src_accessor;
    02090
    02091     src_accessor(src_upperleft);
    02092
    02093     BackInsertable edgels;
    02094     edgels.push_back(Edgel());
    02095     \endcode
    02096
    02097     SrcAccessor::value_type must be a type convertible to float
    02098
    02099     <b> Preconditions:</b>
    02100
    02101     \code
    02102     scale > 0
    02103     \endcode
    02104 */
    02105 doxygen_overloaded_function(template <...> void cannyEdgelList3x3)
    02106
    02107 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    02108 void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
    02109                         BackInsertable & edgels, double scale)
    02110 {
    02111     int w = lr.x - ul.x;
    02112     int h = lr.y - ul.y;
    02113
    02114     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    02115     BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
    02116     gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
    02117
    02118     UInt8Image edges(lr-ul);
    02119     cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
    02120                                        0.0, 1, false);
    02121
    02122     // find edgels
    02123     internalCannyFindEdgels3x3(grad, edges, edgels);
    02124 }
    02125
    02126 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    02127 inline void
    02128 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    02129                BackInsertable & edgels, double scale)
    02130 {
    02131     cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale);
    02132 }
    02133
    02134
    02135
    02136 //@}
    02137
    02138 /** \page CrackEdgeImage Crack Edge Image
    02139
    02140 Crack edges are marked <i>between</i> the pixels of an image.
    02141 A Crack Edge Image is an image that represents these edges. In order
    02142 to accomodate the cracks, the Crack Edge Image must be twice as large
    02143 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
    02144 can easily be derived from a binary image or from the signs of the
    02145 response of a Laplacean filter. Consider the following sketch, where
    02146 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
    02147 <TT>*</TT> the resulting crack edges.
    02148
    02149     \code
    02150 sign of difference image         insert cracks         resulting CrackEdgeImage
    02151
    02152                                    + . - . -              . * . . .
    02153       + - -                        . . . . .              . * * * .
    02154       + + -               =>       + . + . -      =>      . . . * .
    02155       + + +                        . . . . .              . . . * *
    02156                                    + . + . +              . . . . .
    02157     \endcode
    02158
    02159 Starting from the original binary image (left), we insert crack pixels
    02160 to get to the double-sized image (center). Finally, we mark all
    02161 crack pixels whose non-crack neighbors have different signs as
    02162 crack edge points, while all other pixels (crack and non-crack) become
    02163 region pixels.
    02164
    02165 <b>Requirements on a Crack Edge Image:</b>
    02166
    02167 <ul>
    02168     <li>Crack Edge Images have odd width and height.
    02169     <li>Crack pixels have at least one odd coordinate.
    02170     <li>Only crack pixels may be marked as edge points.
    02171     <li>Crack pixels with two odd coordinates must be marked as edge points
    02172         whenever any of their neighboring crack pixels was marked.
    02173 </ul>
    02174
    02175 The last two requirements ensure that both edges and regions are 4-connected.
    02176 Thus, 4-connectivity and 8-connectivity yield identical connected
    02177 components in a Crack Edge Image (so called <i>well-composedness</i>).
    02178 This ensures that Crack Edge Images have nice topological properties
    02179 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000).
    02180 */
    02181
    02182
    02183 } // namespace vigra
    02184
    02185 #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.6.0 (13 Aug 2008)