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

vigra/edgedetection.hxx

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)