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

vigra/edgedetection.hxx |

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