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

vigra/edgedetection.hxx

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

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

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)