[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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) |
html generated using doxygen and Python
|