[ 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.4.0, Dec 21 2005 )                                    */
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 "vigra/utilities.hxx"
00046 #include "vigra/numerictraits.hxx"
00047 #include "vigra/stdimage.hxx"
00048 #include "vigra/stdimagefunctions.hxx"
00049 #include "vigra/recursiveconvolution.hxx"
00050 #include "vigra/separableconvolution.hxx"
00051 #include "vigra/labelimage.hxx"
00052 #include "vigra/mathutil.hxx"
00053 #include "vigra/pixelneighborhood.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         inline 
00125         void differenceOfExponentialEdgeImage(
00126                triple<SrcIterator, SrcIterator, SrcAccessor> src,
00127                pair<DestIterator, DestAccessor> dest,
00128                double scale, GradValue gradient_threshold,
00129                DestValue edge_marker = NumericTraits<DestValue>::one())
00130     }
00131     \endcode
00132     
00133     <b> Usage:</b>
00134     
00135         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00136     Namespace: vigra
00137     
00138     \code
00139     vigra::BImage src(w,h), edges(w,h);
00140     
00141     // empty edge image
00142     edges = 0;
00143     ...
00144     
00145     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00146     vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 
00147                                      0.8, 4.0, 1);
00148     \endcode
00149 
00150     <b> Required Interface:</b>
00151     
00152     \code
00153     SrcImageIterator src_upperleft, src_lowerright;
00154     DestImageIterator dest_upperleft;
00155     
00156     SrcAccessor src_accessor;
00157     DestAccessor dest_accessor;
00158     
00159     SrcAccessor::value_type u = src_accessor(src_upperleft);
00160     double d;
00161     GradValue gradient_threshold;
00162     
00163     u = u + u
00164     u = u - u
00165     u = u * u
00166     u = d * u
00167     u < gradient_threshold
00168     
00169     DestValue edge_marker;
00170     dest_accessor.set(edge_marker, dest_upperleft);
00171     \endcode
00172     
00173     <b> Preconditions:</b>
00174     
00175     \code
00176     scale > 0
00177     gradient_threshold > 0
00178     \endcode
00179 */
00180 template <class SrcIterator, class SrcAccessor,
00181           class DestIterator, class DestAccessor, 
00182           class GradValue, class DestValue>
00183 void differenceOfExponentialEdgeImage(
00184                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00185            DestIterator dul, DestAccessor da,
00186            double scale, GradValue gradient_threshold, DestValue edge_marker)
00187 {
00188     vigra_precondition(scale > 0,
00189                  "differenceOfExponentialEdgeImage(): scale > 0 required.");
00190          
00191     vigra_precondition(gradient_threshold > 0,
00192                  "differenceOfExponentialEdgeImage(): "
00193          "gradient_threshold > 0 required.");
00194          
00195     int w = slr.x - sul.x;
00196     int h = slr.y - sul.y;
00197     int x,y;
00198 
00199     typedef typename 
00200         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00201     TMPTYPE;
00202     typedef BasicImage<TMPTYPE> TMPIMG;
00203 
00204     TMPIMG tmp(w,h);
00205     TMPIMG smooth(w,h);
00206     
00207     recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
00208     recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
00209 
00210     recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
00211     recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
00212 
00213     typename TMPIMG::Iterator iy = smooth.upperLeft();
00214     typename TMPIMG::Iterator ty = tmp.upperLeft();
00215     DestIterator              dy = dul;
00216     
00217     static const Diff2D right(1, 0);
00218     static const Diff2D bottom(0, 1);
00219     
00220     
00221     TMPTYPE thresh = (gradient_threshold * gradient_threshold) * 
00222                      NumericTraits<TMPTYPE>::one();
00223     TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
00224 
00225     for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, ++dy.y)
00226     {
00227         typename TMPIMG::Iterator ix = iy;
00228         typename TMPIMG::Iterator tx = ty;
00229         DestIterator              dx = dy;
00230 
00231         for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
00232         {
00233             TMPTYPE diff = *tx - *ix;
00234             TMPTYPE gx = tx[right] - *tx;
00235             TMPTYPE gy = tx[bottom] - *tx;
00236 
00237             if((gx * gx > thresh) &&
00238                 (diff * (tx[right] - ix[right]) < zero))
00239             {
00240                 if(gx < zero)
00241                 {
00242                     da.set(edge_marker, dx, right);
00243                 }
00244                 else
00245                 {
00246                     da.set(edge_marker, dx);
00247                 }
00248             }
00249             if(((gy * gy > thresh) &&
00250                 (diff * (tx[bottom] - ix[bottom]) < zero)))
00251             {
00252                 if(gy < zero)
00253                 {
00254                     da.set(edge_marker, dx, bottom);
00255                 }
00256                 else
00257                 {
00258                     da.set(edge_marker, dx);
00259                 }
00260             }
00261         }
00262     }
00263     
00264     typename TMPIMG::Iterator ix = iy;
00265     typename TMPIMG::Iterator tx = ty;
00266     DestIterator              dx = dy;
00267     
00268     for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, ++dx.x)
00269     {
00270         TMPTYPE diff = *tx - *ix;
00271         TMPTYPE gx = tx[right] - *tx;
00272 
00273         if((gx * gx > thresh) &&
00274            (diff * (tx[right] - ix[right]) < zero))
00275         {
00276             if(gx < zero)
00277             {
00278                 da.set(edge_marker, dx, right);
00279             }
00280             else
00281             {
00282                 da.set(edge_marker, dx);
00283             }
00284         }
00285     }
00286 }
00287 
00288 template <class SrcIterator, class SrcAccessor,
00289           class DestIterator, class DestAccessor,
00290           class GradValue>
00291 inline 
00292 void differenceOfExponentialEdgeImage(
00293                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00294            DestIterator dul, DestAccessor da,
00295            double scale, GradValue gradient_threshold)
00296 {
00297     differenceOfExponentialEdgeImage(sul, slr, sa, dul, da, 
00298                                         scale, gradient_threshold, 1);
00299 }
00300 
00301 template <class SrcIterator, class SrcAccessor,
00302           class DestIterator, class DestAccessor, 
00303       class GradValue, class DestValue>
00304 inline 
00305 void differenceOfExponentialEdgeImage(
00306            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00307        pair<DestIterator, DestAccessor> dest,
00308        double scale, GradValue gradient_threshold,
00309        DestValue edge_marker)
00310 {
00311     differenceOfExponentialEdgeImage(src.first, src.second, src.third,
00312                                         dest.first, dest.second,
00313                     scale, gradient_threshold,
00314                     edge_marker);
00315 }
00316 
00317 template <class SrcIterator, class SrcAccessor,
00318           class DestIterator, class DestAccessor,
00319           class GradValue>
00320 inline 
00321 void differenceOfExponentialEdgeImage(
00322            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00323        pair<DestIterator, DestAccessor> dest,
00324        double scale, GradValue gradient_threshold)
00325 {
00326     differenceOfExponentialEdgeImage(src.first, src.second, src.third,
00327                                         dest.first, dest.second,
00328                     scale, gradient_threshold, 1);
00329 }
00330 
00331 /********************************************************/
00332 /*                                                      */
00333 /*        differenceOfExponentialCrackEdgeImage         */
00334 /*                                                      */
00335 /********************************************************/
00336 
00337 /** \brief Detect and mark edges in a crack edge image using the Shen/Castan zero-crossing detector.
00338 
00339     This operator applies an exponential filter to the source image 
00340     at the given <TT>scale</TT> and subtracts the result from the original image. 
00341     Zero crossings are detected in the resulting difference image. Whenever the
00342     gradient at a zero crossing is greater than the given <TT>gradient_threshold</TT>,
00343     an edge point is marked (using <TT>edge_marker</TT>) in the destination image 
00344     <i>between</i> the corresponding original pixels. Topologically, this means we 
00345     must insert additional pixels between the original ones to represent the
00346     boundaries between the pixels (the so called zero- and one-cells, with the original
00347     pixels being two-cells). Within VIGRA, such an image is called \ref CrackEdgeImage.
00348     To allow insertion of the zero- and one-cells, the destination image must have twice the 
00349     size of the original (precisely, <TT>(2*w-1)</TT> by <TT>(2*h-1)</TT> pixels). Then the algorithm 
00350     proceeds as follows:
00351     
00352     \code
00353 sign of difference image     insert zero- and one-cells     resulting edge points (*)
00354 
00355                                      + . - . -                   . * . . .
00356       + - -                          . . . . .                   . * * * .
00357       + + -               =>         + . + . -           =>      . . . * .
00358       + + +                          . . . . .                   . . . * *
00359                                      + . + . +                   . . . . .
00360     \endcode
00361     
00362     Thus the edge points are marked where they actually are - in between the pixels. 
00363     An important property of the resulting edge image is that it conforms to the notion 
00364     of well-composedness as defined by Latecki et al., i.e. connected regions and edges 
00365     obtained by a subsequent \ref Labeling do not depend on 
00366     whether 4- or 8-connectivity is used.
00367     The non-edge pixels (<TT>.</TT>) in the destination image remain unchanged. 
00368     The result conformes to the requirements of a \ref CrackEdgeImage. It can be further
00369     improved by the post-processing operations \ref removeShortEdges() and
00370     \ref closeGapsInCrackEdgeImage().
00371     
00372     The source value type (<TT>SrcAccessor::value_type</TT>) must be a linear algebra, i.e. addition, 
00373     subtraction and multiplication of the type with itself, and multiplication 
00374     with double and 
00375     \ref NumericTraits "NumericTraits" must 
00376     be defined. In addition, this type must be less-comparable.
00377     
00378     <b> Declarations:</b>
00379     
00380     pass arguments explicitly:
00381     \code
00382     namespace vigra {
00383         template <class SrcIterator, class SrcAccessor,
00384               class DestIterator, class DestAccessor,
00385               class GradValue,
00386               class DestValue = DestAccessor::value_type>
00387         void differenceOfExponentialCrackEdgeImage(
00388                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00389                DestIterator dul, DestAccessor da,
00390                double scale, GradValue gradient_threshold, 
00391                DestValue edge_marker = NumericTraits<DestValue>::one())
00392     }
00393     \endcode
00394     
00395     use argument objects in conjunction with \ref ArgumentObjectFactories:
00396     \code
00397     namespace vigra {
00398         template <class SrcIterator, class SrcAccessor,
00399               class DestIterator, class DestAccessor, 
00400               class GradValue,
00401               class DestValue = DestAccessor::value_type>
00402         inline 
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 template <class SrcIterator, class SrcAccessor,
00465           class DestIterator, class DestAccessor,
00466           class GradValue, class DestValue>
00467 void differenceOfExponentialCrackEdgeImage(
00468                SrcIterator sul, SrcIterator slr, SrcAccessor sa,
00469            DestIterator dul, DestAccessor da,
00470            double scale, GradValue gradient_threshold, 
00471            DestValue edge_marker)
00472 {
00473     vigra_precondition(scale > 0,
00474                  "differenceOfExponentialCrackEdgeImage(): scale > 0 required.");
00475          
00476     vigra_precondition(gradient_threshold > 0,
00477                  "differenceOfExponentialCrackEdgeImage(): "
00478          "gradient_threshold > 0 required.");
00479          
00480     int w = slr.x - sul.x;
00481     int h = slr.y - sul.y;
00482     int x, y;
00483 
00484     typedef typename 
00485         NumericTraits<typename SrcAccessor::value_type>::RealPromote
00486     TMPTYPE;
00487     typedef BasicImage<TMPTYPE> TMPIMG;
00488 
00489     TMPIMG tmp(w,h);
00490     TMPIMG smooth(w,h);
00491     
00492     TMPTYPE zero = NumericTraits<TMPTYPE>::zero();
00493     
00494     static const Diff2D right(1,0);
00495     static const Diff2D bottom(0,1);
00496     static const Diff2D left(-1,0);
00497     static const Diff2D top(0,-1);
00498     
00499     recursiveSmoothX(srcIterRange(sul, slr, sa), destImage(tmp), scale / 2.0);
00500     recursiveSmoothY(srcImageRange(tmp), destImage(tmp), scale / 2.0);
00501 
00502     recursiveSmoothX(srcImageRange(tmp), destImage(smooth), scale);
00503     recursiveSmoothY(srcImageRange(smooth), destImage(smooth), scale);
00504 
00505     typename TMPIMG::Iterator iy = smooth.upperLeft();
00506     typename TMPIMG::Iterator ty = tmp.upperLeft();
00507     DestIterator              dy = dul;
00508     
00509     TMPTYPE thresh = (gradient_threshold * gradient_threshold) * 
00510                      NumericTraits<TMPTYPE>::one();
00511 
00512     // find zero crossings above threshold
00513     for(y=0; y<h-1; ++y, ++iy.y, ++ty.y, dy.y+=2)
00514     {
00515         typename TMPIMG::Iterator ix = iy;
00516         typename TMPIMG::Iterator tx = ty;
00517         DestIterator              dx = dy;
00518 
00519         for(int x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
00520         {
00521             TMPTYPE diff = *tx - *ix;
00522             TMPTYPE gx = tx[right] - *tx;
00523             TMPTYPE gy = tx[bottom] - *tx;
00524 
00525             if((gx * gx > thresh) &&
00526                (diff * (tx[right] - ix[right]) < zero))
00527             {
00528                 da.set(edge_marker, dx, right);
00529             }
00530             if((gy * gy > thresh) &&
00531                (diff * (tx[bottom] - ix[bottom]) < zero))
00532             {
00533                 da.set(edge_marker, dx, bottom);
00534             }
00535         }
00536 
00537         TMPTYPE diff = *tx - *ix;
00538         TMPTYPE gy = tx[bottom] - *tx;
00539 
00540         if((gy * gy > thresh) &&
00541            (diff * (tx[bottom] - ix[bottom]) < zero))
00542         {
00543             da.set(edge_marker, dx, bottom);
00544         }
00545     }
00546     
00547     typename TMPIMG::Iterator ix = iy;
00548     typename TMPIMG::Iterator tx = ty;
00549     DestIterator              dx = dy;
00550     
00551     for(x=0; x<w-1; ++x, ++ix.x, ++tx.x, dx.x+=2)
00552     {
00553         TMPTYPE diff = *tx - *ix;
00554         TMPTYPE gx = tx[right] - *tx;
00555 
00556         if((gx * gx > thresh) &&
00557            (diff * (tx[right] - ix[right]) < zero))
00558         {
00559             da.set(edge_marker, dx, right);
00560         }
00561     }
00562 
00563     iy = smooth.upperLeft() + Diff2D(0,1);
00564     ty = tmp.upperLeft() + Diff2D(0,1);
00565     dy = dul + Diff2D(1,2);
00566     
00567     static const Diff2D topleft(-1,-1);
00568     static const Diff2D topright(1,-1);
00569     static const Diff2D bottomleft(-1,1);
00570     static const Diff2D bottomright(1,1);
00571 
00572     // find missing 1-cells below threshold (x-direction)
00573     for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
00574     {
00575         typename TMPIMG::Iterator ix = iy;
00576         typename TMPIMG::Iterator tx = ty;
00577         DestIterator              dx = dy;
00578 
00579         for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
00580         {
00581             if(da(dx) == edge_marker) continue;
00582 
00583             TMPTYPE diff = *tx - *ix;
00584 
00585             if((diff * (tx[right] - ix[right]) < zero) &&
00586                (((da(dx, bottomright) == edge_marker) && 
00587                  (da(dx, topleft) == edge_marker)) ||
00588                 ((da(dx, bottomleft) == edge_marker) && 
00589                  (da(dx, topright) == edge_marker))))
00590 
00591             {
00592                 da.set(edge_marker, dx);
00593             }
00594         }
00595     }
00596     
00597     iy = smooth.upperLeft() + Diff2D(1,0);
00598     ty = tmp.upperLeft() + Diff2D(1,0);
00599     dy = dul + Diff2D(2,1);
00600 
00601     // find missing 1-cells below threshold (y-direction)
00602     for(y=0; y<h-2; ++y, ++iy.y, ++ty.y, dy.y+=2)
00603     {
00604         typename TMPIMG::Iterator ix = iy;
00605         typename TMPIMG::Iterator tx = ty;
00606         DestIterator              dx = dy;
00607 
00608         for(int x=0; x<w-2; ++x, ++ix.x, ++tx.x, dx.x+=2)
00609         {
00610             if(da(dx) == edge_marker) continue;
00611 
00612             TMPTYPE diff = *tx - *ix;
00613 
00614             if((diff * (tx[bottom] - ix[bottom]) < zero) &&
00615                (((da(dx, bottomright) == edge_marker) && 
00616                  (da(dx, topleft) == edge_marker)) ||
00617                 ((da(dx, bottomleft) == edge_marker) && 
00618                  (da(dx, topright) == edge_marker))))
00619 
00620             {
00621                 da.set(edge_marker, dx);
00622             }
00623         }
00624     }
00625     
00626     dy = dul + Diff2D(1,1);
00627 
00628     // find missing 0-cells 
00629     for(y=0; y<h-1; ++y, dy.y+=2)
00630     {
00631         DestIterator              dx = dy;
00632 
00633         for(int x=0; x<w-1; ++x, dx.x+=2)
00634         {
00635             static const Diff2D dist[] = {right, top, left, bottom };
00636 
00637             int i;
00638             for(i=0; i<4; ++i)
00639             {
00640             if(da(dx, dist[i]) == edge_marker) break;
00641             }
00642 
00643             if(i < 4) da.set(edge_marker, dx);
00644         }
00645     }
00646 }
00647 
00648 template <class SrcIterator, class SrcAccessor,
00649           class DestIterator, class DestAccessor,
00650           class GradValue, class DestValue>
00651 inline 
00652 void differenceOfExponentialCrackEdgeImage(
00653            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00654        pair<DestIterator, DestAccessor> dest,
00655        double scale, GradValue gradient_threshold,
00656        DestValue edge_marker)
00657 {
00658     differenceOfExponentialCrackEdgeImage(src.first, src.second, src.third,
00659                                         dest.first, dest.second,
00660                     scale, gradient_threshold,
00661                     edge_marker);
00662 }
00663 
00664 /********************************************************/
00665 /*                                                      */
00666 /*                  removeShortEdges                    */
00667 /*                                                      */
00668 /********************************************************/
00669 
00670 /** \brief Remove short edges from an edge image.
00671 
00672     This algorithm can be applied as a post-processing operation of 
00673     \ref differenceOfExponentialEdgeImage() and \ref differenceOfExponentialCrackEdgeImage(). 
00674     It removes all edges that are shorter than <TT>min_edge_length</TT>. The corresponding
00675     pixels are set to the <TT>non_edge_marker</TT>. The idea behind this algorithms is
00676     that very short edges are probably caused by noise and don't represent interesting
00677     image structure. Technically, the algorithms executes a connected components labeling,
00678     so the image's value type must be equality comparable. 
00679     
00680     If the source image fulfills the requirements of a \ref CrackEdgeImage,
00681     it will still do so after application of this algorithm.
00682     
00683     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 
00684     i.e. on only one image. Also, the algorithm assumes that all non-edges pixels are already
00685     marked with the given <TT>non_edge_marker</TT> value.
00686     
00687     <b> Declarations:</b>
00688     
00689     pass arguments explicitly:
00690     \code
00691     namespace vigra {
00692         template <class Iterator, class Accessor, class SrcValue>
00693         void removeShortEdges(
00694                Iterator sul, Iterator slr, Accessor sa,
00695                int min_edge_length, SrcValue non_edge_marker)
00696     }
00697     \endcode
00698     
00699     use argument objects in conjunction with \ref ArgumentObjectFactories:
00700     \code
00701     namespace vigra {
00702         template <class Iterator, class Accessor, class SrcValue>
00703         inline 
00704         void removeShortEdges(
00705                triple<Iterator, Iterator, Accessor> src,
00706                int min_edge_length, SrcValue non_edge_marker)
00707     }
00708     \endcode
00709     
00710     <b> Usage:</b>
00711     
00712         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00713     Namespace: vigra
00714     
00715     \code
00716     vigra::BImage src(w,h), edges(w,h);
00717     
00718     // empty edge image
00719     edges = 0;
00720     ...
00721     
00722     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00723     vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 
00724                                      0.8, 4.0, 1);
00725                     
00726     // zero edges shorter than 10 pixels
00727     vigra::removeShortEdges(srcImageRange(edges), 10, 0);
00728     \endcode
00729 
00730     <b> Required Interface:</b>
00731     
00732     \code
00733     SrcImageIterator src_upperleft, src_lowerright;
00734     DestImageIterator dest_upperleft;
00735     
00736     SrcAccessor src_accessor;
00737     DestAccessor dest_accessor;
00738     
00739     SrcAccessor::value_type u = src_accessor(src_upperleft);
00740 
00741     u == u
00742     
00743     SrcValue non_edge_marker;
00744     src_accessor.set(non_edge_marker, src_upperleft);
00745     \endcode
00746 */
00747 template <class Iterator, class Accessor, class Value>
00748 void removeShortEdges(
00749                Iterator sul, Iterator slr, Accessor sa,
00750            unsigned int min_edge_length, Value non_edge_marker)
00751 {
00752     int w = slr.x - sul.x;
00753     int h = slr.y - sul.y;
00754     int x,y;
00755 
00756     IImage labels(w, h);
00757     labels = 0;
00758 
00759     int number_of_regions = 
00760                 labelImageWithBackground(srcIterRange(sul,slr,sa), 
00761                                      destImage(labels), true, non_edge_marker);
00762     
00763     ArrayOfRegionStatistics<FindROISize<int> > 
00764                                          region_stats(number_of_regions);
00765     
00766     inspectTwoImages(srcImageRange(labels), srcImage(labels), region_stats);
00767              
00768     IImage::Iterator ly = labels.upperLeft();
00769     Iterator oy = sul;
00770     
00771     for(y=0; y<h; ++y, ++oy.y, ++ly.y)
00772     {
00773         Iterator ox(oy);
00774         IImage::Iterator lx(ly);
00775 
00776         for(x=0; x<w; ++x, ++ox.x, ++lx.x)
00777         {
00778             if(sa(ox) == non_edge_marker) continue;
00779             if((region_stats[*lx].count) < min_edge_length)
00780             {
00781                  sa.set(non_edge_marker, ox);
00782             }
00783         }
00784     }
00785 }
00786 
00787 template <class Iterator, class Accessor, class Value>
00788 inline 
00789 void removeShortEdges(
00790            triple<Iterator, Iterator, Accessor> src,
00791        unsigned int min_edge_length, Value non_edge_marker)
00792 {
00793     removeShortEdges(src.first, src.second, src.third,
00794                      min_edge_length, non_edge_marker);
00795 }
00796 
00797 /********************************************************/
00798 /*                                                      */
00799 /*             closeGapsInCrackEdgeImage                */
00800 /*                                                      */
00801 /********************************************************/
00802 
00803 /** \brief Close one-pixel wide gaps in a cell grid edge image.
00804 
00805     This algorithm is typically applied as a post-processing operation of 
00806     \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
00807     the requirements of a \ref CrackEdgeImage, and will still do so after 
00808     application of this algorithm.
00809 
00810     It closes one pixel wide gaps in the edges resulting from this algorithm. 
00811     Since these gaps are usually caused by zero crossing slightly below the gradient 
00812     threshold used in edge detection, this algorithms acts like a weak hysteresis 
00813     thresholding. The newly found edge pixels are marked with the given <TT>edge_marker</TT>.
00814     The image's value type must be equality comparable. 
00815     
00816     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 
00817     i.e. on only one image.
00818     
00819     <b> Declarations:</b>
00820     
00821     pass arguments explicitly:
00822     \code
00823     namespace vigra {
00824         template <class SrcIterator, class SrcAccessor, class SrcValue>
00825         void closeGapsInCrackEdgeImage(
00826                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
00827                SrcValue edge_marker)
00828     }
00829     \endcode
00830     
00831     use argument objects in conjunction with \ref ArgumentObjectFactories:
00832     \code
00833     namespace vigra {
00834         template <class SrcIterator, class SrcAccessor, class SrcValue>
00835         inline
00836         void closeGapsInCrackEdgeImage(
00837                triple<SrcIterator, SrcIterator, SrcAccessor> src,
00838                SrcValue edge_marker)
00839     }
00840     \endcode
00841     
00842     <b> Usage:</b>
00843     
00844         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
00845     Namespace: vigra
00846     
00847     \code
00848     vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
00849     
00850     // empty edge image
00851     edges = 0;
00852     ...
00853     
00854     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
00855     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 
00856                                          0.8, 4.0, 1);
00857                     
00858     // close gaps, mark with 1
00859     vigra::closeGapsInCrackEdgeImage(srcImageRange(edges), 1);
00860                     
00861     // zero edges shorter than 20 pixels
00862     vigra::removeShortEdges(srcImageRange(edges), 10, 0);
00863     \endcode
00864 
00865     <b> Required Interface:</b>
00866     
00867     \code
00868     SrcImageIterator src_upperleft, src_lowerright;
00869     
00870     SrcAccessor src_accessor;
00871     DestAccessor dest_accessor;
00872     
00873     SrcAccessor::value_type u = src_accessor(src_upperleft);
00874 
00875     u == u
00876     u != u
00877     
00878     SrcValue edge_marker;
00879     src_accessor.set(edge_marker, src_upperleft);
00880     \endcode
00881 */
00882 template <class SrcIterator, class SrcAccessor, class SrcValue>
00883 void closeGapsInCrackEdgeImage(
00884                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
00885            SrcValue edge_marker)
00886 {
00887     int w = (slr.x - sul.x) / 2;
00888     int h = (slr.y - sul.y) / 2;
00889     int x, y;
00890 
00891     int count1, count2, count3;
00892 
00893     static const Diff2D right(1,0);
00894     static const Diff2D bottom(0,1);
00895     static const Diff2D left(-1,0);
00896     static const Diff2D top(0,-1);
00897     
00898     static const Diff2D leftdist[] = { 
00899         Diff2D(0, 0), Diff2D(-1, 1), Diff2D(-2, 0), Diff2D(-1, -1)};
00900     static const Diff2D rightdist[] = { 
00901         Diff2D(2, 0), Diff2D(1, 1), Diff2D(0, 0), Diff2D(1, -1)};
00902     static const Diff2D topdist[] = { 
00903         Diff2D(1, -1), Diff2D(0, 0), Diff2D(-1, -1), Diff2D(0, -2)};
00904     static const Diff2D bottomdist[] = { 
00905         Diff2D(1, 1), Diff2D(0, 2), Diff2D(-1, 1), Diff2D(0, 0)};
00906 
00907     int i;
00908 
00909     SrcIterator sy = sul + Diff2D(0,1);
00910     SrcIterator sx;
00911     
00912     // close 1-pixel wide gaps (x-direction)
00913     for(y=0; y<h; ++y, sy.y+=2)
00914     {
00915         sx = sy + Diff2D(2,0);
00916 
00917         for(x=2; x<w; ++x, sx.x+=2)
00918         {
00919             if(sa(sx) == edge_marker) continue;
00920 
00921             if(sa(sx, left) != edge_marker) continue;
00922             if(sa(sx, right) != edge_marker) continue;
00923 
00924             count1 = 0;
00925             count2 = 0;
00926             count3 = 0;
00927 
00928             for(i=0; i<4; ++i)
00929             {
00930                 if(sa(sx, leftdist[i]) == edge_marker) 
00931                 {
00932                     ++count1;
00933                     count3 ^= 1 << i;
00934                 }
00935                 if(sa(sx, rightdist[i]) == edge_marker) 
00936                 {
00937                     ++count2;
00938                     count3 ^= 1 << i;
00939                 }
00940             }
00941 
00942             if(count1 <= 1 || count2 <= 1 || count3 == 15) 
00943             {
00944                 sa.set(edge_marker, sx);
00945             }
00946         }
00947    }
00948     
00949     sy = sul + Diff2D(1,2);
00950 
00951     // close 1-pixel wide gaps (y-direction)
00952     for(y=2; y<h; ++y, sy.y+=2)
00953     {
00954         sx = sy;
00955 
00956         for(x=0; x<w; ++x, sx.x+=2)
00957         {
00958             if(sa(sx) == edge_marker) continue;
00959 
00960             if(sa(sx, top) != edge_marker) continue;
00961             if(sa(sx, bottom) != edge_marker) continue;
00962 
00963             count1 = 0;
00964             count2 = 0;
00965             count3 = 0;
00966 
00967             for(i=0; i<4; ++i)
00968             {
00969                 if(sa(sx, topdist[i]) == edge_marker)
00970                 {
00971                     ++count1;
00972                     count3 ^= 1 << i;
00973                 }
00974                 if(sa(sx, bottomdist[i]) == edge_marker)
00975                 {
00976                     ++count2;
00977                     count3 ^= 1 << i;
00978                 }
00979             }
00980 
00981             if(count1 <= 1 || count2 <= 1 || count3 == 15)
00982             {
00983                 sa.set(edge_marker, sx);
00984             }
00985         }
00986     }
00987 }
00988 
00989 template <class SrcIterator, class SrcAccessor, class SrcValue>
00990 inline
00991 void closeGapsInCrackEdgeImage(
00992            triple<SrcIterator, SrcIterator, SrcAccessor> src,
00993        SrcValue edge_marker)
00994 {
00995     closeGapsInCrackEdgeImage(src.first, src.second, src.third,
00996                     edge_marker);
00997 }
00998 
00999 /********************************************************/
01000 /*                                                      */
01001 /*              beautifyCrackEdgeImage                  */
01002 /*                                                      */
01003 /********************************************************/
01004 
01005 /** \brief Beautify crack edge image for visualization.
01006 
01007     This algorithm is applied as a post-processing operation of 
01008     \ref differenceOfExponentialCrackEdgeImage(). The source image must fulfill
01009     the requirements of a \ref CrackEdgeImage, but will <b> not</b> do so after 
01010     application of this algorithm. In particular, the algorithm removes zero-cells 
01011     marked as edges to avoid staircase effects on diagonal lines like this:
01012     
01013     \code
01014     original edge points (*)     resulting edge points
01015     
01016           . * . . .                   . * . . .
01017           . * * * .                   . . * . .
01018           . . . * .           =>      . . . * .
01019           . . . * *                   . . . . *
01020           . . . . .                   . . . . .
01021     \endcode
01022     
01023     Therfore, this algorithm should only be applied as a vizualization aid, i.e. 
01024     for human inspection. The algorithm assumes that edges are marked with <TT>edge_marker</TT>, 
01025     and background pixels with <TT>background_marker</TT>. The image's value type must be 
01026     equality comparable. 
01027     
01028     Note that this algorithm, unlike most other algorithms in VIGRA, operates in-place, 
01029     i.e. on only one image.
01030     
01031     <b> Declarations:</b>
01032     
01033     pass arguments explicitly:
01034     \code
01035     namespace vigra {
01036         template <class SrcIterator, class SrcAccessor, class SrcValue>
01037         void beautifyCrackEdgeImage(
01038                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
01039                SrcValue edge_marker, SrcValue background_marker)
01040     }
01041     \endcode
01042     
01043     use argument objects in conjunction with \ref ArgumentObjectFactories:
01044     \code
01045     namespace vigra {
01046         template <class SrcIterator, class SrcAccessor, class SrcValue>
01047         inline
01048         void beautifyCrackEdgeImage(
01049                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01050                SrcValue edge_marker, SrcValue background_marker)
01051     }
01052     \endcode
01053     
01054     <b> Usage:</b>
01055     
01056         <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01057     Namespace: vigra
01058     
01059     \code
01060     vigra::BImage src(w,h), edges(2*w-1, 2*h-1);
01061     
01062     // empty edge image
01063     edges = 0;
01064     ...
01065     
01066     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
01067     vigra::differenceOfExponentialCrackEdgeImage(srcImageRange(src), destImage(edges), 
01068                                          0.8, 4.0, 1);
01069                     
01070     // beautify edge image for visualization
01071     vigra::beautifyCrackEdgeImage(destImageRange(edges), 1, 0);
01072     
01073     // show to the user
01074     window.open(edges);
01075     \endcode
01076 
01077     <b> Required Interface:</b>
01078     
01079     \code
01080     SrcImageIterator src_upperleft, src_lowerright;
01081     
01082     SrcAccessor src_accessor;
01083     DestAccessor dest_accessor;
01084     
01085     SrcAccessor::value_type u = src_accessor(src_upperleft);
01086 
01087     u == u
01088     u != u
01089     
01090     SrcValue background_marker;
01091     src_accessor.set(background_marker, src_upperleft);
01092     \endcode
01093 */
01094 template <class SrcIterator, class SrcAccessor, class SrcValue>
01095 void beautifyCrackEdgeImage(
01096                SrcIterator sul, SrcIterator slr, SrcAccessor sa, 
01097            SrcValue edge_marker, SrcValue background_marker)
01098 {
01099     int w = (slr.x - sul.x) / 2;
01100     int h = (slr.y - sul.y) / 2;
01101     int x, y;
01102 
01103     SrcIterator sy = sul + Diff2D(1,1);
01104     SrcIterator sx;
01105     
01106     static const Diff2D right(1,0);
01107     static const Diff2D bottom(0,1);
01108     static const Diff2D left(-1,0);
01109     static const Diff2D top(0,-1);
01110     
01111     //  delete 0-cells at corners
01112     for(y=0; y<h; ++y, sy.y+=2)
01113     {
01114         sx = sy;
01115 
01116         for(x=0; x<w; ++x, sx.x+=2)
01117         {
01118             if(sa(sx) != edge_marker) continue;
01119 
01120             if(sa(sx, right) == edge_marker && sa(sx, left) == edge_marker) continue;
01121             if(sa(sx, bottom) == edge_marker && sa(sx, top) == edge_marker) continue;
01122 
01123             sa.set(background_marker, sx);
01124         }
01125     }
01126 }
01127 
01128 template <class SrcIterator, class SrcAccessor, class SrcValue>
01129 inline
01130 void beautifyCrackEdgeImage(
01131            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01132            SrcValue edge_marker, SrcValue background_marker)
01133 {
01134     beautifyCrackEdgeImage(src.first, src.second, src.third,
01135                     edge_marker, background_marker);
01136 }
01137 
01138 
01139 /** Helper class that stores edgel attributes.
01140 */
01141 class Edgel
01142 {
01143   public:
01144         /** The edgel's sub-pixel x coordinate.
01145         */    
01146     float x; 
01147 
01148         /** The edgel's sub-pixel y coordinate.
01149         */    
01150     float y;
01151 
01152         /** The edgel's strength (magnitude of the gradient vector).
01153         */    
01154     float strength;
01155     
01156         /**
01157         The edgel's orientation. This is the angle 
01158         between the x-axis and the edge, so that the bright side of the 
01159         edge is on the right. The angle is measured
01160         counter-clockwise in radians like this:
01161 
01162 
01163         \code
01164 
01165   edgel axis
01166       \  (bright side)
01167  (dark \
01168  side)  \ /__ 
01169          \\  \ orientation angle
01170           \  |
01171            +------------> x-axis
01172            |
01173            |
01174            |
01175            |
01176     y-axis V
01177         \endcode
01178 
01179         So, for example a vertical edge with its dark side on the left
01180         has orientation PI/2, and a horizontal edge with dark side on top
01181         has orientation 0. Obviously, the edge's orientation changes
01182         by PI if the contrast is reversed.
01183 
01184         */
01185     float orientation;
01186     
01187     Edgel()
01188     : x(0.0f), y(0.0f), strength(0.0f), orientation(0.0f)
01189     {}
01190     
01191     Edgel(float ix, float iy, float is, float io)
01192     : x(ix), y(iy), strength(is), orientation(io)
01193     {}
01194 };
01195 
01196 template <class Image1, class Image2>
01197 void internalCannyFindEdgels(Image1 const & dx,
01198                              Image1 const & dy,
01199                              Image2 const & magnitude,
01200                              std::vector<Edgel> & edgels)
01201 {
01202     typedef typename Image1::value_type PixelType;
01203     
01204     PixelType zero = NumericTraits<PixelType>::zero();
01205     double tan22_5 = M_SQRT2 - 1.0;
01206     
01207     for(int y=1; y<dx.height()-1; ++y)
01208     {
01209         for(int x=1; x<dx.width()-1; ++x)
01210         {
01211             bool maximum_found = false;
01212             Edgel edgel;
01213             
01214             PixelType gradx = dx(x,y);
01215             PixelType grady = dy(x,y);
01216             
01217             // find out quadrant
01218             if(abs(grady) < tan22_5*abs(gradx))
01219             {
01220                 // north-south edge
01221                 PixelType m1 = magnitude(x-1, y);
01222                 PixelType m2 = magnitude(x, y);
01223                 PixelType m3 = magnitude(x+1, y);
01224                 
01225                 if(m1 < m2 && m3 <= m2)
01226                 {
01227                     edgel.y = y;
01228                 
01229                     // local maximum => quadratic interpolation of sub-pixel location
01230                     PixelType del = (m1 - m3) / 2.0;
01231                     del /= (m1 + m3 - 2.0*m2);
01232                     edgel.x = x + del;
01233                     edgel.strength = m2;
01234                     
01235                     maximum_found = true;                    
01236                 }
01237             }
01238             else if(abs(gradx) < tan22_5*abs(grady))
01239             {
01240                 // west-east edge
01241                 PixelType m1 = magnitude(x, y-1);
01242                 PixelType m2 = magnitude(x, y);
01243                 PixelType m3 = magnitude(x, y+1);
01244                 
01245                 if(m1 < m2 && m3 <= m2)
01246                 {
01247                     edgel.x = x;
01248                 
01249                     // local maximum => quadratic interpolation of sub-pixel location
01250                     PixelType del = (m1 - m3) / 2.0;
01251                     del /= (m1 + m3 - 2.0*m2);
01252                     edgel.y = y + del;
01253                     edgel.strength = m2;
01254                     
01255                     maximum_found = true;                    
01256                 }
01257             }
01258             else if(gradx*grady < zero)
01259             {
01260                 // north-west-south-east edge
01261                 PixelType m1 = magnitude(x+1, y-1);
01262                 PixelType m2 = magnitude(x, y);
01263                 PixelType m3 = magnitude(x-1, y+1);
01264                 
01265                 if(m1 < m2 && m3 <= m2)
01266                 {
01267                     // local maximum => quadratic interpolation of sub-pixel location
01268                     PixelType del = (m1 - m3) / 2.0;
01269                     del /= (m1 + m3 - 2.0*m2);
01270                     edgel.x = x - del;
01271                     edgel.y = y + del;
01272                     edgel.strength = m2;
01273                     
01274                     maximum_found = true;                    
01275                 }
01276             }
01277             else
01278             {
01279                 // north-east-south-west edge
01280                 PixelType m1 = magnitude(x-1, y-1);
01281                 PixelType m2 = magnitude(x, y);
01282                 PixelType m3 = magnitude(x+1, y+1);
01283                 
01284                 if(m1 < m2 && m3 <= m2)
01285                 {
01286                     // local maximum => quadratic interpolation of sub-pixel location
01287                     PixelType del = (m1 - m3) / 2.0;
01288                     del /= (m1 + m3 - 2.0*m2);
01289                     edgel.x = x + del;
01290                     edgel.y = y + del;
01291                     edgel.strength = m2;
01292                     
01293                     maximum_found = true;                    
01294                 }
01295             }
01296             
01297             if(maximum_found)
01298             {
01299                 double orientation = VIGRA_CSTD::atan2(-grady, gradx) - M_PI * 1.5;
01300                 if(orientation < 0.0)
01301                     orientation += 2.0*M_PI;
01302                 edgel.orientation = orientation;
01303                 edgels.push_back(edgel);
01304             }
01305         }
01306     }
01307 }
01308 
01309 /********************************************************/
01310 /*                                                      */
01311 /*                      cannyEdgelList                  */
01312 /*                                                      */
01313 /********************************************************/
01314 
01315 /** \brief Simple implementation of Canny's edge detector.
01316 
01317     This operator first calculates the gradient vector for each
01318     pixel of the image using first derivatives of a Gaussian at the 
01319     given scale. Then a very simple non-maxima supression is performed: 
01320     for each 3x3 neighborhood, it is determined whether the center pixel has 
01321     larger gradient magnitude than its two neighbors in gradient direction
01322     (where the direction is rounded into octands). If this is the case,
01323     a new \ref Edgel is appended to the given vector of <TT>edgels</TT>. The subpixel
01324     edgel position is determined by fitting a parabola 
01325     to the three gradient magnitude values 
01326     mentioned above. The sub-pixel location of the parabola's tip 
01327     and the gradient magnitude and direction are written in the newly created edgel.
01328     
01329     <b> Declarations:</b>
01330     
01331     pass arguments explicitly:
01332     \code
01333     namespace vigra {
01334         template <class SrcIterator, class SrcAccessor>
01335         void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01336                                 std::vector<Edgel> & edgels, double scale);
01337     }
01338     \endcode
01339     
01340     use argument objects in conjunction with \ref ArgumentObjectFactories:
01341     \code
01342     namespace vigra {
01343         template <class SrcIterator, class SrcAccessor>
01344         void 
01345         cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01346                        std::vector<Edgel> & edgels, double scale);
01347     }
01348     \endcode
01349     
01350     <b> Usage:</b>
01351     
01352     <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01353     Namespace: vigra
01354     
01355     \code
01356     vigra::BImage src(w,h);
01357     
01358     // empty edgel list
01359     std::vector<vigra::Edgel> edgels;
01360     ...
01361     
01362     // find edgels at scale 0.8  
01363     vigra::cannyEdgelList(srcImageRange(src), edgels, 0.8);
01364     \endcode
01365 
01366     <b> Required Interface:</b>
01367     
01368     \code
01369     SrcImageIterator src_upperleft;
01370     SrcAccessor src_accessor;
01371     
01372     src_accessor(src_upperleft);
01373     \endcode
01374     
01375     SrcAccessor::value_type must be a type convertible to float
01376     
01377     <b> Preconditions:</b>
01378     
01379     \code
01380     scale > 0
01381     \endcode
01382 */
01383 template <class SrcIterator, class SrcAccessor>
01384 void cannyEdgelList(SrcIterator ul, SrcIterator lr, SrcAccessor src,
01385                         std::vector<Edgel> & edgels, double scale)
01386 {
01387     int w = lr.x - ul.x;
01388     int h = lr.y - ul.y;
01389     
01390     // calculate image gradients
01391     typedef typename 
01392         NumericTraits<typename SrcAccessor::value_type>::RealPromote
01393         TmpType;
01394         
01395     BasicImage<TmpType> tmp(w,h), dx(w,h), dy(w,h);
01396     
01397     Kernel1D<double> smooth, grad;
01398     
01399     smooth.initGaussian(scale);
01400     grad.initGaussianDerivative(scale, 1);
01401     
01402     separableConvolveX(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
01403     separableConvolveY(srcImageRange(tmp), destImage(dx), kernel1d(smooth));
01404     
01405     separableConvolveY(srcIterRange(ul, lr, src), destImage(tmp), kernel1d(grad));
01406     separableConvolveX(srcImageRange(tmp), destImage(dy), kernel1d(smooth));
01407     
01408     combineTwoImages(srcImageRange(dx), srcImage(dy), destImage(tmp),
01409                      MagnitudeFunctor<TmpType>());
01410     
01411     
01412     // find edgels
01413     internalCannyFindEdgels(dx, dy, tmp, edgels);
01414 }
01415 
01416 template <class SrcIterator, class SrcAccessor>
01417 inline void 
01418 cannyEdgelList(triple<SrcIterator, SrcIterator, SrcAccessor> src,
01419                std::vector<Edgel> & edgels, double scale)
01420 {
01421     cannyEdgelList(src.first, src.second, src.third, edgels, scale);
01422 }
01423 
01424 /********************************************************/
01425 /*                                                      */
01426 /*                       cannyEdgeImage                 */
01427 /*                                                      */
01428 /********************************************************/
01429 
01430 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
01431 
01432     This operator first calls \ref cannyEdgelList() to generate an 
01433     edgel list for the given image. Then it scans this list and selects edgels
01434     whose strength is above the given <TT>gradient_threshold</TT>. For each of these 
01435     edgels, the edgel's location is rounded to the nearest pixel, and that
01436     pixel marked with the given <TT>edge_marker</TT>.
01437     
01438     <b> Declarations:</b>
01439     
01440     pass arguments explicitly:
01441     \code
01442     namespace vigra {
01443         template <class SrcIterator, class SrcAccessor,
01444                   class DestIterator, class DestAccessor, 
01445                   class GradValue, class DestValue>
01446         void cannyEdgeImage(
01447                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01448                    DestIterator dul, DestAccessor da,
01449                    double scale, GradValue gradient_threshold, DestValue edge_marker);
01450     }
01451     \endcode
01452     
01453     use argument objects in conjunction with \ref ArgumentObjectFactories:
01454     \code
01455     namespace vigra {
01456         template <class SrcIterator, class SrcAccessor,
01457                   class DestIterator, class DestAccessor, 
01458                   class GradValue, class DestValue>
01459         inline void cannyEdgeImage(
01460                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01461                    pair<DestIterator, DestAccessor> dest,
01462                    double scale, GradValue gradient_threshold, DestValue edge_marker);
01463     }
01464     \endcode
01465     
01466     <b> Usage:</b>
01467     
01468     <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01469     Namespace: vigra
01470     
01471     \code
01472     vigra::BImage src(w,h), edges(w,h);
01473     
01474     // empty edge image
01475     edges = 0;
01476     ...
01477     
01478     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1 
01479     vigra::cannyEdgeImage(srcImageRange(src), destImage(edges), 
01480                                      0.8, 4.0, 1);
01481     \endcode
01482 
01483     <b> Required Interface:</b>
01484     
01485     see also: \ref cannyEdgelList().
01486     
01487     \code
01488     DestImageIterator dest_upperleft;
01489     DestAccessor dest_accessor;
01490     DestValue edge_marker;
01491     
01492     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
01493     \endcode
01494     
01495     <b> Preconditions:</b>
01496     
01497     \code
01498     scale > 0
01499     gradient_threshold > 0
01500     \endcode
01501 */
01502 template <class SrcIterator, class SrcAccessor,
01503           class DestIterator, class DestAccessor, 
01504           class GradValue, class DestValue>
01505 void cannyEdgeImage(
01506            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01507            DestIterator dul, DestAccessor da,
01508            double scale, GradValue gradient_threshold, DestValue edge_marker)
01509 {
01510     std::vector<Edgel> edgels;
01511     
01512     cannyEdgelList(sul, slr, sa, edgels, scale);
01513     
01514     for(unsigned int i=0; i<edgels.size(); ++i)
01515     {
01516         if(gradient_threshold < edgels[i].strength)
01517         {
01518             Diff2D pix((int)(edgels[i].x + 0.5), (int)(edgels[i].y + 0.5));
01519             
01520             da.set(edge_marker, dul, pix);
01521         }
01522     }
01523 }
01524 
01525 template <class SrcIterator, class SrcAccessor,
01526           class DestIterator, class DestAccessor, 
01527           class GradValue, class DestValue>
01528 inline void cannyEdgeImage(
01529            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01530            pair<DestIterator, DestAccessor> dest,
01531            double scale, GradValue gradient_threshold, DestValue edge_marker)
01532 {
01533     cannyEdgeImage(src.first, src.second, src.third,
01534                    dest.first, dest.second,
01535                    scale, gradient_threshold, edge_marker);
01536 }
01537 
01538 /********************************************************/
01539 
01540 namespace detail {
01541 
01542 template <class DestIterator>
01543 int neighborhoodConfiguration(DestIterator dul)
01544 {
01545     int v = 0;
01546     NeighborhoodCirculator<DestIterator, EightNeighborCode> c(dul, EightNeighborCode::SouthEast);
01547     for(int i=0; i<8; ++i, --c)
01548     {
01549         v = (v << 1) | ((*c != 0) ? 1 : 0);
01550     }
01551     
01552     return v;
01553 }
01554 
01555 template <class GradValue>
01556 struct SimplePoint
01557 {
01558     Diff2D point;
01559     GradValue grad;
01560     
01561     SimplePoint(Diff2D const & p, GradValue g)
01562     : point(p), grad(g)
01563     {}
01564     
01565     bool operator<(SimplePoint const & o) const
01566     {
01567         return grad < o.grad; 
01568     }
01569 };
01570 
01571 } // namespace detail
01572 
01573 /********************************************************/
01574 /*                                                      */
01575 /*              cannyEdgeImageWithThinning              */
01576 /*                                                      */
01577 /********************************************************/
01578 
01579 /** \brief Detect and mark edges in an edge image using Canny's algorithm.
01580 
01581     This operator first calls \ref cannyEdgeImage() to generate an 
01582     edge image. The resulting edge pixels are then subjected to topological thinning
01583     so that the remaining edge pixels can be linked into edgel chains with a provable,
01584     non-heuristic algorithm. Optionally, the outermost pixels are marked as edge pixels
01585     as well when <tt>addBorder</tt> is true.
01586     
01587     <b> Declarations:</b>
01588     
01589     pass arguments explicitly:
01590     \code
01591     namespace vigra {
01592         template <class SrcIterator, class SrcAccessor,
01593                   class DestIterator, class DestAccessor, 
01594                   class GradValue, class DestValue>
01595         void cannyEdgeImageWithThinning(
01596                    SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01597                    DestIterator dul, DestAccessor da,
01598                    double scale, GradValue gradient_threshold, 
01599                    DestValue edge_marker, bool addBorder = true);
01600     }
01601     \endcode
01602     
01603     use argument objects in conjunction with \ref ArgumentObjectFactories:
01604     \code
01605     namespace vigra {
01606         template <class SrcIterator, class SrcAccessor,
01607                   class DestIterator, class DestAccessor, 
01608                   class GradValue, class DestValue>
01609         void cannyEdgeImageWithThinning(
01610                    triple<SrcIterator, SrcIterator, SrcAccessor> src,
01611                    pair<DestIterator, DestAccessor> dest,
01612                    double scale, GradValue gradient_threshold, 
01613                    DestValue edge_marker, bool addBorder = true);
01614     }
01615     \endcode
01616     
01617     <b> Usage:</b>
01618     
01619     <b>\#include</b> "<a href="edgedetection_8hxx-source.html">vigra/edgedetection.hxx</a>"<br>
01620     Namespace: vigra
01621     
01622     \code
01623     vigra::BImage src(w,h), edges(w,h);
01624     
01625     // empty edge image
01626     edges = 0;
01627     ...
01628     
01629     // find edges at scale 0.8 with gradient larger than 4.0, mark with 1, annd add border
01630     vigra::cannyEdgeImageWithThinning(srcImageRange(src), destImage(edges), 
01631                                      0.8, 4.0, 1, true);
01632     \endcode
01633 
01634     <b> Required Interface:</b>
01635     
01636     see also: \ref cannyEdgelList().
01637     
01638     \code
01639     DestImageIterator dest_upperleft;
01640     DestAccessor dest_accessor;
01641     DestValue edge_marker;
01642     
01643     dest_accessor.set(edge_marker, dest_upperleft, vigra::Diff2D(1,1));
01644     \endcode
01645     
01646     <b> Preconditions:</b>
01647     
01648     \code
01649     scale > 0
01650     gradient_threshold > 0
01651     \endcode
01652 */
01653 template <class SrcIterator, class SrcAccessor,
01654           class DestIterator, class DestAccessor, 
01655           class GradValue, class DestValue>
01656 void cannyEdgeImageWithThinning(
01657            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01658            DestIterator dul, DestAccessor da,
01659            double scale, GradValue gradient_threshold, 
01660            DestValue edge_marker, bool addBorder)
01661 {
01662     int w = slr.x - sul.x;
01663     int h = slr.y - sul.y;
01664     
01665     BImage edgeImage(w, h, BImage::value_type(0));
01666     BImage::traverser eul = edgeImage.upperLeft();
01667     BImage::Accessor ea = edgeImage.accessor();
01668     if(addBorder)
01669         initImageBorder(destImageRange(edgeImage), 1, 1);
01670     cannyEdgeImage(sul, slr, sa, eul, ea, 
01671                    scale, gradient_threshold, 1);
01672     
01673     static bool isSimplePoint[256] = {
01674         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 
01675         0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 
01676         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 
01677         1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 
01678         0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 
01679         0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 
01680         0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 
01681         1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 
01682         0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 
01683         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
01684         0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 
01685         0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 
01686         1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 
01687         0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 
01688         1, 0, 1, 0 };
01689         
01690     eul += Diff2D(1,1);
01691     sul += Diff2D(1,1);
01692     int w2 = w-2;
01693     int h2 = h-2;
01694     
01695     typedef detail::SimplePoint<GradValue> SP;
01696     std::priority_queue<SP, std::vector<SP> >  pqueue;
01697     
01698     Diff2D p(0,0);
01699     for(; p.y < h2; ++p.y)
01700     {
01701         for(p.x = 0; p.x < w2; ++p.x)
01702         {
01703             BImage::traverser e = eul + p;
01704             if(*e == 0)
01705                 continue;
01706             int v = detail::neighborhoodConfiguration(e);
01707             if(isSimplePoint[v])
01708             {
01709                 pqueue.push(SP(p, norm(sa(sul+p))));
01710                 *e = 2; // remember that it is already in queue
01711             }
01712         }
01713     }
01714     
01715     static const Diff2D dist[] = { Diff2D(-1,0), Diff2D(0,-1),
01716                                    Diff2D(1,0),  Diff2D(0,1) };
01717 
01718     while(pqueue.size())
01719     {
01720         p = pqueue.top().point;
01721         pqueue.pop();
01722         
01723         BImage::traverser e = eul + p;
01724         int v = detail::neighborhoodConfiguration(e);
01725         if(!isSimplePoint[v])
01726             continue; // point may no longer be simple because its neighbors changed
01727             
01728         *e = 0; // delete simple point
01729         
01730         for(int i=0; i<4; ++i)
01731         {
01732             Diff2D pneu = p + dist[i];
01733             if(pneu.x == -1 || pneu.y == -1 || pneu.x == w2 || pneu.y == h2)
01734                 continue; // do not remove points at the border
01735                 
01736             BImage::traverser eneu = eul + pneu;
01737             if(*eneu == 1) // point is boundary and not yet in the queue
01738             {
01739                 int v = detail::neighborhoodConfiguration(eneu);
01740                 if(isSimplePoint[v])
01741                 {
01742                     pqueue.push(SP(pneu, norm(sa(sul+pneu))));
01743                     *eneu = 2; // remember that it is already in queue
01744                 }
01745             }
01746         }
01747     }
01748     
01749     initImageIf(destIterRange(dul, dul+Diff2D(w,h), da),
01750                 maskImage(edgeImage), edge_marker);
01751 }
01752 
01753 template <class SrcIterator, class SrcAccessor,
01754           class DestIterator, class DestAccessor, 
01755           class GradValue, class DestValue>
01756 inline void cannyEdgeImageWithThinning(
01757            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01758            pair<DestIterator, DestAccessor> dest,
01759            double scale, GradValue gradient_threshold, 
01760            DestValue edge_marker, bool addBorder)
01761 {
01762     cannyEdgeImageWithThinning(src.first, src.second, src.third,
01763                                dest.first, dest.second,
01764                                scale, gradient_threshold, edge_marker, addBorder);
01765 }
01766 
01767 template <class SrcIterator, class SrcAccessor,
01768           class DestIterator, class DestAccessor, 
01769           class GradValue, class DestValue>
01770 inline void cannyEdgeImageWithThinning(
01771            SrcIterator sul, SrcIterator slr, SrcAccessor sa,
01772            DestIterator dul, DestAccessor da,
01773            double scale, GradValue gradient_threshold, DestValue edge_marker)
01774 {
01775     cannyEdgeImageWithThinning(sul, slr, sa,
01776                                dul, da,
01777                                scale, gradient_threshold, edge_marker, true);
01778 }
01779 
01780 template <class SrcIterator, class SrcAccessor,
01781           class DestIterator, class DestAccessor, 
01782           class GradValue, class DestValue>
01783 inline void cannyEdgeImageWithThinning(
01784            triple<SrcIterator, SrcIterator, SrcAccessor> src,
01785            pair<DestIterator, DestAccessor> dest,
01786            double scale, GradValue gradient_threshold, DestValue edge_marker)
01787 {
01788     cannyEdgeImageWithThinning(src.first, src.second, src.third,
01789                                dest.first, dest.second,
01790                                scale, gradient_threshold, edge_marker, true);
01791 }
01792 
01793 //@}
01794 
01795 /** \page CrackEdgeImage Crack Edge Image
01796 
01797 Crack edges are marked <i>between</i> the pixels of an image. 
01798 A Crack Edge Image is an image that represents these edges. In order
01799 to accomodate the cracks, the Crack Edge Image must be twice as large
01800 as the original image (precisely (2*w - 1) by (2*h - 1)). A Crack Edge Image
01801 can easily be derived from a binary image or from the signs of the 
01802 response of a Laplacean filter. Consider the following sketch, where
01803 <TT>+</TT> encodes the foreground, <TT>-</TT> the background, and
01804 <TT>*</TT> the resulting crack edges.
01805 
01806     \code
01807 sign of difference image         insert cracks         resulting CrackEdgeImage
01808 
01809                                    + . - . -              . * . . .
01810       + - -                        . . . . .              . * * * .
01811       + + -               =>       + . + . -      =>      . . . * .
01812       + + +                        . . . . .              . . . * *
01813                                    + . + . +              . . . . .
01814     \endcode
01815 
01816 Starting from the original binary image (left), we insert crack pixels
01817 to get to the double-sized image (center). Finally, we mark all 
01818 crack pixels whose non-crack neighbors have different signs as 
01819 crack edge points, while all other pixels (crack and non-crack) become 
01820 region pixels.
01821 
01822 <b>Requirements on a Crack Edge Image:</b>
01823 
01824 <ul>
01825     <li>Crack Edge Images have odd width and height.
01826     <li>Crack pixels have at least one odd coordinate.
01827     <li>Only crack pixels may be marked as edge points.
01828     <li>Crack pixels with two odd coordinates must be marked as edge points
01829         whenever any of their neighboring crack pixels was marked.  
01830 </ul>
01831 
01832 The last two requirements ensure that both edges and regions are 4-connected. 
01833 Thus, 4-connectivity and 8-connectivity yield identical connected 
01834 components in a Crack Edge Image (so called <i>well-composedness</i>).
01835 This ensures that Crack Edge Images have nice topological properties
01836 (cf. L. J. Latecki: "Well-Composed Sets", Academic Press, 2000). 
01837 */
01838 
01839 
01840 } // namespace vigra
01841 
01842 #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.4.0 (21 Dec 2005)