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

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

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

html generated using doxygen and Python
vigra 1.7.1 (Mon Apr 16 2012)