[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

multi_distance.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn */
4 /* and Ullrich Koethe */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* The VIGRA Website is */
8 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
9 /* Please direct questions, bug reports, and contributions to */
10 /* ullrich.koethe@iwr.uni-heidelberg.de or */
11 /* vigra@informatik.uni-hamburg.de */
12 /* */
13 /* Permission is hereby granted, free of charge, to any person */
14 /* obtaining a copy of this software and associated documentation */
15 /* files (the "Software"), to deal in the Software without */
16 /* restriction, including without limitation the rights to use, */
17 /* copy, modify, merge, publish, distribute, sublicense, and/or */
18 /* sell copies of the Software, and to permit persons to whom the */
19 /* Software is furnished to do so, subject to the following */
20 /* conditions: */
21 /* */
22 /* The above copyright notice and this permission notice shall be */
23 /* included in all copies or substantial portions of the */
24 /* Software. */
25 /* */
26 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
27 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
28 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
29 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
30 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
31 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
32 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
33 /* OTHER DEALINGS IN THE SOFTWARE. */
34 /* */
35 /************************************************************************/
36 
37 #ifndef VIGRA_MULTI_DISTANCE_HXX
38 #define VIGRA_MULTI_DISTANCE_HXX
39 
40 #include <vector>
41 #include <functional>
42 #include "array_vector.hxx"
43 #include "multi_array.hxx"
44 #include "accessor.hxx"
45 #include "numerictraits.hxx"
46 #include "navigator.hxx"
47 #include "metaprogramming.hxx"
48 #include "multi_pointoperators.hxx"
49 #include "functorexpression.hxx"
50 
51 namespace vigra
52 {
53 
54 namespace detail
55 {
56 
57 template <class Value>
58 struct DistParabolaStackEntry
59 {
60  double left, center, right;
61  Value prevVal;
62 
63  DistParabolaStackEntry(Value const & p, double l, double c, double r)
64  : left(l), center(c), right(r), prevVal(p)
65  {}
66 };
67 
68 /********************************************************/
69 /* */
70 /* distParabola */
71 /* */
72 /* Version with sigma (parabola spread) for morphology */
73 /* */
74 /********************************************************/
75 
76 template <class SrcIterator, class SrcAccessor,
77  class DestIterator, class DestAccessor >
78 void distParabola(SrcIterator is, SrcIterator iend, SrcAccessor sa,
79  DestIterator id, DestAccessor da, double sigma )
80 {
81  // We assume that the data in the input is distance squared and treat it as such
82  double w = iend - is;
83  double sigma2 = sigma * sigma;
84  double sigma22 = 2.0 * sigma2;
85 
86  typedef typename SrcAccessor::value_type SrcType;
87  typedef DistParabolaStackEntry<SrcType> Influence;
88  std::vector<Influence> _stack;
89  _stack.push_back(Influence(sa(is), 0.0, 0.0, w));
90 
91  ++is;
92  double current = 1.0;
93  while(current < w )
94  {
95  Influence & s = _stack.back();
96  double diff = current - s.center;
97  double intersection = current + (sa(is) - s.prevVal - sigma2*sq(diff)) / (sigma22 * diff);
98 
99  if( intersection < s.left) // previous point has no influence
100  {
101  _stack.pop_back();
102  if(_stack.empty())
103  {
104  _stack.push_back(Influence(sa(is), 0.0, current, w));
105  }
106  else
107  {
108  continue; // try new top of stack without advancing current
109  }
110  }
111  else if(intersection < s.right)
112  {
113  s.right = intersection;
114  _stack.push_back(Influence(sa(is), intersection, current, w));
115  }
116  ++is;
117  ++current;
118  }
119 
120  // Now we have the stack indicating which rows are influenced by (and therefore
121  // closest to) which row. We can go through the stack and calculate the
122  // distance squared for each element of the column.
123  typename std::vector<Influence>::iterator it = _stack.begin();
124  for(current = 0.0; current < w; ++current, ++id)
125  {
126  while( current >= it->right)
127  ++it;
128  da.set(sigma2 * sq(current - it->center) + it->prevVal, id);
129  }
130 }
131 
132 template <class SrcIterator, class SrcAccessor,
133  class DestIterator, class DestAccessor>
134 inline void distParabola(triple<SrcIterator, SrcIterator, SrcAccessor> src,
135  pair<DestIterator, DestAccessor> dest, double sigma)
136 {
137  distParabola(src.first, src.second, src.third,
138  dest.first, dest.second, sigma);
139 }
140 
141 /********************************************************/
142 /* */
143 /* internalSeparableMultiArrayDistTmp */
144 /* */
145 /********************************************************/
146 
147 template <class SrcIterator, class SrcShape, class SrcAccessor,
148  class DestIterator, class DestAccessor, class Array>
149 void internalSeparableMultiArrayDistTmp(
150  SrcIterator si, SrcShape const & shape, SrcAccessor src,
151  DestIterator di, DestAccessor dest, Array const & sigmas, bool invert)
152 {
153  // Sigma is the spread of the parabolas. It determines the structuring element size
154  // for ND morphology. When calculating the distance transforms, sigma is usually set to 1,
155  // unless one wants to account for anisotropic pixel pitch
156  enum { N = SrcShape::static_size};
157 
158  // we need the Promote type here if we want to invert the image (dilation)
159  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
160 
161  // temporary array to hold the current line to enable in-place operation
162  ArrayVector<TmpType> tmp( shape[0] );
163 
164  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
165  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
166 
167 
168  // only operate on first dimension here
169  SNavigator snav( si, shape, 0 );
170  DNavigator dnav( di, shape, 0 );
171 
172  using namespace vigra::functor;
173 
174  for( ; snav.hasMore(); snav++, dnav++ )
175  {
176  // first copy source to temp for maximum cache efficiency
177  // Invert the values if necessary. Only needed for grayscale morphology
178  if(invert)
179  transformLine( snav.begin(), snav.end(), src, tmp.begin(),
180  typename AccessorTraits<TmpType>::default_accessor(),
181  Param(NumericTraits<TmpType>::zero())-Arg1());
182  else
183  copyLine( snav.begin(), snav.end(), src, tmp.begin(),
184  typename AccessorTraits<TmpType>::default_accessor() );
185 
186  detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
187  typename AccessorTraits<TmpType>::default_const_accessor()),
188  destIter( dnav.begin(), dest ), sigmas[0] );
189  }
190 
191  // operate on further dimensions
192  for( int d = 1; d < N; ++d )
193  {
194  DNavigator dnav( di, shape, d );
195 
196  tmp.resize( shape[d] );
197 
198 
199  for( ; dnav.hasMore(); dnav++ )
200  {
201  // first copy source to temp for maximum cache efficiency
202  copyLine( dnav.begin(), dnav.end(), dest,
203  tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
204 
205  detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
206  typename AccessorTraits<TmpType>::default_const_accessor()),
207  destIter( dnav.begin(), dest ), sigmas[d] );
208  }
209  }
210  if(invert) transformMultiArray( di, shape, dest, di, dest, -Arg1());
211 }
212 
213 template <class SrcIterator, class SrcShape, class SrcAccessor,
214  class DestIterator, class DestAccessor, class Array>
215 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
216  DestIterator di, DestAccessor dest, Array const & sigmas)
217 {
218  internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas, false );
219 }
220 
221 template <class SrcIterator, class SrcShape, class SrcAccessor,
222  class DestIterator, class DestAccessor>
223 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
224  DestIterator di, DestAccessor dest)
225 {
226  ArrayVector<double> sigmas(shape.size(), 1.0);
227  internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas, false );
228 }
229 
230 } // namespace detail
231 
232 /** \addtogroup MultiArrayDistanceTransform Euclidean distance transform for multi-dimensional arrays.
233 
234  These functions perform the Euclidean distance transform an arbitrary dimensional
235  array that is specified by iterators (compatible to \ref MultiIteratorPage)
236  and shape objects. It can therefore be applied to a wide range of data structures
237  (\ref vigra::MultiArrayView, \ref vigra::MultiArray etc.).
238 */
239 //@{
240 
241 /********************************************************/
242 /* */
243 /* separableMultiDistSquared */
244 /* */
245 /********************************************************/
246 
247 /** \brief Euclidean distance squared on multi-dimensional arrays.
248 
249  The algorithm is taken from Donald Bailey: "An Efficient Euclidean Distance Transform",
250  Proc. IWCIA'04, Springer LNCS 3322, 2004.
251 
252  <b> Declarations:</b>
253 
254  pass arguments explicitly:
255  \code
256  namespace vigra {
257  // explicitly specify pixel pitch for each coordinate
258  template <class SrcIterator, class SrcShape, class SrcAccessor,
259  class DestIterator, class DestAccessor, class Array>
260  void
261  separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
262  DestIterator d, DestAccessor dest,
263  bool background,
264  Array const & pixelPitch);
265 
266  // use default pixel pitch = 1.0 for each coordinate
267  template <class SrcIterator, class SrcShape, class SrcAccessor,
268  class DestIterator, class DestAccessor>
269  void
270  separableMultiDistSquared(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
271  DestIterator diter, DestAccessor dest,
272  bool background);
273 
274  }
275  \endcode
276 
277  use argument objects in conjunction with \ref ArgumentObjectFactories :
278  \code
279  namespace vigra {
280  // explicitly specify pixel pitch for each coordinate
281  template <class SrcIterator, class SrcShape, class SrcAccessor,
282  class DestIterator, class DestAccessor, class Array>
283  void
284  separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
285  pair<DestIterator, DestAccessor> const & dest,
286  bool background,
287  Array const & pixelPitch);
288 
289  // use default pixel pitch = 1.0 for each coordinate
290  template <class SrcIterator, class SrcShape, class SrcAccessor,
291  class DestIterator, class DestAccessor>
292  void
293  separableMultiDistSquared(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
294  pair<DestIterator, DestAccessor> const & dest,
295  bool background);
296 
297  }
298  \endcode
299 
300  This function performs a squared Euclidean squared distance transform on the given
301  multi-dimensional array. Both source and destination
302  arrays are represented by iterators, shape objects and accessors.
303  The destination array is required to already have the correct size.
304 
305  This function expects a mask as its source, where background pixels are
306  marked as zero, and non-background pixels as non-zero. If the parameter
307  <i>background</i> is true, then the squared distance of all background
308  pixels to the nearest object is calculated. Otherwise, the distance of all
309  object pixels to the nearest background pixel is calculated.
310 
311  Optionally, one can pass an array that specifies the pixel pitch in each direction.
312  This is necessary when the data have non-uniform resolution (as is common in confocal
313  microscopy, for example).
314 
315  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
316  A full-sized internal array is only allocated if working on the destination
317  array directly would cause overflow errors (i.e. if
318  <tt> NumericTraits<typename DestAccessor::value_type>::max() < N * M*M</tt>, where M is the
319  size of the largest dimension of the array.
320 
321  <b> Usage:</b>
322 
323  <b>\#include</b> <vigra/multi_distance.hxx>
324 
325  \code
326  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
327  MultiArray<3, unsigned char> source(shape);
328  MultiArray<3, unsigned int> dest(shape);
329  ...
330 
331  // Calculate Euclidean distance squared for all background pixels
332  separableMultiDistSquared(srcMultiArrayRange(source), destMultiArray(dest), true);
333  \endcode
334 
335  \see vigra::distanceTransform(), vigra::separableMultiDistance()
336 */
337 doxygen_overloaded_function(template <...> void separableMultiDistSquared)
338 
339 template <class SrcIterator, class SrcShape, class SrcAccessor,
340  class DestIterator, class DestAccessor, class Array>
341 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
342  DestIterator d, DestAccessor dest, bool background,
343  Array const & pixelPitch)
344 {
345  int N = shape.size();
346 
347  typedef typename SrcAccessor::value_type SrcType;
348  typedef typename DestAccessor::value_type DestType;
349  typedef typename NumericTraits<DestType>::RealPromote Real;
350 
351  SrcType zero = NumericTraits<SrcType>::zero();
352 
353  double dmax = 0.0;
354  bool pixelPitchIsReal = false;
355  for( int k=0; k<N; ++k)
356  {
357  if(int(pixelPitch[k]) != pixelPitch[k])
358  pixelPitchIsReal = true;
359  dmax += sq(pixelPitch[k]*shape[k]);
360  }
361 
362  using namespace vigra::functor;
363 
364  if(dmax > NumericTraits<DestType>::toRealPromote(NumericTraits<DestType>::max())
365  || pixelPitchIsReal) // need a temporary array to avoid overflows
366  {
367  // Threshold the values so all objects have infinity value in the beginning
368  Real maxDist = (Real)dmax, rzero = (Real)0.0;
369  MultiArray<SrcShape::static_size, Real> tmpArray(shape);
370  if(background == true)
371  transformMultiArray( s, shape, src,
372  tmpArray.traverser_begin(), typename AccessorTraits<Real>::default_accessor(),
373  ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
374  else
375  transformMultiArray( s, shape, src,
376  tmpArray.traverser_begin(), typename AccessorTraits<Real>::default_accessor(),
377  ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
378 
379  detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(),
380  shape, typename AccessorTraits<Real>::default_accessor(),
381  tmpArray.traverser_begin(),
382  typename AccessorTraits<Real>::default_accessor(), pixelPitch);
383 
384  copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
385  }
386  else // work directly on the destination array
387  {
388  // Threshold the values so all objects have infinity value in the beginning
389  DestType maxDist = DestType(std::ceil(dmax)), rzero = (DestType)0;
390  if(background == true)
391  transformMultiArray( s, shape, src, d, dest,
392  ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
393  else
394  transformMultiArray( s, shape, src, d, dest,
395  ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
396 
397  detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest, pixelPitch);
398  }
399 }
400 
401 template <class SrcIterator, class SrcShape, class SrcAccessor,
402  class DestIterator, class DestAccessor, class Array>
403 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
404  pair<DestIterator, DestAccessor> const & dest, bool background,
405  Array const & pixelPitch)
406 {
407  separableMultiDistSquared( source.first, source.second, source.third,
408  dest.first, dest.second, background, pixelPitch );
409 }
410 
411 template <class SrcIterator, class SrcShape, class SrcAccessor,
412  class DestIterator, class DestAccessor>
413 inline
414 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
415  DestIterator d, DestAccessor dest, bool background)
416 {
417  ArrayVector<double> pixelPitch(shape.size(), 1.0);
418  separableMultiDistSquared( s, shape, src, d, dest, background, pixelPitch );
419 }
420 
421 template <class SrcIterator, class SrcShape, class SrcAccessor,
422  class DestIterator, class DestAccessor>
423 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
424  pair<DestIterator, DestAccessor> const & dest, bool background)
425 {
426  separableMultiDistSquared( source.first, source.second, source.third,
427  dest.first, dest.second, background );
428 }
429 
430 /********************************************************/
431 /* */
432 /* separableMultiDistance */
433 /* */
434 /********************************************************/
435 
436 /** \brief Euclidean distance on multi-dimensional arrays.
437 
438  <b> Declarations:</b>
439 
440  pass arguments explicitly:
441  \code
442  namespace vigra {
443  // explicitly specify pixel pitch for each coordinate
444  template <class SrcIterator, class SrcShape, class SrcAccessor,
445  class DestIterator, class DestAccessor, class Array>
446  void
447  separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
448  DestIterator d, DestAccessor dest,
449  bool background,
450  Array const & pixelPitch);
451 
452  // use default pixel pitch = 1.0 for each coordinate
453  template <class SrcIterator, class SrcShape, class SrcAccessor,
454  class DestIterator, class DestAccessor>
455  void
456  separableMultiDistance(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
457  DestIterator diter, DestAccessor dest,
458  bool background);
459 
460  }
461  \endcode
462 
463  use argument objects in conjunction with \ref ArgumentObjectFactories :
464  \code
465  namespace vigra {
466  // explicitly specify pixel pitch for each coordinate
467  template <class SrcIterator, class SrcShape, class SrcAccessor,
468  class DestIterator, class DestAccessor, class Array>
469  void
470  separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
471  pair<DestIterator, DestAccessor> const & dest,
472  bool background,
473  Array const & pixelPitch);
474 
475  // use default pixel pitch = 1.0 for each coordinate
476  template <class SrcIterator, class SrcShape, class SrcAccessor,
477  class DestIterator, class DestAccessor>
478  void
479  separableMultiDistance(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
480  pair<DestIterator, DestAccessor> const & dest,
481  bool background);
482 
483  }
484  \endcode
485 
486  This function performs a Euclidean distance transform on the given
487  multi-dimensional array. It simply calls \ref separableMultiDistSquared()
488  and takes the pixel-wise square root of the result. See \ref separableMultiDistSquared()
489  for more documentation.
490 
491  <b> Usage:</b>
492 
493  <b>\#include</b> <vigra/multi_distance.hxx>
494 
495  \code
496  MultiArray<3, unsigned char>::size_type shape(width, height, depth);
497  MultiArray<3, unsigned char> source(shape);
498  MultiArray<3, float> dest(shape);
499  ...
500 
501  // Calculate Euclidean distance squared for all background pixels
502  separableMultiDistance(srcMultiArrayRange(source), destMultiArray(dest), true);
503  \endcode
504 
505  \see vigra::distanceTransform(), vigra::separableMultiDistSquared()
506 */
507 doxygen_overloaded_function(template <...> void separableMultiDistance)
508 
509 template <class SrcIterator, class SrcShape, class SrcAccessor,
510  class DestIterator, class DestAccessor, class Array>
511 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
512  DestIterator d, DestAccessor dest, bool background,
513  Array const & pixelPitch)
514 {
515  separableMultiDistSquared( s, shape, src, d, dest, background, pixelPitch);
516 
517  // Finally, calculate the square root of the distances
518  using namespace vigra::functor;
519 
520  transformMultiArray( d, shape, dest, d, dest, sqrt(Arg1()) );
521 }
522 
523 template <class SrcIterator, class SrcShape, class SrcAccessor,
524  class DestIterator, class DestAccessor>
525 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
526  DestIterator d, DestAccessor dest, bool background)
527 {
528  separableMultiDistSquared( s, shape, src, d, dest, background);
529 
530  // Finally, calculate the square root of the distances
531  using namespace vigra::functor;
532 
533  transformMultiArray( d, shape, dest, d, dest, sqrt(Arg1()) );
534 }
535 
536 template <class SrcIterator, class SrcShape, class SrcAccessor,
537  class DestIterator, class DestAccessor, class Array>
538 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
539  pair<DestIterator, DestAccessor> const & dest, bool background,
540  Array const & pixelPitch)
541 {
542  separableMultiDistance( source.first, source.second, source.third,
543  dest.first, dest.second, background, pixelPitch );
544 }
545 
546 template <class SrcIterator, class SrcShape, class SrcAccessor,
547  class DestIterator, class DestAccessor>
548 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
549  pair<DestIterator, DestAccessor> const & dest, bool background)
550 {
551  separableMultiDistance( source.first, source.second, source.third,
552  dest.first, dest.second, background );
553 }
554 
555 //@}
556 
557 } //-- namespace vigra
558 
559 
560 #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.8.0 (Wed Sep 26 2012)