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