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

vigra/edgedetection.hxx VIGRA

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

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

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)