[ 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 /* ( 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) |
html generated using doxygen and Python
|