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