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