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

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

© 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.7.0 (Thu Aug 25 2011)