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

vigra/multi_distance.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*     Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn       */
00004 /*                        and Ullrich Koethe                            */
00005 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00006 /*                                                                      */
00007 /*    This file is part of the VIGRA computer vision library.           */
00008 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
00009 /*    The VIGRA Website is                                              */
00010 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00011 /*    Please direct questions, bug reports, and contributions to        */
00012 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00013 /*        vigra@informatik.uni-hamburg.de                               */
00014 /*                                                                      */
00015 /*    Permission is hereby granted, free of charge, to any person       */
00016 /*    obtaining a copy of this software and associated documentation    */
00017 /*    files (the "Software"), to deal in the Software without           */
00018 /*    restriction, including without limitation the rights to use,      */
00019 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00020 /*    sell copies of the Software, and to permit persons to whom the    */
00021 /*    Software is furnished to do so, subject to the following          */
00022 /*    conditions:                                                       */
00023 /*                                                                      */
00024 /*    The above copyright notice and this permission notice shall be    */
00025 /*    included in all copies or substantial portions of the             */
00026 /*    Software.                                                         */
00027 /*                                                                      */
00028 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00029 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00030 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00031 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00032 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00033 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00034 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00035 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00036 /*                                                                      */
00037 /************************************************************************/
00038 
00039 #ifndef VIGRA_MULTI_DISTANCE_HXX
00040 #define VIGRA_MULTI_DISTANCE_HXX
00041 
00042 #include <vector>
00043 #include <functional>
00044 #include "array_vector.hxx"
00045 #include "multi_array.hxx"
00046 #include "accessor.hxx"
00047 #include "numerictraits.hxx"
00048 #include "navigator.hxx"
00049 #include "metaprogramming.hxx"
00050 #include "multi_pointoperators.hxx"
00051 #include "functorexpression.hxx"
00052 
00053 namespace vigra
00054 {
00055 
00056 
00057 namespace detail
00058 {
00059 
00060 /********************************************************/
00061 /*                                                      */
00062 /*                distParabola                          */
00063 /*                                                      */
00064 /*  Version with sigma (parabola spread) for morphology */
00065 /*                                                      */
00066 /********************************************************/
00067 
00068 template <class SrcIterator, class SrcAccessor,
00069           class DestIterator, class DestAccessor >
00070 void distParabola(SrcIterator is, SrcIterator iend, SrcAccessor sa,
00071                               DestIterator id, DestAccessor da, float sigma )
00072 {
00073     // We assume that the data in the input is distance squared and treat it as such
00074     int w = iend - is;
00075     
00076     typedef typename NumericTraits<typename DestAccessor::value_type>::ValueType ValueType;
00077     typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote SumType;
00078     
00079     // Define the stack we use to determine the nearest background row 
00080     // (from previous dimension), the items on the stack will separate this column into
00081     // separate regions of influence. Each region of influence is closest to the same 
00082     // background row from the previous dimension.
00083     typedef triple<int, ValueType, int> influence;
00084     std::vector<influence> _stack;
00085 
00086     SrcIterator ibegin = is;
00087     _stack.push_back(influence(0, sa(is), w));
00088     
00089     ++is;
00090     int current = 1;
00091     
00092     int y0, y1, y2, y_dash, delta_y;
00093     sigma = sigma * sigma;
00094     bool nosigma = closeAtTolerance( sigma, 1.0 );
00095     
00096     y0 = 0;   // The beginning of the influence of row y1
00097    
00098     while( is != iend && current < w )
00099     {
00100         y1 = _stack.back().first;
00101         y2 = current;
00102         delta_y = y2 - y1;
00103 
00104         // If sigma is 1 (common case) avoid float multiplication here.
00105         if(nosigma)
00106             y_dash = (int)(sa(is) - _stack.back().second) - delta_y*delta_y;
00107         else
00108             y_dash = (int)(sigma * (sa(is) - _stack.back().second)) - delta_y*delta_y;
00109         y_dash = y_dash / (delta_y + delta_y);
00110         y_dash += y2;
00111 
00112         if( y_dash > y0)      
00113         {
00114             if( y_dash <= w )   // CASE 2 -- A new region of influence
00115             {
00116                 y0 = y_dash;
00117                 
00118                 _stack.back().third = y_dash; 
00119                 
00120                 _stack.push_back(influence(current, sa(is), w));
00121             }
00122 
00123             // CASE 1 -- This parabola is never active
00124             ++is;
00125             ++current;
00126             continue;
00127         } 
00128         else    // CASE 3 -- Parabola shadows the previous one completely
00129         {
00130             _stack.pop_back();
00131 
00132             if(_stack.size() < 2) 
00133                 y0=0;
00134             else    
00135                 y0=_stack[_stack.size()-2].third;
00136             
00137             if(_stack.empty())  // This row influences all previous rows.
00138             {
00139                 _stack.push_back(influence(current, sa(is), w));
00140 
00141                 ++is;
00142                 ++current;
00143                 continue;
00144             }
00145         }
00146     }
00147 
00148     // Now we have the stack indicating which rows are influenced by (and therefore
00149     // closest to) which row. We can go through the stack and calculate the
00150     // distance squared for each element of the column.
00151 
00152     typename std::vector<influence>::iterator it = _stack.begin();
00153 
00154     ValueType distance = 0;   // The distance squared
00155     current = 0;
00156     delta_y = 0;
00157     is = ibegin;
00158 
00159     for(; is != iend; ++current, ++id, ++is)
00160     {
00161         // FIXME FIXME Bound checking incorrect here? vvv
00162         if( current >= (*it).third && it != _stack.end()) ++it; 
00163        
00164         // FIXME FIXME The following check speeds things up for distance, but completely
00165         // messes up the grayscale morphology. Use an extra flag???
00166   /*      if( *is == 0 ) // Skip background pixels
00167         {
00168             *id = 0;
00169             continue;
00170         }
00171   */      
00172         delta_y = current - (*it).first;
00173         distance = delta_y * delta_y + (*it).second;
00174         *id = distance;
00175     }
00176 
00177 }
00178 
00179 template <class SrcIterator, class SrcAccessor,
00180           class DestIterator, class DestAccessor>
00181 inline void distParabola(triple<SrcIterator, SrcIterator, SrcAccessor> src,
00182                          pair<DestIterator, DestAccessor> dest, float sigma )
00183 {
00184     distParabola(src.first, src.second, src.third,
00185                  dest.first, dest.second, sigma);
00186 }
00187 
00188 /********************************************************/
00189 /*                                                      */
00190 /*        internalSeparableMultiArrayDistTmp            */
00191 /*                                                      */
00192 /********************************************************/
00193 
00194 template <class SrcIterator, class SrcShape, class SrcAccessor,
00195           class DestIterator, class DestAccessor>
00196 void internalSeparableMultiArrayDistTmp(
00197                       SrcIterator si, SrcShape const & shape, SrcAccessor src,
00198                       DestIterator di, DestAccessor dest, float sigma, bool invert)
00199 {
00200     // Sigma is the spread of the parabolas and is only used for ND morphology. When
00201     // calculating the distance transform, it is set to 1
00202     enum { N = 1 + SrcIterator::level };
00203 
00204     // we need the Promote type here if we want to invert the image (dilation)
00205     typedef typename NumericTraits<typename DestAccessor::value_type>::Promote TmpType;
00206     
00207     // temporary array to hold the current line to enable in-place operation
00208     ArrayVector<TmpType> tmp( shape[0] );
00209 
00210     typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
00211     typedef MultiArrayNavigator<DestIterator, N> DNavigator;
00212     
00213     
00214     // only operate on first dimension here
00215     SNavigator snav( si, shape, 0 );
00216     DNavigator dnav( di, shape, 0 );
00217 
00218     using namespace vigra::functor;
00219 
00220     for( ; snav.hasMore(); snav++, dnav++ )
00221     {
00222             // first copy source to temp for maximum cache efficiency
00223             // Invert the values if necessary. Only needed for grayscale morphology
00224             if(invert)
00225                 transformLine( snav.begin(), snav.end(), src, tmp.begin(),
00226                                typename AccessorTraits<TmpType>::default_accessor(), 
00227                                Param(NumericTraits<TmpType>::zero())-Arg1());
00228             else
00229                 copyLine( snav.begin(), snav.end(), src, tmp.begin(),
00230                           typename AccessorTraits<TmpType>::default_accessor() );
00231 
00232             detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
00233                           typename AccessorTraits<TmpType>::default_const_accessor()),
00234                           destIter( dnav.begin(), dest ), sigma );
00235     }
00236 
00237     // operate on further dimensions
00238     for( int d = 1; d < N; ++d )
00239     {
00240         DNavigator dnav( di, shape, d );
00241 
00242         tmp.resize( shape[d] );
00243 
00244         for( ; dnav.hasMore(); dnav++ )
00245         {
00246              // first copy source to temp for maximum cache efficiency
00247              copyLine( dnav.begin(), dnav.end(), dest,
00248                        tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
00249 
00250              detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
00251                            typename AccessorTraits<TmpType>::default_const_accessor()),
00252                            destIter( dnav.begin(), dest ), sigma );
00253         }
00254     }
00255     if(invert) transformMultiArray( di, shape, dest, di, dest, -Arg1());
00256 }
00257 
00258 template <class SrcIterator, class SrcShape, class SrcAccessor,
00259           class DestIterator, class DestAccessor>
00260 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
00261                                                 DestIterator di, DestAccessor dest, float sigma)
00262 {
00263     internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigma, false );
00264 }
00265 
00266 template <class SrcIterator, class SrcShape, class SrcAccessor,
00267           class DestIterator, class DestAccessor>
00268 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
00269                                                 DestIterator di, DestAccessor dest)
00270 {
00271     internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, 1, false );
00272 }
00273 
00274 } // namespace detail
00275 
00276 /** \addtogroup MultiArrayDistanceTransform Euclidean distance transform for multi-dimensional arrays.
00277 
00278     These functions perform the Euclidean distance transform an arbitrary dimensional
00279     array that is specified by iterators (compatible to \ref MultiIteratorPage)
00280     and shape objects. It can therefore be applied to a wide range of data structures
00281     (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
00282 */
00283 //@{
00284 
00285 /********************************************************/
00286 /*                                                      */
00287 /*             separableMultiDistSquared                */
00288 /*                                                      */
00289 /********************************************************/
00290 
00291 /** \brief Euclidean distance squared on multi-dimensional arrays.
00292 
00293     <b> Declarations:</b>
00294 
00295     pass arguments explicitly:
00296     \code
00297     namespace vigra {
00298         // apply the same kernel to all dimensions
00299         template <class SrcIterator, class SrcShape, class SrcAccessor,
00300                   class DestIterator, class DestAccessor>
00301         void
00302         separableMultiDistSquared(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00303                                     DestIterator diter, DestAccessor dest, bool background);
00304 
00305     }
00306     \endcode
00307 
00308     use argument objects in conjunction with \ref ArgumentObjectFactories :
00309     \code
00310     namespace vigra {
00311         template <class SrcIterator, class SrcShape, class SrcAccessor,
00312                   class DestIterator, class DestAccessor>
00313         void
00314         separableMultiDistance(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00315                                     pair<DestIterator, DestAccessor> const & dest,
00316                                     bool background);
00317 
00318     }
00319     \endcode
00320 
00321     This function performs a Euclidean distance squared transform on the given
00322     multi-dimensional array. Both source and destination
00323     arrays are represented by iterators, shape objects and accessors.
00324     The destination array is required to already have the correct size.
00325 
00326     This function expects a mask as its source, where background pixels are 
00327     marked as zero, and non-background pixels as non-zero. If the parameter 
00328     <i>background</i> is true, then the squared distance of all background
00329     pixels to the nearest object is calculated. Otherwise, the distance of all
00330     object pixels to the nearest background pixel is calculated.
00331 
00332     This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
00333     A full-sized internal array is only allocated if working on the destination
00334     array directly would cause overflow errors (i.e. if
00335     <tt> typeid(typename DestAccessor::value_type) < N * M*M</tt>, where M is the
00336     size of the largest dimension of the array.
00337 
00338     <b> Usage:</b>
00339 
00340     <b>\#include</b> <<a href="multi__distance_8hxx-source.html">vigra/multi_distance.hxx</a>>
00341 
00342     \code
00343     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00344     MultiArray<3, unsigned char> source(shape);
00345     MultiArray<3, unsigned int> dest(shape);
00346     ...
00347 
00348     // Calculate Euclidean distance squared for all background pixels 
00349     separableMultiDistSquared(srcMultiArrayRange(source), destMultiArray(dest), true);
00350     \endcode
00351 
00352     \see vigra::distanceTransform()
00353 */
00354 doxygen_overloaded_function(template <...> void separableMultiDistSquared)
00355 
00356 template <class SrcIterator, class SrcShape, class SrcAccessor,
00357           class DestIterator, class DestAccessor>
00358 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00359                                 DestIterator d, DestAccessor dest, bool background)
00360 {
00361     typedef typename NumericTraits<typename DestAccessor::value_type>::ValueType DestType;
00362     typedef typename NumericTraits<typename DestAccessor::value_type>::Promote TmpType;
00363     DestType MaxValue = NumericTraits<DestType>::max();
00364     enum { N = 1 + SrcIterator::level };
00365 
00366     int MaxDim = 0;
00367     for( int i=0; i<N; i++)
00368         if(MaxDim < shape[i]) MaxDim = shape[i];
00369     int MaxDist = MaxDim*MaxDim;
00370 
00371     using namespace vigra::functor;
00372    
00373     if(N*MaxDim*MaxDim > MaxValue) // need a temporary array to avoid overflows
00374     {
00375         // Threshold the values so all objects have infinity value in the beginning
00376         MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
00377         if(background == true)
00378             transformMultiArray( s, shape, src, tmpArray.traverser_begin(),
00379                                  typename AccessorTraits<TmpType>::default_accessor(),
00380                                  ifThenElse( Arg1() == Param(0), Param(MaxDist), Param(0) ));
00381         else
00382             transformMultiArray( s, shape, src, tmpArray.traverser_begin(),
00383                                  typename AccessorTraits<TmpType>::default_accessor(),
00384                                  ifThenElse( Arg1() != Param(0), Param(MaxDist), Param(0) ));
00385         
00386         detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(), 
00387                 shape, typename AccessorTraits<TmpType>::default_accessor(),
00388                 tmpArray.traverser_begin(), 
00389                 typename AccessorTraits<TmpType>::default_accessor());
00390         
00391         //copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
00392         transformMultiArray( tmpArray.traverser_begin(), shape,
00393                              typename AccessorTraits<TmpType>::default_accessor(), d, dest,
00394                              ifThenElse( Arg1() > Param(MaxValue), Param(MaxValue), Arg1() ) );
00395               
00396     }
00397     else        // work directly on the destination array    
00398     {
00399         // Threshold the values so all objects have infinity value in the beginning
00400         if(background == true)
00401             transformMultiArray( s, shape, src, d, dest,
00402                                  ifThenElse( Arg1() == Param(0), Param(MaxDist), Param(0) ));
00403         else
00404             transformMultiArray( s, shape, src, d, dest, 
00405                                  ifThenElse( Arg1() != Param(0), Param(MaxDist), Param(0) ));
00406      
00407         detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest);
00408     }
00409 }
00410 
00411 template <class SrcIterator, class SrcShape, class SrcAccessor,
00412           class DestIterator, class DestAccessor>
00413 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00414                                        pair<DestIterator, DestAccessor> const & dest, bool background)
00415 {
00416     separableMultiDistSquared( source.first, source.second, source.third,
00417                                dest.first, dest.second, background );
00418 }
00419 
00420 /********************************************************/
00421 /*                                                      */
00422 /*             separableMultiDistance                   */
00423 /*                                                      */
00424 /********************************************************/
00425 
00426 /** \brief Euclidean distance on multi-dimensional arrays.
00427 
00428     <b> Declarations:</b>
00429 
00430     pass arguments explicitly:
00431     \code
00432     namespace vigra {
00433         // apply the same kernel to all dimensions
00434         template <class SrcIterator, class SrcShape, class SrcAccessor,
00435                   class DestIterator, class DestAccessor>
00436         void
00437         separableMultiDistance(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
00438                                     DestIterator diter, DestAccessor dest, bool background);
00439 
00440     }
00441     \endcode
00442 
00443     use argument objects in conjunction with \ref ArgumentObjectFactories :
00444     \code
00445     namespace vigra {
00446         template <class SrcIterator, class SrcShape, class SrcAccessor,
00447                   class DestIterator, class DestAccessor>
00448         void
00449         separableMultiDistance(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00450                                     pair<DestIterator, DestAccessor> const & dest,
00451                                     bool background);
00452 
00453     }
00454     \endcode
00455 
00456     This function performs a Euclidean distance transform on the given
00457     multi-dimensional array. Both source and destination
00458     arrays are represented by iterators, shape objects and accessors.
00459     The destination array is required to already have the correct size.
00460 
00461     This function expects a mask as its source, where background pixels are 
00462     marked as zero, and non-background pixels as non-zero. If the parameter 
00463     <i>background</i> is true, then the squared distance of all background
00464     pixels to the nearest object is calculated. Otherwise, the distance of all
00465     object pixels to the nearest background pixel is calculated.
00466 
00467     This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
00468     A full-sized internal array is only allocated if working on the destination
00469     array directly would cause overflow errors (i.e. if
00470     <tt> typeid(typename DestAccessor::value_type) < N * M*M</tt>, where M is the
00471     size of the largest dimension of the array.
00472 
00473     <b> Usage:</b>
00474 
00475     <b>\#include</b> <<a href="multi__distance_8hxx-source.html">vigra/multi_distance.hxx</a>>
00476 
00477     \code
00478     MultiArray<3, unsigned char>::size_type shape(width, height, depth);
00479     MultiArray<3, unsigned char> source(shape);
00480     MultiArray<3, unsigned float> dest(shape);
00481     ...
00482 
00483     // Calculate Euclidean distance squared for all background pixels 
00484     separableMultiDistance(srcMultiArrayRange(source), destMultiArray(dest), true);
00485     \endcode
00486 
00487     \see vigra::distanceTransform()
00488 */
00489 doxygen_overloaded_function(template <...> void separableMultiDistance)
00490 
00491 template <class SrcIterator, class SrcShape, class SrcAccessor,
00492           class DestIterator, class DestAccessor>
00493 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
00494                              DestIterator d, DestAccessor dest, bool background)
00495 {
00496     separableMultiDistSquared( s, shape, src, d, dest, background);
00497     
00498     // Finally, calculate the square root of the distances
00499     transformMultiArray( d, shape, dest, d, dest, (double(*)(double))&std::sqrt );
00500 }
00501 
00502 template <class SrcIterator, class SrcShape, class SrcAccessor,
00503           class DestIterator, class DestAccessor>
00504 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
00505                                     pair<DestIterator, DestAccessor> const & dest, bool background)
00506 {
00507     separableMultiDistance( source.first, source.second, source.third,
00508                             dest.first, dest.second, background );
00509 }
00510 
00511 //@}
00512 
00513 } //-- namespace vigra
00514 
00515 
00516 #endif        //-- VIGRA_MULTI_DISTANCE_HXX

© 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)