[ 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) / 2;
00893     int h = (slr.y - sul.y) / 2;
00894     int x, y;
00895 
00896     int count1, count2, count3;
00897 
00898     static const Diff2D right(1,0);
00899     static const Diff2D bottom(0,1);
00900     static const Diff2D left(-1,0);
00901     static const Diff2D top(0,-1);
00902 
00903     static const Diff2D leftdist[] = {
00904         Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
00905     static const Diff2D rightdist[] = {
00906         Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
00907     static const Diff2D topdist[] = {
00908         Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
00909     static const Diff2D bottomdist[] = {
00910         Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};
00911 
00912     int i;
00913 
00914     SrcIterator sy = sul + Diff2D(0,1);
00915     SrcIterator sx;
00916 
00917     // close 1-pixel wide gaps (x-direction)
00918     for(y=0; y<h; ++y, sy.y+=2)
00919     {
00920         sx = sy + Diff2D(2,0);
00921 
00922         for(x=2; x<w; ++x, sx.x+=2)
00923         {
00924             if(sa(sx) == edge_marker) continue;
00925 
00926             if(sa(sx, left) != edge_marker) continue;
00927             if(sa(sx, right) != edge_marker) continue;
00928 
00929             count1 = 0;
00930             count2 = 0;
00931             count3 = 0;
00932 
00933             for(i=0; i<4; ++i)
00934             {
00935                 if(sa(sx, leftdist[i]) == edge_marker)
00936                 {
00937                     ++count1;
00938                     count3 ^= 1 << i;
00939                 }
00940                 if(sa(sx, rightdist[i]) == edge_marker)
00941                 {
00942                     ++count2;
00943                     count3 ^= 1 << i;
00944                 }
00945             }
00946 
00947             if(count1 <= 1 || count2 <= 1 || count3 == 15)
00948             {
00949                 sa.set(edge_marker, sx);
00950             }
00951         }
00952    }
00953 
00954     sy = sul + Diff2D(1,2);
00955 
00956     // close 1-pixel wide gaps (y-direction)
00957     for(y=2; y<h; ++y, sy.y+=2)
00958     {
00959         sx = sy;
00960 
00961         for(x=0; x<w; ++x, sx.x+=2)
00962         {
00963             if(sa(sx) == edge_marker) continue;
00964 
00965             if(sa(sx, top) != edge_marker) continue;
00966             if(sa(sx, bottom) != edge_marker) continue;
00967 
00968             count1 = 0;
00969             count2 = 0;
00970             count3 = 0;
00971 
00972             for(i=0; i<4; ++i)
00973             {
00974                 if(sa(sx, topdist[i]) == edge_marker)
00975                 {
00976                     ++count1;
00977                     count3 ^= 1 << i;
00978                 }
00979                 if(sa(sx, bottomdist[i]) == edge_marker)
00980                 {
00981                     ++count2;
00982                     count3 ^= 1 << i;
00983                 }
00984             }
00985 
00986             if(count1 <= 1 || count2 <= 1 || count3 == 15)
00987             {
00988                 sa.set(edge_marker, sx);
00989             }
00990         }
00991     }
00992 }
00993 
00994 template <class SrcIterator, class SrcAccessor, class SrcValue>
00995 inline
00996 void closeGapsInCrackEdgeImage(
00997            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00998        SrcValue edge_marker)
00999 {
01000     closeGapsInCrackEdgeImage(src.first, src.second, src.third,
01001                     edge_marker);
01002 }
01003 
01004 /********************************************************/
01005 /*                                                      */
01006 /*              beautifyCrackEdgeImage                  */
01007 /*                                                      */
01008 /********************************************************/
01009 
01010 /** \brief Beautify crack edge image for visualization.
01011 
01012     This algorithm is applied as a post-processing operation of
01013     \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
01014     the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after
01015     application of this algorithm. In particular, the algorithm removes zero-cells
01016     marked as edges to avoid staircase effects on diagonal lines like this:
01017 
01018     \code
01019     original edge points (*)     resulting edge points
01020 
01021           . * . . .                   . * . . .
01022           . * * * .                   . . * . .
01023           . . . * .           =>      . . . * .
01024           . . . * *                   . . . . *
01025           . . . . .                   . . . . .
01026     \endcode
01027 
01028     Therfore, this algorithm should only be applied as a vizualization aid, i.e.
01029     for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>,
01030     and background pixels with <TT>background_marker</TT>. The image's value type must be
01031     equality comparable.
01032 
01033     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place,
01034     i.e. on only one image.
01035 
01036     <b> Declarations:</b>
01037 
01038     pass arguments explicitly:
01039     \code
01040     namespace vigra {
01041         template <class SrcIterator, class SrcAccessor, class SrcValue>
01042         void beautifyCrackEdgeImage(
01043                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01044                SrcValue edge_marker, SrcValue background_marker)
01045     }
01046     \endcode
01047 
01048     use argument objects in conjunction with \ref ArgumentObjectFactories :
01049     \code
01050     namespace vigra {
01051         template <class SrcIterator, class SrcAccessor, class SrcValue>
01052         void beautifyCrackEdgeImage(
01053                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01054                SrcValue edge_marker, SrcValue background_marker)
01055     }
01056     \endcode
01057 
01058     <b> Usage:</b>
01059 
01060         <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
01061     Namespace: vigra
01062 
01063     \code
01064     vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
01065 
01066     // empty edge image
01067     edges = 0;
01068     ...
01069 
01070     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
01071     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges),
01072                                          0.8, 4.0, 1);
01073 
01074     // beautify edge image for visualization
01075     vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
01076 
01077     // show to the user
01078     window.open(edges);
01079     \endcode
01080 
01081     <b> Required Interface:</b>
01082 
01083     \code
01084     SrcImageIterator src_upperleft, src_lowerright;
01085 
01086     SrcAccessor src_accessor;
01087     DestAccessor dest_accessor;
01088 
01089     SrcAccessor::value_type u = src_accessor(src_upperleft);
01090 
01091     u == u
01092     u != u
01093 
01094     SrcValue background_marker;
01095     src_accessor.set(background_marker, src_upperleft);
01096     \endcode
01097 */
01098 doxygen_overloaded_function(template <...> void beautifyCrackEdgeImage)
01099 
01100 template <class SrcIterator, class SrcAccessor, class SrcValue>
01101 void beautifyCrackEdgeImage(
01102                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01103            SrcValue edge_marker, SrcValue background_marker)
01104 {
01105     int w = (slr.x - sul.x) / 2;
01106     int h = (slr.y - sul.y) / 2;
01107     int x, y;
01108 
01109     SrcIterator sy = sul + Diff2D(1,1);
01110     SrcIterator sx;
01111 
01112     static const Diff2D right(1,0);
01113     static const Diff2D bottom(0,1);
01114     static const Diff2D left(-1,0);
01115     static const Diff2D top(0,-1);
01116 
01117     //  delete 0-cells at corners
01118     for(y=0; y<h; ++y, sy.y+=2)
01119     {
01120         sx = sy;
01121 
01122         for(x=0; x<w; ++x, sx.x+=2)
01123         {
01124             if(sa(sx) != edge_marker) continue;
01125 
01126             if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
01127             if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;
01128 
01129             sa.set(background_marker, sx);
01130         }
01131     }
01132 }
01133 
01134 template <class SrcIterator, class SrcAccessor, class SrcValue>
01135 inline
01136 void beautifyCrackEdgeImage(
01137            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01138            SrcValue edge_marker, SrcValue background_marker)
01139 {
01140     beautifyCrackEdgeImage(src.first, src.second, src.third,
01141                     edge_marker, background_marker);
01142 }
01143 
01144 
01145 /** Helper class that stores edgel attributes.
01146 */
01147 class Edgel
01148 {
01149   public:
01150   
01151         /** The type of an Edgel's members.
01152         */
01153     typedef float value_type;
01154 
01155         /** The edgel's sub-pixel x coordinate.
01156         */
01157     value_type x;
01158 
01159         /** The edgel's sub-pixel y coordinate.
01160         */
01161     value_type y;
01162 
01163         /** The edgel's strength (magnitude of the gradient vector).
01164         */
01165     value_type strength;
01166 
01167         /**
01168         The edgel's orientation. This is the clockwise angle in radians
01169         between the x-axis and the edge, so that the bright side of the
01170         edge is on the left when one looks along the orientation vector. 
01171         The angle is measured clockwise because the y-axis increases 
01172         downwards (left-handed coordinate system):
01173 
01174         \code
01175 
01176   edgel axis
01177        \  
01178   (dark \  (bright side)
01179   side)  \ 
01180           \ 
01181            +------------> x-axis
01182            |\    |
01183            | \ /_/  orientation angle
01184            |  \\
01185            |   \
01186            |   
01187     y-axis V
01188         \endcode
01189 
01190         So, for example a vertical edge with its dark side on the left
01191         has orientation PI/2, and a horizontal edge with dark side on top
01192         has orientation PI. Obviously, the edge's orientation changes
01193         by PI if the contrast is reversed.
01194         
01195         Note that this convention changed as of VIGRA version 1.7.0.
01196 
01197         */
01198     value_type orientation;
01199 
01200     Edgel()
01201     : x(0.0), y(0.0), strength(0.0), orientation(0.0)
01202     {}
01203 
01204     Edgel(value_type ix, value_type iy, value_type is, value_type io)
01205     : x(ix), y(iy), strength(is), orientation(io)
01206     {}
01207 };
01208 
01209 template <class SrcIterator, class SrcAccessor, 
01210           class MagnitudeImage, class BackInsertable>
01211 void internalCannyFindEdgels(SrcIterator ul, SrcAccessor grad,
01212                              MagnitudeImage const & magnitude,
01213                              BackInsertable & edgels)
01214 {
01215     typedef typename SrcAccessor::value_type PixelType;
01216     typedef typename PixelType::value_type ValueType;
01217 
01218     double t = 0.5 / VIGRA_CSTD::sin(M_PI/8.0);
01219 
01220     ul += Diff2D(1,1);
01221     for(int y=1; y<magnitude.height()-1; ++y, ++ul.y)
01222     {
01223         SrcIterator ix = ul;
01224         for(int x=1; x<magnitude.width()-1; ++x, ++ix.x)
01225         {
01226             double mag = magnitude(x, y);
01227             if(mag == 0.0)
01228                    continue;
01229             ValueType gradx = grad.getComponent(ix, 0);
01230             ValueType grady = grad.getComponent(ix, 1);
01231 
01232             int dx = (int)VIGRA_CSTD::floor(gradx*t/mag + 0.5);
01233             int dy = (int)VIGRA_CSTD::floor(grady*t/mag + 0.5);
01234 
01235             int x1 = x - dx,
01236                 x2 = x + dx,
01237                 y1 = y - dy,
01238                 y2 = y + dy;
01239 
01240             double m1 = magnitude(x1, y1);
01241             double m3 = magnitude(x2, y2);
01242 
01243             if(m1 < mag && m3 <= mag)
01244             {
01245                 Edgel edgel;
01246 
01247                 // local maximum => quadratic interpolation of sub-pixel location
01248                 double del = 0.5 * (m1 - m3) / (m1 + m3 - 2.0*mag);
01249                 edgel.x = Edgel::value_type(x + dx*del);
01250                 edgel.y = Edgel::value_type(y + dy*del);
01251                 edgel.strength = Edgel::value_type(mag);
01252                 double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
01253                 if(orientation < 0.0)
01254                     orientation += 2.0*M_PI;
01255                 edgel.orientation = Edgel::value_type(orientation);
01256                 edgels.push_back(edgel);
01257             }
01258         }
01259     }
01260 }
01261 
01262 /********************************************************/
01263 /*                                                      */
01264 /*                      cannyEdgelList                  */
01265 /*                                                      */
01266 /********************************************************/
01267 
01268 /** \brief Simple implementation of Canny's edge detector.
01269     
01270     The function can be called in two modes: If you pass a 'scale', it is assumed that the 
01271     original image is scalar, and the Gaussian gradient is internally computed at the
01272     given 'scale'. If the function is called without scale parameter, it is assumed that
01273     the given image already contains the gradient (i.e. its value_type must be 
01274     a vector of length 2).
01275 
01276     On the basis of the gradient image, a simple non-maxima supression is performed:
01277     for each 3x3 neighborhood, it is determined whether the center pixel has
01278     larger gradient magnitude than its two neighbors in gradient direction
01279     (where the direction is rounded into octands). If this is the case,
01280     a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
01281     edgel position is determined by fitting a parabola to the three gradient 
01282     magnitude values mentioned above. The sub-pixel location of the parabola's tip
01283     and the gradient magnitude and direction (from the pixel center)
01284     are written in the newly created edgel.
01285 
01286     <b> Declarations:</b>
01287 
01288     pass arguments explicitly:
01289     \code
01290     namespace vigra {
01291         // compute edgels from a gradient image
01292         template <class SrcIterator, class SrcAccessor, class BackInsertable>
01293         void 
01294         cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01295                        BackInsertable & edgels);
01296 
01297         // compute edgels from a scaler image (determine gradient internally at 'scale')
01298         template <class SrcIterator, class SrcAccessor, class BackInsertable>
01299         void
01300         cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01301                        BackInsertable & edgels, double scale);
01302     }
01303     \endcode
01304 
01305     use argument objects in conjunction with \ref ArgumentObjectFactories :
01306     \code
01307     namespace vigra {
01308         // compute edgels from a gradient image
01309         template <class SrcIterator, class SrcAccessor, class BackInsertable>
01310         void
01311         cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01312                        BackInsertable & edgels);
01313 
01314         // compute edgels from a scaler image (determine gradient internally at 'scale')
01315         template <class SrcIterator, class SrcAccessor, class BackInsertable>
01316         void
01317         cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01318                        BackInsertable & edgels, double scale);
01319     }
01320     \endcode
01321 
01322     <b> Usage:</b>
01323 
01324     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
01325     Namespace: vigra
01326 
01327     \code
01328     vigra::BImage src(w,h);
01329 
01330     // empty edgel list
01331     std::vector<vigra::Edgel> edgels;
01332     ...
01333 
01334     // find edgels at scale 0.8
01335     vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
01336     \endcode
01337 
01338     <b> Required Interface:</b>
01339 
01340     \code
01341     SrcImageIterator src_upperleft;
01342     SrcAccessor src_accessor;
01343 
01344     src_accessor(src_upperleft);
01345 
01346     BackInsertable edgels;
01347     edgels.push_back(Edgel());
01348     \endcode
01349 
01350     SrcAccessor::value_type must be a type convertible to float
01351 
01352     <b> Preconditions:</b>
01353 
01354     \code
01355     scale > 0
01356     \endcode
01357 */
01358 doxygen_overloaded_function(template <...> void cannyEdgelList)
01359 
01360 template <class SrcIterator, class SrcAccessor, class BackInsertable>
01361 void 
01362 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01363                BackInsertable & edgels, double scale)
01364 {
01365     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
01366     BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
01367     gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
01368 
01369     cannyEdgelList(srcImageRange(grad), edgels);
01370 }
01371 
01372 template <class SrcIterator, class SrcAccessor, class BackInsertable>
01373 inline void
01374 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01375                BackInsertable & edgels, double scale)
01376 {
01377     cannyEdgelList(src.first, src.second, src.third, edgels, scale);
01378 }
01379 
01380 template <class SrcIterator, class SrcAccessor, class BackInsertable>
01381 void 
01382 cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01383                BackInsertable & edgels)
01384 {
01385     using namespace functor;
01386     
01387     typedef typename SrcAccessor::value_type SrcType;
01388     typedef typename NumericTraits<typename SrcType::value_type>::RealPromote TmpType;
01389     BasicImage<TmpType> magnitude(lr-ul);
01390     transformImage(srcIterRange(ul, lr, src), destImage(magnitude), norm(Arg1()));
01391 
01392     // find edgels
01393     internalCannyFindEdgels(ul, src, magnitude, edgels);
01394 }
01395 
01396 template <class SrcIterator, class SrcAccessor, class BackInsertable>
01397 inline void
01398 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01399                BackInsertable & edgels)
01400 {
01401     cannyEdgelList(src.first, src.second, src.third, edgels);
01402 }
01403 
01404 
01405 /********************************************************/
01406 /*                                                      */
01407 /*                       cannyEdgeImage                 */
01408 /*                                                      */
01409 /********************************************************/
01410 
01411 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
01412 
01413     This operator first calls \ref cannyEdgelList() to generate an
01414     edgel list for the given image. Then it scans this list and selects edgels
01415     whose strength is above the given <TT>gradient_threshold</TT>. For each of these
01416     edgels, the edgel's location is rounded to the nearest pixel, and that
01417     pixel marked with the given <TT>edge_marker</TT>.
01418 
01419     <b> Declarations:</b>
01420 
01421     pass arguments explicitly:
01422     \code
01423     namespace vigra {
01424         template <class SrcIterator, class SrcAccessor,
01425                   class DestIterator, class DestAccessor,
01426                   class GradValue, class DestValue>
01427         void cannyEdgeImage(
01428                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01429                    DestIterator dul, DestAccessor da,
01430                    double scale, GradValue gradient_threshold, DestValue edge_marker);
01431     }
01432     \endcode
01433 
01434     use argument objects in conjunction with \ref ArgumentObjectFactories :
01435     \code
01436     namespace vigra {
01437         template <class SrcIterator, class SrcAccessor,
01438                   class DestIterator, class DestAccessor,
01439                   class GradValue, class DestValue>
01440         void cannyEdgeImage(
01441                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01442                    pair<DestIterator, DestAccessor> dest,
01443                    double scale, GradValue gradient_threshold, DestValue edge_marker);
01444     }
01445     \endcode
01446 
01447     <b> Usage:</b>
01448 
01449     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
01450     Namespace: vigra
01451 
01452     \code
01453     vigra::BImage src(w,h), edges(w,h);
01454 
01455     // empty edge image
01456     edges = 0;
01457     ...
01458 
01459     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1
01460     vigra::cannyEdgeImage(srcImageRange(src), destImage(edges),
01461                                      0.8, 4.0, 1);
01462     \endcode
01463 
01464     <b> Required Interface:</b>
01465 
01466     see also: \ref cannyEdgelList().
01467 
01468     \code
01469     DestImageIterator dest_upperleft;
01470     DestAccessor dest_accessor;
01471     DestValue edge_marker;
01472 
01473     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
01474     \endcode
01475 
01476     <b> Preconditions:</b>
01477 
01478     \code
01479     scale > 0
01480     gradient_threshold > 0
01481     \endcode
01482 */
01483 doxygen_overloaded_function(template <...> void cannyEdgeImage)
01484 
01485 template <class SrcIterator, class SrcAccessor,
01486           class DestIterator, class DestAccessor,
01487           class GradValue, class DestValue>
01488 void cannyEdgeImage(
01489            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01490            DestIterator dul, DestAccessor da,
01491            double scale, GradValue gradient_threshold, DestValue edge_marker)
01492 {
01493     std::vector<Edgel> edgels;
01494 
01495     cannyEdgelList(sul, slr, sa, edgels, scale);
01496     
01497     int w = slr.x - sul.x;
01498     int h = slr.y - sul.y;
01499 
01500     for(unsigned int i=0; i<edgels.size(); ++i)
01501     {
01502         if(gradient_threshold < edgels[i].strength)
01503         {
01504             Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
01505             
01506             if(pix.x < 0 || pix.x >= w || pix.y < 0 || pix.y >= h)
01507                 continue;
01508 
01509             da.set(edge_marker, dul, pix);
01510         }
01511     }
01512 }
01513 
01514 template <class SrcIterator, class SrcAccessor,
01515           class DestIterator, class DestAccessor,
01516           class GradValue, class DestValue>
01517 inline void cannyEdgeImage(
01518            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01519            pair<DestIterator, DestAccessor> dest,
01520            double scale, GradValue gradient_threshold, DestValue edge_marker)
01521 {
01522     cannyEdgeImage(src.first, src.second, src.third,
01523                    dest.first, dest.second,
01524                    scale, gradient_threshold, edge_marker);
01525 }
01526 
01527 /********************************************************/
01528 
01529 namespace detail {
01530 
01531 template <class DestIterator>
01532 int neighborhoodConfiguration(DestIterator dul)
01533 {
01534     int v = 0;
01535     NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
01536     for(int i=0; i<8; ++i, --c)
01537     {
01538         v = (v << 1) | ((*c != 0) ? 1 : 0);
01539     }
01540 
01541     return v;
01542 }
01543 
01544 template <class GradValue>
01545 struct SimplePoint
01546 {
01547     Diff2D point;
01548     GradValue grad;
01549 
01550     SimplePoint(Diff2D const & p, GradValue g)
01551     : point(p), grad(g)
01552     {}
01553 
01554     bool operator<(SimplePoint const & o) const
01555     {
01556         return grad < o.grad;
01557     }
01558 
01559     bool operator>(SimplePoint const & o) const
01560     {
01561         return grad > o.grad;
01562     }
01563 };
01564 
01565 template <class SrcIterator, class SrcAccessor,
01566           class DestIterator, class DestAccessor,
01567           class GradValue, class DestValue>
01568 void cannyEdgeImageFromGrad(
01569            SrcIterator sul, SrcIterator slr, SrcAccessor grad,
01570            DestIterator dul, DestAccessor da,
01571            GradValue gradient_threshold, DestValue edge_marker)
01572 {
01573     typedef typename SrcAccessor::value_type PixelType;
01574     typedef typename NormTraits<PixelType>::SquaredNormType NormType;
01575 
01576     NormType zero = NumericTraits<NormType>::zero();
01577     double tan22_5 = M_SQRT2 - 1.0;
01578     typename NormTraits<GradValue>::SquaredNormType g2thresh = squaredNorm(gradient_threshold);
01579 
01580     int w = slr.x - sul.x;
01581     int h = slr.y - sul.y;
01582 
01583     sul += Diff2D(1,1);
01584     dul += Diff2D(1,1);
01585     Diff2D p(0,0);
01586 
01587     for(int y = 1; y < h-1; ++y, ++sul.y, ++dul.y)
01588     {
01589         SrcIterator sx = sul;
01590         DestIterator dx = dul;
01591         for(int x = 1; x < w-1; ++x, ++sx.x, ++dx.x)
01592         {
01593             PixelType g = grad(sx);
01594             NormType g2n = squaredNorm(g);
01595             if(g2n < g2thresh)
01596                 continue;
01597 
01598             NormType g2n1, g2n3;
01599             // find out quadrant
01600             if(abs(g[1]) < tan22_5*abs(g[0]))
01601             {
01602                 // north-south edge
01603                 g2n1 = squaredNorm(grad(sx, Diff2D(-1, 0)));
01604                 g2n3 = squaredNorm(grad(sx, Diff2D(1, 0)));
01605             }
01606             else if(abs(g[0]) < tan22_5*abs(g[1]))
01607             {
01608                 // west-east edge
01609                 g2n1 = squaredNorm(grad(sx, Diff2D(0, -1)));
01610                 g2n3 = squaredNorm(grad(sx, Diff2D(0, 1)));
01611             }
01612             else if(g[0]*g[1] < zero)
01613             {
01614                 // north-west-south-east edge
01615                 g2n1 = squaredNorm(grad(sx, Diff2D(1, -1)));
01616                 g2n3 = squaredNorm(grad(sx, Diff2D(-1, 1)));
01617             }
01618             else
01619             {
01620                 // north-east-south-west edge
01621                 g2n1 = squaredNorm(grad(sx, Diff2D(-1, -1)));
01622                 g2n3 = squaredNorm(grad(sx, Diff2D(1, 1)));
01623             }
01624 
01625             if(g2n1 < g2n && g2n3 <= g2n)
01626             {
01627                 da.set(edge_marker, dx);
01628             }
01629         }
01630     }
01631 }
01632 
01633 } // namespace detail
01634 
01635 /********************************************************/
01636 /*                                                      */
01637 /*              cannyEdgeImageWithThinning              */
01638 /*                                                      */
01639 /********************************************************/
01640 
01641 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
01642 
01643     The input pixels of this algorithms must be vectors of length 2 (see Required Interface below).
01644     It first searches for all pixels whose gradient magnitude is larger
01645     than the given <tt>gradient_threshold</tt> and larger than the magnitude of its two neighbors
01646     in gradient direction (where these neighbors are determined by nearest neighbor
01647     interpolation, i.e. according to the octant where the gradient points into).
01648     The resulting edge pixel candidates are then subjected to topological thinning
01649     so that the remaining edge pixels can be linked into edgel chains with a provable,
01650     non-heuristic algorithm. Thinning is performed so that the pixels with highest gradient
01651     magnitude survive. Optionally, the outermost pixels are marked as edge pixels
01652     as well when <tt>addBorder</tt> is true. The remaining pixels will be marked in the destination
01653     image with the value of <tt>edge_marker</tt> (all non-edge pixels remain untouched).
01654 
01655     <b> Declarations:</b>
01656 
01657     pass arguments explicitly:
01658     \code
01659     namespace vigra {
01660         template <class SrcIterator, class SrcAccessor,
01661                   class DestIterator, class DestAccessor,
01662                   class GradValue, class DestValue>
01663         void cannyEdgeImageFromGradWithThinning(
01664                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01665                    DestIterator dul, DestAccessor da,
01666                    GradValue gradient_threshold,
01667                    DestValue edge_marker, bool addBorder = true);
01668     }
01669     \endcode
01670 
01671     use argument objects in conjunction with \ref ArgumentObjectFactories :
01672     \code
01673     namespace vigra {
01674         template <class SrcIterator, class SrcAccessor,
01675                   class DestIterator, class DestAccessor,
01676                   class GradValue, class DestValue>
01677         void cannyEdgeImageFromGradWithThinning(
01678                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01679                    pair<DestIterator, DestAccessor> dest,
01680                    GradValue gradient_threshold,
01681                    DestValue edge_marker, bool addBorder = true);
01682     }
01683     \endcode
01684 
01685     <b> Usage:</b>
01686 
01687     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
01688     Namespace: vigra
01689 
01690     \code
01691     vigra::BImage src(w,h), edges(w,h);
01692 
01693     vigra::FVector2Image grad(w,h);
01694     // compute the image gradient at scale 0.8
01695     vigra::gaussianGradient(srcImageRange(src), destImage(grad), 0.8);
01696 
01697     // empty edge image
01698     edges = 0;
01699     // find edges gradient larger than 4.0, mark with 1, and add border
01700     vigra::cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destImage(edges),
01701                                               4.0, 1, true);
01702     \endcode
01703 
01704     <b> Required Interface:</b>
01705 
01706     \code
01707     // the input pixel type must be a vector with two elements
01708     SrcImageIterator src_upperleft;
01709     SrcAccessor src_accessor;
01710     typedef SrcAccessor::value_type SrcPixel;
01711     typedef NormTraits<SrcPixel>::SquaredNormType SrcSquaredNormType;
01712 
01713     SrcPixel g = src_accessor(src_upperleft);
01714     SrcPixel::value_type g0 = g[0];
01715     SrcSquaredNormType gn = squaredNorm(g);
01716 
01717     DestImageIterator dest_upperleft;
01718     DestAccessor dest_accessor;
01719     DestValue edge_marker;
01720 
01721     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
01722     \endcode
01723 
01724     <b> Preconditions:</b>
01725 
01726     \code
01727     gradient_threshold > 0
01728     \endcode
01729 */
01730 doxygen_overloaded_function(template <...> void cannyEdgeImageFromGradWithThinning)
01731 
01732 template <class SrcIterator, class SrcAccessor,
01733           class DestIterator, class DestAccessor,
01734           class GradValue, class DestValue>
01735 void cannyEdgeImageFromGradWithThinning(
01736            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01737            DestIterator dul, DestAccessor da,
01738            GradValue gradient_threshold,
01739            DestValue edge_marker, bool addBorder)
01740 {
01741     int w = slr.x - sul.x;
01742     int h = slr.y - sul.y;
01743 
01744     BImage edgeImage(w, h, BImage::value_type(0));
01745     BImage::traverser eul = edgeImage.upperLeft();
01746     BImage::Accessor ea = edgeImage.accessor();
01747     if(addBorder)
01748         initImageBorder(destImageRange(edgeImage), 1, 1);
01749     detail::cannyEdgeImageFromGrad(sul, slr, sa, eul, ea, gradient_threshold, 1);
01750 
01751     static bool isSimplePoint[256] = {
01752         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0,
01753         0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0,
01754         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1,
01755         1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
01756         0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,
01757         0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0,
01758         0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0,
01759         1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,
01760         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
01761         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
01762         0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,
01763         0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0,
01764         1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0,
01765         0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1,
01766         1, 0, 1, 0 };
01767 
01768     eul += Diff2D(1,1);
01769     sul += Diff2D(1,1);
01770     int w2 = w-2;
01771     int h2 = h-2;
01772 
01773     typedef detail::SimplePoint<GradValue> SP;
01774     // use std::greater becaus we need the smallest gradients at the top of the queue
01775     std::priority_queue<SP, std::vector<SP>, std::greater<SP> >  pqueue;
01776 
01777     Diff2D p(0,0);
01778     for(; p.y < h2; ++p.y)
01779     {
01780         for(p.x = 0; p.x < w2; ++p.x)
01781         {
01782             BImage::traverser e = eul + p;
01783             if(*e == 0)
01784                 continue;
01785             int v = detail::neighborhoodConfiguration(e);
01786             if(isSimplePoint[v])
01787             {
01788                 pqueue.push(SP(p, norm(sa(sul+p))));
01789                 *e = 2; // remember that it is already in queue
01790             }
01791         }
01792     }
01793 
01794     static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
01795                                    Diff2D(1,0),  Diff2D(0,1) };
01796 
01797     while(pqueue.size())
01798     {
01799         p = pqueue.top().point;
01800         pqueue.pop();
01801 
01802         BImage::traverser e = eul + p;
01803         int v = detail::neighborhoodConfiguration(e);
01804         if(!isSimplePoint[v])
01805             continue; // point may no longer be simple because its neighbors changed
01806 
01807         *e = 0; // delete simple point
01808 
01809         for(int i=0; i<4; ++i)
01810         {
01811             Diff2D pneu = p + dist[i];
01812             if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
01813                 continue; // do not remove points at the border
01814 
01815             BImage::traverser eneu = eul + pneu;
01816             if(*eneu == 1) // point is boundary and not yet in the queue
01817             {
01818                 int v = detail::neighborhoodConfiguration(eneu);
01819                 if(isSimplePoint[v])
01820                 {
01821                     pqueue.push(SP(pneu, norm(sa(sul+pneu))));
01822                     *eneu = 2; // remember that it is already in queue
01823                 }
01824             }
01825         }
01826     }
01827 
01828     initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
01829                 maskImage(edgeImage), edge_marker);
01830 }
01831 
01832 template <class SrcIterator, class SrcAccessor,
01833           class DestIterator, class DestAccessor,
01834           class GradValue, class DestValue>
01835 inline void cannyEdgeImageFromGradWithThinning(
01836            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01837            pair<DestIterator, DestAccessor> dest,
01838            GradValue gradient_threshold,
01839            DestValue edge_marker, bool addBorder)
01840 {
01841     cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
01842                                dest.first, dest.second,
01843                                gradient_threshold, edge_marker, addBorder);
01844 }
01845 
01846 template <class SrcIterator, class SrcAccessor,
01847           class DestIterator, class DestAccessor,
01848           class GradValue, class DestValue>
01849 inline void cannyEdgeImageFromGradWithThinning(
01850            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01851            DestIterator dul, DestAccessor da,
01852            GradValue gradient_threshold, DestValue edge_marker)
01853 {
01854     cannyEdgeImageFromGradWithThinning(sul, slr, sa,
01855                                dul, da,
01856                                gradient_threshold, edge_marker, true);
01857 }
01858 
01859 template <class SrcIterator, class SrcAccessor,
01860           class DestIterator, class DestAccessor,
01861           class GradValue, class DestValue>
01862 inline void cannyEdgeImageFromGradWithThinning(
01863            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01864            pair<DestIterator, DestAccessor> dest,
01865            GradValue gradient_threshold, DestValue edge_marker)
01866 {
01867     cannyEdgeImageFromGradWithThinning(src.first, src.second, src.third,
01868                                dest.first, dest.second,
01869                                gradient_threshold, edge_marker, true);
01870 }
01871 
01872 /********************************************************/
01873 /*                                                      */
01874 /*              cannyEdgeImageWithThinning              */
01875 /*                                                      */
01876 /********************************************************/
01877 
01878 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
01879 
01880     This operator first calls \ref gaussianGradient() to compute the gradient of the input
01881     image, ad then \ref cannyEdgeImageFromGradWithThinning() to generate an
01882     edge image. See there for more detailed documentation.
01883 
01884     <b> Declarations:</b>
01885 
01886     pass arguments explicitly:
01887     \code
01888     namespace vigra {
01889         template <class SrcIterator, class SrcAccessor,
01890                   class DestIterator, class DestAccessor,
01891                   class GradValue, class DestValue>
01892         void cannyEdgeImageWithThinning(
01893                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01894                    DestIterator dul, DestAccessor da,
01895                    double scale, GradValue gradient_threshold,
01896                    DestValue edge_marker, bool addBorder = true);
01897     }
01898     \endcode
01899 
01900     use argument objects in conjunction with \ref ArgumentObjectFactories :
01901     \code
01902     namespace vigra {
01903         template <class SrcIterator, class SrcAccessor,
01904                   class DestIterator, class DestAccessor,
01905                   class GradValue, class DestValue>
01906         void cannyEdgeImageWithThinning(
01907                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01908                    pair<DestIterator, DestAccessor> dest,
01909                    double scale, GradValue gradient_threshold,
01910                    DestValue edge_marker, bool addBorder = true);
01911     }
01912     \endcode
01913 
01914     <b> Usage:</b>
01915 
01916     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
01917     Namespace: vigra
01918 
01919     \code
01920     vigra::BImage src(w,h), edges(w,h);
01921 
01922     // empty edge image
01923     edges = 0;
01924     ...
01925 
01926     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
01927     vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges),
01928                                      0.8, 4.0, 1, true);
01929     \endcode
01930 
01931     <b> Required Interface:</b>
01932 
01933     see also: \ref cannyEdgelList().
01934 
01935     \code
01936     DestImageIterator dest_upperleft;
01937     DestAccessor dest_accessor;
01938     DestValue edge_marker;
01939 
01940     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
01941     \endcode
01942 
01943     <b> Preconditions:</b>
01944 
01945     \code
01946     scale > 0
01947     gradient_threshold > 0
01948     \endcode
01949 */
01950 doxygen_overloaded_function(template <...> void cannyEdgeImageWithThinning)
01951 
01952 template <class SrcIterator, class SrcAccessor,
01953           class DestIterator, class DestAccessor,
01954           class GradValue, class DestValue>
01955 void cannyEdgeImageWithThinning(
01956            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01957            DestIterator dul, DestAccessor da,
01958            double scale, GradValue gradient_threshold,
01959            DestValue edge_marker, bool addBorder)
01960 {
01961     // mark pixels that are higher than their neighbors in gradient direction
01962     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
01963     BasicImage<TinyVector<TmpType, 2> > grad(slr-sul);
01964     gaussianGradient(srcIterRange(sul, slr, sa), destImage(grad), scale);
01965     cannyEdgeImageFromGradWithThinning(srcImageRange(grad), destIter(dul, da),
01966                                gradient_threshold, edge_marker, addBorder);
01967 }
01968 
01969 template <class SrcIterator, class SrcAccessor,
01970           class DestIterator, class DestAccessor,
01971           class GradValue, class DestValue>
01972 inline void cannyEdgeImageWithThinning(
01973            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01974            pair<DestIterator, DestAccessor> dest,
01975            double scale, GradValue gradient_threshold,
01976            DestValue edge_marker, bool addBorder)
01977 {
01978     cannyEdgeImageWithThinning(src.first, src.second, src.third,
01979                                dest.first, dest.second,
01980                                scale, gradient_threshold, edge_marker, addBorder);
01981 }
01982 
01983 template <class SrcIterator, class SrcAccessor,
01984           class DestIterator, class DestAccessor,
01985           class GradValue, class DestValue>
01986 inline void cannyEdgeImageWithThinning(
01987            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01988            DestIterator dul, DestAccessor da,
01989            double scale, GradValue gradient_threshold, DestValue edge_marker)
01990 {
01991     cannyEdgeImageWithThinning(sul, slr, sa,
01992                                dul, da,
01993                                scale, gradient_threshold, edge_marker, true);
01994 }
01995 
01996 template <class SrcIterator, class SrcAccessor,
01997           class DestIterator, class DestAccessor,
01998           class GradValue, class DestValue>
01999 inline void cannyEdgeImageWithThinning(
02000            triple<SrcIterator, SrcIterator, SrcAccessor> src,
02001            pair<DestIterator, DestAccessor> dest,
02002            double scale, GradValue gradient_threshold, DestValue edge_marker)
02003 {
02004     cannyEdgeImageWithThinning(src.first, src.second, src.third,
02005                                dest.first, dest.second,
02006                                scale, gradient_threshold, edge_marker, true);
02007 }
02008 
02009 /********************************************************/
02010 
02011 template <class SrcIterator, class SrcAccessor, 
02012           class MaskImage, class BackInsertable>
02013 void internalCannyFindEdgels3x3(SrcIterator ul, SrcAccessor grad,
02014                                 MaskImage const & mask,
02015                                 BackInsertable & edgels)
02016 {
02017     typedef typename SrcAccessor::value_type PixelType;
02018     typedef typename PixelType::value_type ValueType;
02019 
02020     ul += Diff2D(1,1);
02021     for(int y=1; y<mask.height()-1; ++y, ++ul.y)
02022     {
02023         SrcIterator ix = ul;
02024         for(int x=1; x<mask.width()-1; ++x, ++ix.x)
02025         {
02026             if(!mask(x,y))
02027                 continue;
02028 
02029             ValueType gradx = grad.getComponent(ix, 0);
02030             ValueType grady = grad.getComponent(ix, 1);
02031             double mag = hypot(gradx, grady);
02032             if(mag == 0.0)
02033                    continue;
02034             double c = gradx / mag,
02035                    s = grady / mag;
02036 
02037             Matrix<double> ml(3,3), mr(3,1), l(3,1), r(3,1);
02038             l(0,0) = 1.0;
02039 
02040             for(int yy = -1; yy <= 1; ++yy)
02041             {
02042                 for(int xx = -1; xx <= 1; ++xx)
02043                 {
02044                     double u = c*xx + s*yy;
02045                     double v = norm(grad(ix, Diff2D(xx, yy)));
02046                     l(1,0) = u;
02047                     l(2,0) = u*u;
02048                     ml += outer(l);
02049                     mr += v*l;
02050                 }
02051             }
02052 
02053             linearSolve(ml, mr, r);
02054 
02055             Edgel edgel;
02056 
02057             // local maximum => quadratic interpolation of sub-pixel location
02058             double del = -r(1,0) / 2.0 / r(2,0);
02059             if(std::fabs(del) > 1.5)  // don't move by more than about a pixel diameter
02060                 del = 0.0;
02061             edgel.x = Edgel::value_type(x + c*del);
02062             edgel.y = Edgel::value_type(y + s*del);
02063             edgel.strength = Edgel::value_type(mag);
02064             double orientation = VIGRA_CSTD::atan2(grady, gradx) + 0.5*M_PI;
02065             if(orientation < 0.0)
02066                 orientation += 2.0*M_PI;
02067             edgel.orientation = Edgel::value_type(orientation);
02068             edgels.push_back(edgel);
02069         }
02070     }
02071 }
02072 
02073 
02074 /********************************************************/
02075 /*                                                      */
02076 /*                   cannyEdgelList3x3                  */
02077 /*                                                      */
02078 /********************************************************/
02079 
02080 /** \brief Improved implementation of Canny's edge detector.
02081 
02082     This operator first computes pixels which are crossed by the edge using
02083     cannyEdgeImageWithThinning(). The gradient magnitudes in the 3x3 neighborhood of these
02084     pixels are then projected onto the normal of the edge (as determined
02085     by the gradient direction). The edgel's subpixel location is found by fitting a
02086     parabola through the 9 gradient values and determining the parabola's tip.
02087     A new \ref Edgel is appended to the given vector of <TT>edgels</TT>. Since the parabola
02088     is fitted to 9 points rather than 3 points as in cannyEdgelList(), the accuracy is higher.
02089     
02090     The function can be called in two modes: If you pass a 'scale', it is assumed that the 
02091     original image is scalar, and the Gaussian gradient is internally computed at the
02092     given 'scale'. If the function is called without scale parameter, it is assumed that
02093     the given image already contains the gradient (i.e. its value_type must be 
02094     a vector of length 2).
02095 
02096     <b> Declarations:</b>
02097 
02098     pass arguments explicitly:
02099     \code
02100     namespace vigra {
02101         // compute edgels from a gradient image
02102         template <class SrcIterator, class SrcAccessor, class BackInsertable>
02103         void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
02104                                BackInsertable & edgels);
02105 
02106         // compute edgels from a scaler image (determine gradient internally at 'scale')
02107         template <class SrcIterator, class SrcAccessor, class BackInsertable>
02108         void cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
02109                                BackInsertable & edgels, double scale);
02110     }
02111     \endcode
02112 
02113     use argument objects in conjunction with \ref ArgumentObjectFactories :
02114     \code
02115     namespace vigra {
02116         // compute edgels from a gradient image
02117         template <class SrcIterator, class SrcAccessor, class BackInsertable>
02118         void
02119         cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
02120                           BackInsertable & edgels);
02121 
02122         // compute edgels from a scaler image (determine gradient internally at 'scale')
02123         template <class SrcIterator, class SrcAccessor, class BackInsertable>
02124         void
02125         cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
02126                           BackInsertable & edgels, double scale);
02127     }
02128     \endcode
02129 
02130     <b> Usage:</b>
02131 
02132     <b>\#include</b> <<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>><br>
02133     Namespace: vigra
02134 
02135     \code
02136     vigra::BImage src(w,h);
02137 
02138     // empty edgel list
02139     std::vector<vigra::Edgel> edgels;
02140     ...
02141 
02142     // find edgels at scale 0.8
02143     vigra::cannyEdgelList3x3(srcImageRange(src), edgels, 0.8);
02144     \endcode
02145 
02146     <b> Required Interface:</b>
02147 
02148     \code
02149     SrcImageIterator src_upperleft;
02150     SrcAccessor src_accessor;
02151 
02152     src_accessor(src_upperleft);
02153 
02154     BackInsertable edgels;
02155     edgels.push_back(Edgel());
02156     \endcode
02157 
02158     SrcAccessor::value_type must be a type convertible to float
02159 
02160     <b> Preconditions:</b>
02161 
02162     \code
02163     scale > 0
02164     \endcode
02165 */
02166 doxygen_overloaded_function(template <...> void cannyEdgelList3x3)
02167 
02168 template <class SrcIterator, class SrcAccessor, class BackInsertable>
02169 void 
02170 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
02171                   BackInsertable & edgels, double scale)
02172 {
02173     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
02174     BasicImage<TinyVector<TmpType, 2> > grad(lr-ul);
02175     gaussianGradient(srcIterRange(ul, lr, src), destImage(grad), scale);
02176 
02177     cannyEdgelList3x3(srcImageRange(grad), edgels);
02178 }
02179 
02180 template <class SrcIterator, class SrcAccessor, class BackInsertable>
02181 inline void
02182 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
02183                   BackInsertable & edgels, double scale)
02184 {
02185     cannyEdgelList3x3(src.first, src.second, src.third, edgels, scale);
02186 }
02187 
02188 template <class SrcIterator, class SrcAccessor, class BackInsertable>
02189 void 
02190 cannyEdgelList3x3(SrcIterator ul, SrcIterator lr, SrcAccessor src,
02191                   BackInsertable & edgels)
02192 {
02193     UInt8Image edges(lr-ul);
02194     cannyEdgeImageFromGradWithThinning(srcIterRange(ul, lr, src), destImage(edges),
02195                                        0.0, 1, false);
02196 
02197     // find edgels
02198     internalCannyFindEdgels3x3(ul, src, edges, edgels);
02199 }
02200 
02201 template <class SrcIterator, class SrcAccessor, class BackInsertable>
02202 inline void
02203 cannyEdgelList3x3(triple<SrcIterator, SrcIterator, SrcAccessor> src,
02204                   BackInsertable & edgels)
02205 {
02206     cannyEdgelList3x3(src.first, src.second, src.third, edgels);
02207 }
02208 
02209 //@}
02210 
02211 /** \page CrackEdgeImage Crack Edge Image
02212 
02213 Crack edges are marked <i>between</i> the pixels of an image.
02214 A Crack Edge Image is an image that represents these edges. In order
02215 to accomodate the cracks, the Crack Edge Image must be twice as large
02216 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
02217 can easily be derived from a binary image or from the signs of the
02218 response of a Laplacean filter. Consider the following sketch, where
02219 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
02220 <TT>*</TT> the resulting crack edges.
02221 
02222     \code
02223 sign of difference image         insert cracks         resulting CrackEdgeImage
02224 
02225                                    + . - . -              . * . . .
02226       + - -                        . . . . .              . * * * .
02227       + + -               =>       + . + . -      =>      . . . * .
02228       + + +                        . . . . .              . . . * *
02229                                    + . + . +              . . . . .
02230     \endcode
02231 
02232 Starting from the original binary image (left), we insert crack pixels
02233 to get to the double-sized image (center). Finally, we mark all
02234 crack pixels whose non-crack neighbors have different signs as
02235 crack edge points, while all other pixels (crack and non-crack) become
02236 region pixels.
02237 
02238 <b>Requirements on a Crack Edge Image:</b>
02239 
02240 <ul>
02241     <li>Crack Edge Images have odd width and height.
02242     <li>Crack pixels have at least one odd coordinate.
02243     <li>Only crack pixels may be marked as edge points.
02244     <li>Crack pixels with two odd coordinates must be marked as edge points
02245         whenever any of their neighboring crack pixels was marked.
02246 </ul>
02247 
02248 The last two requirements ensure that both edges and regions are 4-connected.
02249 Thus, 4-connectivity and 8-connectivity yield identical connected
02250 components in a Crack Edge Image (so called <i>well-composedness</i>).
02251 This ensures that Crack Edge Images have nice topological properties
02252 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000).
02253 */
02254 
02255 
02256 } // namespace vigra
02257 
02258 #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.0 (Thu Aug 25 2011)