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

vigra/distancetransform.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2002 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00012 /*        vigra@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_DISTANCETRANSFORM_HXX
00040 #define VIGRA_DISTANCETRANSFORM_HXX
00041 
00042 #include <cmath>
00043 #include "stdimage.hxx"
00044 
00045 namespace vigra {
00046 
00047 /*
00048  * functors to determine the distance norm 
00049  * these functors assume that dx and dy are positive
00050  * (this is OK for use in internalDistanceTransform())
00051  */
00052  
00053 // chessboard metric
00054 struct InternalDistanceTransformLInifinityNormFunctor
00055 {
00056     float operator()(float dx, float dy) const
00057     {
00058         return (dx < dy) ? dy : dx;
00059     }
00060 };
00061 
00062 // Manhattan metric
00063 struct InternalDistanceTransformL1NormFunctor
00064 {
00065     float operator()(float dx, float dy) const
00066     {
00067         return dx + dy;
00068     }
00069 };
00070 
00071 // Euclidean metric
00072 struct InternalDistanceTransformL2NormFunctor
00073 {
00074     float operator()(float dx, float dy) const
00075     {
00076         return VIGRA_CSTD::sqrt(dx*dx + dy*dy);
00077     }
00078 };
00079 
00080 
00081 template <class SrcImageIterator, class SrcAccessor,
00082                    class DestImageIterator, class DestAccessor,
00083                    class ValueType, class NormFunctor>
00084 void
00085 internalDistanceTransform(SrcImageIterator src_upperleft, 
00086                 SrcImageIterator src_lowerright, SrcAccessor sa,
00087                 DestImageIterator dest_upperleft, DestAccessor da,
00088                 ValueType background, NormFunctor norm)
00089 {
00090     int w = src_lowerright.x - src_upperleft.x;  
00091     int h = src_lowerright.y - src_upperleft.y;  
00092     
00093     FImage xdist(w,h), ydist(w,h);
00094     
00095     xdist = w;    // init x and
00096     ydist = h;    // y distances with 'large' values
00097 
00098     SrcImageIterator sy = src_upperleft;
00099     DestImageIterator ry = dest_upperleft;
00100     FImage::Iterator xdy = xdist.upperLeft();
00101     FImage::Iterator ydy = ydist.upperLeft();
00102     SrcImageIterator sx = sy;
00103     DestImageIterator rx = ry;
00104     FImage::Iterator xdx = xdy;
00105     FImage::Iterator ydx = ydy;
00106     
00107     static const Diff2D left(-1, 0);
00108     static const Diff2D right(1, 0);
00109     static const Diff2D top(0, -1);
00110     static const Diff2D bottom(0, 1);
00111         
00112     int x,y;
00113     if(sa(sx) != background)    // first pixel
00114     {
00115         *xdx = 0.0;
00116         *ydx = 0.0;
00117         da.set(0.0, rx);
00118     }
00119     else
00120     {
00121         da.set(norm(*xdx, *ydx), rx);
00122     }
00123     
00124     
00125     for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x; 
00126         x<w; 
00127         ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x)   // first row left to right
00128     {
00129         if(sa(sx) != background)
00130         {
00131             *xdx = 0.0;
00132             *ydx = 0.0;
00133             da.set(0.0, rx);
00134         }
00135         else
00136         {
00137             *xdx  = xdx[left] + 1.0;   // propagate x and
00138             *ydx  = ydx[left];         // y components of distance from left pixel
00139             da.set(norm(*xdx, *ydx), rx); // calculate distance from x and y components
00140         }
00141     }
00142     for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2; 
00143         x>=0; 
00144         --x, --xdx.x, --ydx.x, --sx.x, --rx.x)   // first row right to left
00145     {
00146         float d = norm(xdx[right] + 1.0, ydx[right]);
00147         
00148         if(da(rx) < d) continue;
00149         
00150         *xdx = xdx[right] + 1.0;
00151         *ydx = ydx[right];
00152         da.set(d, rx);
00153     }
00154     for(y=1, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y; 
00155         y<h;
00156         ++y, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y)   // top to bottom
00157     {
00158         sx = sy;
00159         rx = ry;
00160         xdx = xdy;
00161         ydx = ydy;
00162         
00163         if(sa(sx) != background)    // first pixel of current row
00164         {
00165             *xdx = 0.0;
00166             *ydx = 0.0;
00167             da.set(0.0, rx);
00168         }
00169         else
00170         {
00171             *xdx = xdx[top];
00172             *ydx = ydx[top] + 1.0;
00173             da.set(norm(*xdx, *ydx), rx);
00174         }
00175         
00176         for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x; 
00177             x<w; 
00178             ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x)  // current row left to right
00179         {
00180             if(sa(sx) != background)
00181             {
00182                 *xdx = 0.0;
00183                 *ydx = 0.0;
00184                 da.set(0.0, rx);
00185             }
00186             else
00187             {
00188                 float d1 = norm(xdx[left] + 1.0, ydx[left]);
00189                 float d2 = norm(xdx[top], ydx[top] + 1.0);
00190                 
00191                 if(d1 < d2)
00192                 {
00193                     *xdx = xdx[left] + 1.0;
00194                     *ydx = ydx[left];
00195                     da.set(d1, rx);
00196                 }
00197                 else
00198                 {
00199                     *xdx = xdx[top];
00200                     *ydx = ydx[top] + 1.0;
00201                     da.set(d2, rx);
00202                 }
00203             }
00204         }
00205         for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2; 
00206             x>=0; 
00207             --x, --xdx.x, --ydx.x, --sx.x, --rx.x)  // current row right to left
00208         {
00209             float d1 = norm(xdx[right] + 1.0, ydx[right]);
00210             
00211             if(da(rx) < d1) continue;
00212             
00213             *xdx = xdx[right] + 1.0;
00214             *ydx = ydx[right];
00215             da.set(d1, rx);
00216         }
00217     }
00218     for(y=h-2, xdy.y -= 2, ydy.y -= 2, sy.y -= 2, ry.y -= 2; 
00219         y>=0;
00220         --y, --xdy.y, --ydy.y, --sy.y, --ry.y)    // bottom to top
00221     {
00222         sx = sy;
00223         rx = ry;
00224         xdx = xdy;
00225         ydx = ydy;
00226         
00227         float d = norm(xdx[bottom], ydx[bottom] + 1.0);
00228         if(d < da(rx))    // first pixel of current row
00229         { 
00230             *xdx = xdx[bottom];
00231             *ydx = ydx[bottom] + 1.0;
00232             da.set(d, rx);
00233         }
00234             
00235         for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x; 
00236             x<w;
00237             ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x)  // current row left to right
00238         {
00239             float d1 = norm(xdx[left] + 1.0, ydx[left]);
00240             float d2 = norm(xdx[bottom], ydx[bottom] + 1.0);
00241             
00242             if(d1 < d2)
00243             {
00244                 if(da(rx) < d1) continue;
00245                 *xdx = xdx[left] + 1.0;
00246                 *ydx = ydx[left];
00247                 da.set(d1, rx);
00248             }
00249             else
00250             {
00251                 if(da(rx) < d2) continue;
00252                 *xdx = xdx[bottom];
00253                 *ydx = ydx[bottom] + 1.0;
00254                 da.set(d2, rx);
00255             }
00256         }
00257         for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2; 
00258             x>=0; 
00259             --x, --xdx.x, --ydx.x, --sx.x, --rx.x)  // current row right to left
00260         {
00261             float d1 = norm(xdx[right] + 1.0, ydx[right]);
00262 
00263             if(da(rx) < d1) continue;
00264             *xdx = xdx[right] + 1.0;
00265             *ydx = ydx[right];
00266             da.set(d1, rx);
00267         }
00268     }
00269 }
00270 
00271 /********************************************************/
00272 /*                                                      */
00273 /*                 distanceTransform                    */
00274 /*                                                      */
00275 /********************************************************/
00276 
00277 /** \addtogroup DistanceTransform Distance Transform
00278     Perform a distance transform using either the Euclidean, Manhattan, 
00279     or chessboard metrics.
00280     
00281     See also: \ref MultiArrayDistanceTransform "multi-dimensional distance transforms"
00282 */
00283 //@{
00284 
00285 /** For all background pixels, calculate the distance to 
00286     the nearest object or contour. The label of the pixels to be considered 
00287     background in the source image is passed in the parameter 'background'.
00288     Source pixels with other labels will be considered objects. In the 
00289     destination image, all pixels corresponding to background will be assigned 
00290     the their distance value, all pixels corresponding to objects will be
00291     assigned 0.
00292     
00293     The parameter 'norm' gives the distance norm to be used:
00294     
00295     <ul>
00296 
00297     <li> norm == 0: use chessboard distance (L-infinity norm)
00298     <li> norm == 1: use Manhattan distance (L1 norm)
00299     <li> norm == 2: use Euclidean distance (L2 norm)
00300     
00301     </ul>
00302     
00303     If you use the L2 norm, the destination pixels must be real valued to give
00304     correct results.
00305     
00306     <b> Declarations:</b>
00307     
00308     pass arguments explicitly:
00309     \code
00310     namespace vigra {
00311         template <class SrcImageIterator, class SrcAccessor,
00312                            class DestImageIterator, class DestAccessor,
00313                            class ValueType>
00314         void distanceTransform(SrcImageIterator src_upperleft, 
00315                         SrcImageIterator src_lowerright, SrcAccessor sa,
00316                         DestImageIterator dest_upperleft, DestAccessor da,
00317                         ValueType background, int norm)
00318     }
00319     \endcode
00320     
00321     
00322     use argument objects in conjunction with \ref ArgumentObjectFactories :
00323     \code
00324     namespace vigra {
00325         template <class SrcImageIterator, class SrcAccessor,
00326                            class DestImageIterator, class DestAccessor,
00327                            class ValueType>
00328         void distanceTransform(
00329             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00330             pair<DestImageIterator, DestAccessor> dest,
00331             ValueType background, int norm)
00332     }
00333     \endcode
00334     
00335     <b> Usage:</b>
00336     
00337     <b>\#include</b> <<a href="distancetransform_8hxx-source.html">vigra/distancetransform.hxx</a>><br>
00338     Namespace: vigra
00339     
00340     
00341     \code
00342     
00343     vigra::BImage src(w,h), edges(w,h);
00344     vigra::FImage distance(w, h);
00345 
00346     // empty edge image
00347     edges = 0;
00348     ...
00349 
00350     // detect edges in src image (edges will be marked 1, background 0)
00351     vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 
00352                                      0.8, 4.0);
00353      
00354     // find distance of all pixels from nearest edge
00355     vigra::distanceTransform(srcImageRange(edges), destImage(distance),
00356                              0,                   2);
00357     //                       ^ background label   ^ norm (Euclidean)
00358     \endcode
00359 
00360     <b> Required Interface:</b>
00361     
00362     \code
00363     
00364     SrcImageIterator src_upperleft, src_lowerright;
00365     DestImageIterator dest_upperleft;
00366     
00367     SrcAccessor sa;
00368     DestAccessor da;
00369     
00370     ValueType background;
00371     float distance;
00372     
00373     sa(src_upperleft) != background;
00374     da(dest_upperleft) < distance;
00375     da.set(distance, dest_upperleft);
00376  
00377     \endcode
00378 */
00379 doxygen_overloaded_function(template <...> void distanceTransform)
00380 
00381 template <class SrcImageIterator, class SrcAccessor,
00382                    class DestImageIterator, class DestAccessor,
00383                    class ValueType>
00384 inline void
00385 distanceTransform(SrcImageIterator src_upperleft, 
00386                 SrcImageIterator src_lowerright, SrcAccessor sa,
00387                 DestImageIterator dest_upperleft, DestAccessor da,
00388                 ValueType background, int norm)
00389 {
00390     if(norm == 1)
00391     {
00392         internalDistanceTransform(src_upperleft, src_lowerright, sa,
00393                                   dest_upperleft, da, background,
00394                                   InternalDistanceTransformL1NormFunctor());
00395     }
00396     else if(norm == 2)
00397     {
00398         internalDistanceTransform(src_upperleft, src_lowerright, sa,
00399                                   dest_upperleft, da, background,
00400                                   InternalDistanceTransformL2NormFunctor());
00401     }
00402     else
00403     {
00404         internalDistanceTransform(src_upperleft, src_lowerright, sa,
00405                                   dest_upperleft, da, background,
00406                                   InternalDistanceTransformLInifinityNormFunctor());
00407     }
00408 }
00409 
00410 template <class SrcImageIterator, class SrcAccessor,
00411                    class DestImageIterator, class DestAccessor,
00412                    class ValueType>
00413 inline void
00414 distanceTransform(
00415     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00416     pair<DestImageIterator, DestAccessor> dest,
00417     ValueType background, int norm)
00418 {
00419     distanceTransform(src.first, src.second, src.third,
00420                       dest.first, dest.second, background, norm);
00421 }
00422 
00423 //@}
00424 
00425 } // namespace vigra
00426 
00427 #endif // VIGRA_DISTANCETRANSFORM_HXX
00428 

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

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)