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

details vigra/edgedetection.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)