[ 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 /*    The VIGRA Website is                                              */
00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */                
00034 /*                                                                      */
00035 /************************************************************************/
00036  
00037  
00038 #ifndef VIGRA_DISTANCETRANSFORM_HXX
00039 #define VIGRA_DISTANCETRANSFORM_HXX
00040 
00041 #include <cmath>
00042 #include "stdimage.hxx"
00043 
00044 namespace vigra {
00045 
00046 /*
00047  * functors to determine the distance norm 
00048  * these functors assume that dx and dy are positive
00049  * (this is OK for use in internalDistanceTransform())
00050  */
00051  
00052 // chessboard metric
00053 struct InternalDistanceTransformLInifinityNormFunctor
00054 {
00055     float operator()(float dx, float dy) const
00056     {
00057         return (dx < dy) ? dy : dx;
00058     }
00059 };
00060 
00061 // Manhattan metric
00062 struct InternalDistanceTransformL1NormFunctor
00063 {
00064     float operator()(float dx, float dy) const
00065     {
00066         return dx + dy;
00067     }
00068 };
00069 
00070 // Euclidean metric
00071 struct InternalDistanceTransformL2NormFunctor
00072 {
00073     float operator()(float dx, float dy) const
00074     {
00075         return VIGRA_CSTD::sqrt(dx*dx + dy*dy);
00076     }
00077 };
00078 
00079 
00080 template <class SrcImageIterator, class SrcAccessor,
00081                    class DestImageIterator, class DestAccessor,
00082                    class ValueType, class NormFunctor>
00083 void
00084 internalDistanceTransform(SrcImageIterator src_upperleft, 
00085                 SrcImageIterator src_lowerright, SrcAccessor sa,
00086                 DestImageIterator dest_upperleft, DestAccessor da,
00087                 ValueType background, NormFunctor norm)
00088 {
00089     int w = src_lowerright.x - src_upperleft.x;  
00090     int h = src_lowerright.y - src_upperleft.y;  
00091     
00092     FImage xdist(w,h), ydist(w,h);
00093     
00094     xdist = w;    // init x and
00095     ydist = h;    // y distances with 'large' values
00096 
00097     SrcImageIterator sy = src_upperleft;
00098     DestImageIterator ry = dest_upperleft;
00099     FImage::Iterator xdy = xdist.upperLeft();
00100     FImage::Iterator ydy = ydist.upperLeft();
00101     SrcImageIterator sx = sy;
00102     DestImageIterator rx = ry;
00103     FImage::Iterator xdx = xdy;
00104     FImage::Iterator ydx = ydy;
00105     
00106     static const Diff2D left(-1, 0);
00107     static const Diff2D right(1, 0);
00108     static const Diff2D top(0, -1);
00109     static const Diff2D bottom(0, 1);
00110         
00111     int x,y;
00112     if(sa(sx) != background)    // first pixel
00113     {
00114         *xdx = 0.0;
00115         *ydx = 0.0;
00116         da.set(0.0, rx);
00117     }
00118     else
00119     {
00120         da.set(norm(*xdx, *ydx), rx);
00121     }
00122     
00123     
00124     for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x; 
00125         x<w; 
00126         ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x)   // first row left to right
00127     {
00128         if(sa(sx) != background)
00129         {
00130             *xdx = 0.0;
00131             *ydx = 0.0;
00132             da.set(0.0, rx);
00133         }
00134         else
00135         {
00136             *xdx  = xdx[left] + 1.0;   // propagate x and
00137             *ydx  = ydx[left];         // y components of distance from left pixel
00138             da.set(norm(*xdx, *ydx), rx); // calculate distance from x and y components
00139         }
00140     }
00141     for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2; 
00142         x>=0; 
00143         --x, --xdx.x, --ydx.x, --sx.x, --rx.x)   // first row right to left
00144     {
00145         float d = norm(xdx[right] + 1.0, ydx[right]);
00146         
00147         if(da(rx) < d) continue;
00148         
00149         *xdx = xdx[right] + 1.0;
00150         *ydx = ydx[right];
00151         da.set(d, rx);
00152     }
00153     for(y=1, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y; 
00154         y<h;
00155         ++y, ++xdy.y, ++ydy.y, ++sy.y, ++ry.y)   // top to bottom
00156     {
00157         sx = sy;
00158         rx = ry;
00159         xdx = xdy;
00160         ydx = ydy;
00161         
00162         if(sa(sx) != background)    // first pixel of current row
00163         {
00164             *xdx = 0.0;
00165             *ydx = 0.0;
00166             da.set(0.0, rx);
00167         }
00168         else
00169         {
00170             *xdx = xdx[top];
00171             *ydx = ydx[top] + 1.0;
00172             da.set(norm(*xdx, *ydx), rx);
00173         }
00174         
00175         for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x; 
00176             x<w; 
00177             ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x)  // current row left to right
00178         {
00179             if(sa(sx) != background)
00180             {
00181                 *xdx = 0.0;
00182                 *ydx = 0.0;
00183                 da.set(0.0, rx);
00184             }
00185             else
00186             {
00187                 float d1 = norm(xdx[left] + 1.0, ydx[left]);
00188                 float d2 = norm(xdx[top], ydx[top] + 1.0);
00189                 
00190                 if(d1 < d2)
00191                 {
00192                     *xdx = xdx[left] + 1.0;
00193                     *ydx = ydx[left];
00194                     da.set(d1, rx);
00195                 }
00196                 else
00197                 {
00198                     *xdx = xdx[top];
00199                     *ydx = ydx[top] + 1.0;
00200                     da.set(d2, rx);
00201                 }
00202             }
00203         }
00204         for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2; 
00205             x>=0; 
00206             --x, --xdx.x, --ydx.x, --sx.x, --rx.x)  // current row right to left
00207         {
00208             float d1 = norm(xdx[right] + 1.0, ydx[right]);
00209             
00210             if(da(rx) < d1) continue;
00211             
00212             *xdx = xdx[right] + 1.0;
00213             *ydx = ydx[right];
00214             da.set(d1, rx);
00215         }
00216     }
00217     for(y=h-2, xdy.y -= 2, ydy.y -= 2, sy.y -= 2, ry.y -= 2; 
00218         y>=0;
00219         --y, --xdy.y, --ydy.y, --sy.y, --ry.y)    // bottom to top
00220     {
00221         sx = sy;
00222         rx = ry;
00223         xdx = xdy;
00224         ydx = ydy;
00225         
00226         float d = norm(xdx[bottom], ydx[bottom] + 1.0);
00227         if(d < da(rx))    // first pixel of current row
00228         { 
00229             *xdx = xdx[bottom];
00230             *ydx = ydx[bottom] + 1.0;
00231             da.set(d, rx);
00232         }
00233             
00234         for(x=1, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x; 
00235             x<w;
00236             ++x, ++xdx.x, ++ydx.x, ++sx.x, ++rx.x)  // current row left to right
00237         {
00238             float d1 = norm(xdx[left] + 1.0, ydx[left]);
00239             float d2 = norm(xdx[bottom], ydx[bottom] + 1.0);
00240             
00241             if(d1 < d2)
00242             {
00243                 if(da(rx) < d1) continue;
00244                 *xdx = xdx[left] + 1.0;
00245                 *ydx = ydx[left];
00246                 da.set(d1, rx);
00247             }
00248             else
00249             {
00250                 if(da(rx) < d2) continue;
00251                 *xdx = xdx[bottom];
00252                 *ydx = ydx[bottom] + 1.0;
00253                 da.set(d2, rx);
00254             }
00255         }
00256         for(x=w-2, xdx.x -= 2, ydx.x -= 2, sx.x -= 2, rx.x -= 2; 
00257             x>=0; 
00258             --x, --xdx.x, --ydx.x, --sx.x, --rx.x)  // current row right to left
00259         {
00260             float d1 = norm(xdx[right] + 1.0, ydx[right]);
00261 
00262             if(da(rx) < d1) continue;
00263             *xdx = xdx[right] + 1.0;
00264             *ydx = ydx[right];
00265             da.set(d1, rx);
00266         }
00267     }
00268 }
00269 
00270 /********************************************************/
00271 /*                                                      */
00272 /*                 distanceTransform                    */
00273 /*                                                      */
00274 /********************************************************/
00275 
00276 /** \addtogroup DistanceTransform Distance Transform
00277     Perform a distance transform using either the Euclidean, Manhattan, 
00278     or chessboard metrics.
00279     
00280     See also: \ref MultiArrayDistanceTransform "multi-dimensional distance transforms"
00281 */
00282 //@{
00283 
00284 /** For all background pixels, calculate the distance to 
00285     the nearest object or contour. The label of the pixels to be considered 
00286     background in the source image is passed in the parameter 'background'.
00287     Source pixels with other labels will be considered objects. In the 
00288     destination image, all pixels corresponding to background will be assigned 
00289     the their distance value, all pixels corresponding to objects will be
00290     assigned 0.
00291     
00292     The parameter 'norm' gives the distance norm to be used:
00293     
00294     <ul>
00295 
00296     <li> norm == 0: use chessboard distance (L-infinity norm)
00297     <li> norm == 1: use Manhattan distance (L1 norm)
00298     <li> norm == 2: use Euclidean distance (L2 norm)
00299     
00300     </ul>
00301     
00302     If you use the L2 norm, the destination pixels must be real valued to give
00303     correct results.
00304     
00305     <b> Declarations:</b>
00306     
00307     pass arguments explicitly:
00308     \code
00309     namespace vigra {
00310         template <class SrcImageIterator, class SrcAccessor,
00311                            class DestImageIterator, class DestAccessor,
00312                            class ValueType>
00313         void distanceTransform(SrcImageIterator src_upperleft, 
00314                         SrcImageIterator src_lowerright, SrcAccessor sa,
00315                         DestImageIterator dest_upperleft, DestAccessor da,
00316                         ValueType background, int norm)
00317     }
00318     \endcode
00319     
00320     
00321     use argument objects in conjunction with \ref ArgumentObjectFactories :
00322     \code
00323     namespace vigra {
00324         template <class SrcImageIterator, class SrcAccessor,
00325                            class DestImageIterator, class DestAccessor,
00326                            class ValueType>
00327         void distanceTransform(
00328             triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00329             pair<DestImageIterator, DestAccessor> dest,
00330             ValueType background, int norm)
00331     }
00332     \endcode
00333     
00334     <b> Usage:</b>
00335     
00336     <b>\#include</b> <<a href="distancetransform_8hxx-source.html">vigra/distancetransform.hxx</a>><br>
00337     Namespace: vigra
00338     
00339     
00340     \code
00341     
00342     vigra::BImage src(w,h), edges(w,h);
00343     vigra::FImage distance(w, h);
00344 
00345     // empty edge image
00346     edges = 0;
00347     ...
00348 
00349     // detect edges in src image (edges will be marked 1, background 0)
00350     vigra::differenceOfExponentialEdgeImage(srcImageRange(src), destImage(edges), 
00351                                      0.8, 4.0);
00352      
00353     // find distance of all pixels from nearest edge
00354     vigra::distanceTransform(srcImageRange(edges), destImage(distance),
00355                              0,                   2);
00356     //                       ^ background label   ^ norm (Euclidean)
00357     \endcode
00358 
00359     <b> Required Interface:</b>
00360     
00361     \code
00362     
00363     SrcImageIterator src_upperleft, src_lowerright;
00364     DestImageIterator dest_upperleft;
00365     
00366     SrcAccessor sa;
00367     DestAccessor da;
00368     
00369     ValueType background;
00370     float distance;
00371     
00372     sa(src_upperleft) != background;
00373     da(dest_upperleft) < distance;
00374     da.set(distance, dest_upperleft);
00375  
00376     \endcode
00377 */
00378 doxygen_overloaded_function(template <...> void distanceTransform)
00379 
00380 template <class SrcImageIterator, class SrcAccessor,
00381                    class DestImageIterator, class DestAccessor,
00382                    class ValueType>
00383 inline void
00384 distanceTransform(SrcImageIterator src_upperleft, 
00385                 SrcImageIterator src_lowerright, SrcAccessor sa,
00386                 DestImageIterator dest_upperleft, DestAccessor da,
00387                 ValueType background, int norm)
00388 {
00389     if(norm == 1)
00390     {
00391         internalDistanceTransform(src_upperleft, src_lowerright, sa,
00392                                   dest_upperleft, da, background,
00393                                   InternalDistanceTransformL1NormFunctor());
00394     }
00395     else if(norm == 2)
00396     {
00397         internalDistanceTransform(src_upperleft, src_lowerright, sa,
00398                                   dest_upperleft, da, background,
00399                                   InternalDistanceTransformL2NormFunctor());
00400     }
00401     else
00402     {
00403         internalDistanceTransform(src_upperleft, src_lowerright, sa,
00404                                   dest_upperleft, da, background,
00405                                   InternalDistanceTransformLInifinityNormFunctor());
00406     }
00407 }
00408 
00409 template <class SrcImageIterator, class SrcAccessor,
00410                    class DestImageIterator, class DestAccessor,
00411                    class ValueType>
00412 inline void
00413 distanceTransform(
00414     triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
00415     pair<DestImageIterator, DestAccessor> dest,
00416     ValueType background, int norm)
00417 {
00418     distanceTransform(src.first, src.second, src.third,
00419                       dest.first, dest.second, background, norm);
00420 }
00421 
00422 //@}
00423 
00424 } // namespace vigra
00425 
00426 #endif // VIGRA_DISTANCETRANSFORM_HXX
00427 

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

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