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

multi_distance.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn, */
4 /* Philipp Schubert 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 #include "multi_gridgraph.hxx" //for boundaryGraph & boundaryMultiDistance
52 #include "union_find.hxx" //for boundaryGraph & boundaryMultiDistance
53 
54 namespace vigra
55 {
56 
57 namespace detail
58 {
59 
60 template <class Value>
61 struct DistParabolaStackEntry
62 {
63  double left, center, right;
64  Value apex_height;
65 
66  DistParabolaStackEntry(Value const & p, double l, double c, double r)
67  : left(l), center(c), right(r), apex_height(p)
68  {}
69 };
70 
71 /********************************************************/
72 /* */
73 /* distParabola */
74 /* */
75 /* Version with sigma (parabola spread) for morphology */
76 /* */
77 /********************************************************/
78 
79 template <class SrcIterator, class SrcAccessor,
80  class DestIterator, class DestAccessor >
81 void distParabola(SrcIterator is, SrcIterator iend, SrcAccessor sa,
82  DestIterator id, DestAccessor da, double sigma )
83 {
84  // We assume that the data in the input is distance squared and treat it as such
85  double w = iend - is;
86  if(w <= 0)
87  return;
88 
89  double sigma2 = sigma * sigma;
90  double sigma22 = 2.0 * sigma2;
91 
92  typedef typename SrcAccessor::value_type SrcType;
93  typedef DistParabolaStackEntry<SrcType> Influence;
94  std::vector<Influence> _stack;
95  _stack.push_back(Influence(sa(is), 0.0, 0.0, w));
96 
97  ++is;
98  double current = 1.0;
99  for(;current < w; ++is, ++current)
100  {
101  double intersection;
102 
103  while(true)
104  {
105  Influence & s = _stack.back();
106  double diff = current - s.center;
107  intersection = current + (sa(is) - s.apex_height - sigma2*sq(diff)) / (sigma22 * diff);
108 
109  if( intersection < s.left) // previous point has no influence
110  {
111  _stack.pop_back();
112  if(!_stack.empty())
113  continue; // try new top of stack without advancing current
114  else
115  intersection = 0.0;
116  }
117  else if(intersection < s.right)
118  {
119  s.right = intersection;
120  }
121  break;
122  }
123  _stack.push_back(Influence(sa(is), intersection, current, w));
124  }
125 
126  // Now we have the stack indicating which rows are influenced by (and therefore
127  // closest to) which row. We can go through the stack and calculate the
128  // distance squared for each element of the column.
129  typename std::vector<Influence>::iterator it = _stack.begin();
130  for(current = 0.0; current < w; ++current, ++id)
131  {
132  while( current >= it->right)
133  ++it;
134  da.set(sigma2 * sq(current - it->center) + it->apex_height, id);
135  }
136 }
137 
138 template <class SrcIterator, class SrcAccessor,
139  class DestIterator, class DestAccessor>
140 inline void distParabola(triple<SrcIterator, SrcIterator, SrcAccessor> src,
141  pair<DestIterator, DestAccessor> dest, double sigma)
142 {
143  distParabola(src.first, src.second, src.third,
144  dest.first, dest.second, sigma);
145 }
146 
147 /********************************************************/
148 /* */
149 /* internalSeparableMultiArrayDistTmp */
150 /* */
151 /********************************************************/
152 
153 template <class SrcIterator, class SrcShape, class SrcAccessor,
154  class DestIterator, class DestAccessor, class Array>
155 void internalSeparableMultiArrayDistTmp(
156  SrcIterator si, SrcShape const & shape, SrcAccessor src,
157  DestIterator di, DestAccessor dest, Array const & sigmas, bool invert)
158 {
159  // Sigma is the spread of the parabolas. It determines the structuring element size
160  // for ND morphology. When calculating the distance transforms, sigma is usually set to 1,
161  // unless one wants to account for anisotropic pixel pitch
162  enum { N = SrcShape::static_size};
163 
164  // we need the Promote type here if we want to invert the image (dilation)
165  typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
166 
167  // temporary array to hold the current line to enable in-place operation
168  ArrayVector<TmpType> tmp( shape[0] );
169 
170  typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
171  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
172 
173 
174  // only operate on first dimension here
175  SNavigator snav( si, shape, 0 );
176  DNavigator dnav( di, shape, 0 );
177 
178  using namespace vigra::functor;
179 
180  for( ; snav.hasMore(); snav++, dnav++ )
181  {
182  // first copy source to temp for maximum cache efficiency
183  // Invert the values if necessary. Only needed for grayscale morphology
184  if(invert)
185  transformLine( snav.begin(), snav.end(), src, tmp.begin(),
186  typename AccessorTraits<TmpType>::default_accessor(),
187  Param(NumericTraits<TmpType>::zero())-Arg1());
188  else
189  copyLine( snav.begin(), snav.end(), src, tmp.begin(),
190  typename AccessorTraits<TmpType>::default_accessor() );
191 
192  detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
193  typename AccessorTraits<TmpType>::default_const_accessor()),
194  destIter( dnav.begin(), dest ), sigmas[0] );
195  }
196 
197  // operate on further dimensions
198  for( int d = 1; d < N; ++d )
199  {
200  DNavigator dnav( di, shape, d );
201 
202  tmp.resize( shape[d] );
203 
204 
205  for( ; dnav.hasMore(); dnav++ )
206  {
207  // first copy source to temp for maximum cache efficiency
208  copyLine( dnav.begin(), dnav.end(), dest,
209  tmp.begin(), typename AccessorTraits<TmpType>::default_accessor() );
210 
211  detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
212  typename AccessorTraits<TmpType>::default_const_accessor()),
213  destIter( dnav.begin(), dest ), sigmas[d] );
214  }
215  }
216  if(invert) transformMultiArray( di, shape, dest, di, dest, -Arg1());
217 }
218 
219 template <class SrcIterator, class SrcShape, class SrcAccessor,
220  class DestIterator, class DestAccessor, class Array>
221 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
222  DestIterator di, DestAccessor dest, Array const & sigmas)
223 {
224  internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas, false );
225 }
226 
227 template <class SrcIterator, class SrcShape, class SrcAccessor,
228  class DestIterator, class DestAccessor>
229 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
230  DestIterator di, DestAccessor dest)
231 {
232  ArrayVector<double> sigmas(shape.size(), 1.0);
233  internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas, false );
234 }
235 
236 } // namespace detail
237 
238 /** \addtogroup DistanceTransform
239 */
240 //@{
241 
242 /********************************************************/
243 /* */
244 /* separableMultiDistSquared */
245 /* */
246 /********************************************************/
247 
248 /** \brief Euclidean distance squared on multi-dimensional arrays.
249 
250  The algorithm is taken from Donald Bailey: "An Efficient Euclidean Distance Transform",
251  Proc. IWCIA'04, Springer LNCS 3322, 2004.
252 
253  <b> Declarations:</b>
254 
255  pass arbitrary-dimensional array views:
256  \code
257  namespace vigra {
258  // explicitly specify pixel pitch for each coordinate
259  template <unsigned int N, class T1, class S1,
260  class T2, class S2,
261  class Array>
262  void
263  separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
264  MultiArrayView<N, T2, S2> dest,
265  bool background,
266  Array const & pixelPitch);
267 
268  // use default pixel pitch = 1.0 for each coordinate
269  template <unsigned int N, class T1, class S1,
270  class T2, class S2>
271  void
272  separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
273  MultiArrayView<N, T2, S2> dest,
274  bool background);
275  }
276  \endcode
277 
278  \deprecatedAPI{separableMultiDistSquared}
279  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
280  \code
281  namespace vigra {
282  // explicitly specify pixel pitch for each coordinate
283  template <class SrcIterator, class SrcShape, class SrcAccessor,
284  class DestIterator, class DestAccessor, class Array>
285  void
286  separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
287  DestIterator d, DestAccessor dest,
288  bool background,
289  Array const & pixelPitch);
290 
291  // use default pixel pitch = 1.0 for each coordinate
292  template <class SrcIterator, class SrcShape, class SrcAccessor,
293  class DestIterator, class DestAccessor>
294  void
295  separableMultiDistSquared(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
296  DestIterator diter, DestAccessor dest,
297  bool background);
298 
299  }
300  \endcode
301  use argument objects in conjunction with \ref ArgumentObjectFactories :
302  \code
303  namespace vigra {
304  // explicitly specify pixel pitch for each coordinate
305  template <class SrcIterator, class SrcShape, class SrcAccessor,
306  class DestIterator, class DestAccessor, class Array>
307  void
308  separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
309  pair<DestIterator, DestAccessor> const & dest,
310  bool background,
311  Array const & pixelPitch);
312 
313  // use default pixel pitch = 1.0 for each coordinate
314  template <class SrcIterator, class SrcShape, class SrcAccessor,
315  class DestIterator, class DestAccessor>
316  void
317  separableMultiDistSquared(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
318  pair<DestIterator, DestAccessor> const & dest,
319  bool background);
320 
321  }
322  \endcode
323  \deprecatedEnd
324 
325  This function performs a squared Euclidean squared distance transform on the given
326  multi-dimensional array. Both source and destination
327  arrays are represented by iterators, shape objects and accessors.
328  The destination array is required to already have the correct size.
329 
330  This function expects a mask as its source, where background pixels are
331  marked as zero, and non-background pixels as non-zero. If the parameter
332  <i>background</i> is true, then the squared distance of all background
333  pixels to the nearest object is calculated. Otherwise, the distance of all
334  object pixels to the nearest background pixel is calculated.
335 
336  Optionally, one can pass an array that specifies the pixel pitch in each direction.
337  This is necessary when the data have non-uniform resolution (as is common in confocal
338  microscopy, for example).
339 
340  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
341  A full-sized internal array is only allocated if working on the destination
342  array directly would cause overflow errors (i.e. if
343  <tt> NumericTraits<typename DestAccessor::value_type>::max() < N * M*M</tt>, where M is the
344  size of the largest dimension of the array.
345 
346  <b> Usage:</b>
347 
348  <b>\#include</b> <vigra/multi_distance.hxx><br/>
349  Namespace: vigra
350 
351  \code
352  Shape3 shape(width, height, depth);
353  MultiArray<3, unsigned char> source(shape);
354  MultiArray<3, unsigned int> dest(shape);
355  ...
356 
357  // Calculate Euclidean distance squared for all background pixels
358  separableMultiDistSquared(source, dest, true);
359  \endcode
360 
361  \see vigra::distanceTransform(), vigra::separableMultiDistance()
362 */
363 doxygen_overloaded_function(template <...> void separableMultiDistSquared)
364 
365 template <class SrcIterator, class SrcShape, class SrcAccessor,
366  class DestIterator, class DestAccessor, class Array>
367 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
368  DestIterator d, DestAccessor dest, bool background,
369  Array const & pixelPitch)
370 {
371  int N = shape.size();
372 
373  typedef typename SrcAccessor::value_type SrcType;
374  typedef typename DestAccessor::value_type DestType;
375  typedef typename NumericTraits<DestType>::RealPromote Real;
376 
377  SrcType zero = NumericTraits<SrcType>::zero();
378 
379  double dmax = 0.0;
380  bool pixelPitchIsReal = false;
381  for( int k=0; k<N; ++k)
382  {
383  if(int(pixelPitch[k]) != pixelPitch[k])
384  pixelPitchIsReal = true;
385  dmax += sq(pixelPitch[k]*shape[k]);
386  }
387 
388  using namespace vigra::functor;
389 
390  if(dmax > NumericTraits<DestType>::toRealPromote(NumericTraits<DestType>::max())
391  || pixelPitchIsReal) // need a temporary array to avoid overflows
392  {
393  // Threshold the values so all objects have infinity value in the beginning
394  Real maxDist = (Real)dmax, rzero = (Real)0.0;
395  MultiArray<SrcShape::static_size, Real> tmpArray(shape);
396  if(background == true)
397  transformMultiArray( s, shape, src,
398  tmpArray.traverser_begin(), typename AccessorTraits<Real>::default_accessor(),
399  ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
400  else
401  transformMultiArray( s, shape, src,
402  tmpArray.traverser_begin(), typename AccessorTraits<Real>::default_accessor(),
403  ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
404 
405  detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(),
406  shape, typename AccessorTraits<Real>::default_accessor(),
407  tmpArray.traverser_begin(),
408  typename AccessorTraits<Real>::default_accessor(), pixelPitch);
409 
410  copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
411  }
412  else // work directly on the destination array
413  {
414  // Threshold the values so all objects have infinity value in the beginning
415  DestType maxDist = DestType(std::ceil(dmax)), rzero = (DestType)0;
416  if(background == true)
417  transformMultiArray( s, shape, src, d, dest,
418  ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
419  else
420  transformMultiArray( s, shape, src, d, dest,
421  ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
422 
423  detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest, pixelPitch);
424  }
425 }
426 
427 template <class SrcIterator, class SrcShape, class SrcAccessor,
428  class DestIterator, class DestAccessor>
429 inline
430 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
431  DestIterator d, DestAccessor dest, bool background)
432 {
433  ArrayVector<double> pixelPitch(shape.size(), 1.0);
434  separableMultiDistSquared( s, shape, src, d, dest, background, pixelPitch );
435 }
436 
437 template <class SrcIterator, class SrcShape, class SrcAccessor,
438  class DestIterator, class DestAccessor, class Array>
439 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
440  pair<DestIterator, DestAccessor> const & dest, bool background,
441  Array const & pixelPitch)
442 {
443  separableMultiDistSquared( source.first, source.second, source.third,
444  dest.first, dest.second, background, pixelPitch );
445 }
446 
447 template <class SrcIterator, class SrcShape, class SrcAccessor,
448  class DestIterator, class DestAccessor>
449 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
450  pair<DestIterator, DestAccessor> const & dest, bool background)
451 {
452  separableMultiDistSquared( source.first, source.second, source.third,
453  dest.first, dest.second, background );
454 }
455 
456 template <unsigned int N, class T1, class S1,
457  class T2, class S2,
458  class Array>
459 inline void
460 separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
461  MultiArrayView<N, T2, S2> dest, bool background,
462  Array const & pixelPitch)
463 {
464  vigra_precondition(source.shape() == dest.shape(),
465  "separableMultiDistSquared(): shape mismatch between input and output.");
466  separableMultiDistSquared( srcMultiArrayRange(source),
467  destMultiArray(dest), background, pixelPitch );
468 }
469 
470 template <unsigned int N, class T1, class S1,
471  class T2, class S2>
472 inline void
473 separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
474  MultiArrayView<N, T2, S2> dest, bool background)
475 {
476  vigra_precondition(source.shape() == dest.shape(),
477  "separableMultiDistSquared(): shape mismatch between input and output.");
478  separableMultiDistSquared( srcMultiArrayRange(source),
479  destMultiArray(dest), background );
480 }
481 
482 /********************************************************/
483 /* */
484 /* separableMultiDistance */
485 /* */
486 /********************************************************/
487 
488 /** \brief Euclidean distance on multi-dimensional arrays.
489 
490  <b> Declarations:</b>
491 
492  pass arbitrary-dimensional array views:
493  \code
494  namespace vigra {
495  // explicitly specify pixel pitch for each coordinate
496  template <unsigned int N, class T1, class S1,
497  class T2, class S2, class Array>
498  void
499  separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
500  MultiArrayView<N, T2, S2> dest,
501  bool background,
502  Array const & pixelPitch);
503 
504  // use default pixel pitch = 1.0 for each coordinate
505  template <unsigned int N, class T1, class S1,
506  class T2, class S2>
507  void
508  separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
509  MultiArrayView<N, T2, S2> dest,
510  bool background);
511  }
512  \endcode
513 
514  \deprecatedAPI{separableMultiDistance}
515  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
516  \code
517  namespace vigra {
518  // explicitly specify pixel pitch for each coordinate
519  template <class SrcIterator, class SrcShape, class SrcAccessor,
520  class DestIterator, class DestAccessor, class Array>
521  void
522  separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
523  DestIterator d, DestAccessor dest,
524  bool background,
525  Array const & pixelPitch);
526 
527  // use default pixel pitch = 1.0 for each coordinate
528  template <class SrcIterator, class SrcShape, class SrcAccessor,
529  class DestIterator, class DestAccessor>
530  void
531  separableMultiDistance(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
532  DestIterator diter, DestAccessor dest,
533  bool background);
534 
535  }
536  \endcode
537  use argument objects in conjunction with \ref ArgumentObjectFactories :
538  \code
539  namespace vigra {
540  // explicitly specify pixel pitch for each coordinate
541  template <class SrcIterator, class SrcShape, class SrcAccessor,
542  class DestIterator, class DestAccessor, class Array>
543  void
544  separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
545  pair<DestIterator, DestAccessor> const & dest,
546  bool background,
547  Array const & pixelPitch);
548 
549  // use default pixel pitch = 1.0 for each coordinate
550  template <class SrcIterator, class SrcShape, class SrcAccessor,
551  class DestIterator, class DestAccessor>
552  void
553  separableMultiDistance(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
554  pair<DestIterator, DestAccessor> const & dest,
555  bool background);
556 
557  }
558  \endcode
559  \deprecatedEnd
560 
561  This function performs a Euclidean distance transform on the given
562  multi-dimensional array. It simply calls \ref separableMultiDistSquared()
563  and takes the pixel-wise square root of the result. See \ref separableMultiDistSquared()
564  for more documentation.
565 
566  <b> Usage:</b>
567 
568  <b>\#include</b> <vigra/multi_distance.hxx><br/>
569  Namespace: vigra
570 
571  \code
572  Shape3 shape(width, height, depth);
573  MultiArray<3, unsigned char> source(shape);
574  MultiArray<3, float> dest(shape);
575  ...
576 
577  // Calculate Euclidean distance for all background pixels
578  separableMultiDistance(source, dest, true);
579  \endcode
580 
581  \see vigra::distanceTransform(), vigra::separableMultiDistSquared()
582 */
583 doxygen_overloaded_function(template <...> void separableMultiDistance)
584 
585 template <class SrcIterator, class SrcShape, class SrcAccessor,
586  class DestIterator, class DestAccessor, class Array>
587 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
588  DestIterator d, DestAccessor dest, bool background,
589  Array const & pixelPitch)
590 {
591  separableMultiDistSquared( s, shape, src, d, dest, background, pixelPitch);
592 
593  // Finally, calculate the square root of the distances
594  using namespace vigra::functor;
595 
596  transformMultiArray( d, shape, dest, d, dest, sqrt(Arg1()) );
597 }
598 
599 template <class SrcIterator, class SrcShape, class SrcAccessor,
600  class DestIterator, class DestAccessor>
601 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
602  DestIterator d, DestAccessor dest, bool background)
603 {
604  separableMultiDistSquared( s, shape, src, d, dest, background);
605 
606  // Finally, calculate the square root of the distances
607  using namespace vigra::functor;
608 
609  transformMultiArray( d, shape, dest, d, dest, sqrt(Arg1()) );
610 }
611 
612 template <class SrcIterator, class SrcShape, class SrcAccessor,
613  class DestIterator, class DestAccessor, class Array>
614 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
615  pair<DestIterator, DestAccessor> const & dest, bool background,
616  Array const & pixelPitch)
617 {
618  separableMultiDistance( source.first, source.second, source.third,
619  dest.first, dest.second, background, pixelPitch );
620 }
621 
622 template <class SrcIterator, class SrcShape, class SrcAccessor,
623  class DestIterator, class DestAccessor>
624 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
625  pair<DestIterator, DestAccessor> const & dest, bool background)
626 {
627  separableMultiDistance( source.first, source.second, source.third,
628  dest.first, dest.second, background );
629 }
630 
631 template <unsigned int N, class T1, class S1,
632  class T2, class S2, class Array>
633 inline void
634 separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
635  MultiArrayView<N, T2, S2> dest,
636  bool background,
637  Array const & pixelPitch)
638 {
639  vigra_precondition(source.shape() == dest.shape(),
640  "separableMultiDistance(): shape mismatch between input and output.");
641  separableMultiDistance( srcMultiArrayRange(source),
642  destMultiArray(dest), background, pixelPitch );
643 }
644 
645 template <unsigned int N, class T1, class S1,
646  class T2, class S2>
647 inline void
648 separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
649  MultiArrayView<N, T2, S2> dest,
650  bool background)
651 {
652  vigra_precondition(source.shape() == dest.shape(),
653  "separableMultiDistance(): shape mismatch between input and output.");
654  separableMultiDistance( srcMultiArrayRange(source),
655  destMultiArray(dest), background );
656 }
657 
658 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BoundaryDistanceTransform %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
659 
660 //rewrite labeled data and work with separableMultiDist
661 namespace lemon_graph {
662 
663 template <class Graph, class T1Map, class T2Map>
664 void
665 markRegionBoundaries(Graph const & g,
666  T1Map const & labels,
667  T2Map & out)
668 {
669  typedef typename Graph::NodeIt graph_scanner;
670  typedef typename Graph::OutBackArcIt neighbor_iterator;
671 
672  //find faces
673  for (graph_scanner node(g); node != INVALID; ++node)
674  {
675  typename T1Map::value_type center = labels[*node];
676 
677  for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
678  {
679  // set adjacent nodes with different labels to 1
680  if(center != labels[g.target(*arc)])
681  {
682  out[*node] = 1;
683  out[g.target(*arc)] = 1;
684  }
685  }
686  }
687 }
688 
689 } //-- namspace lemon_graph
690 
691 doxygen_overloaded_function(template <...> unsigned int markRegionBoundaries)
692 
693 template <unsigned int N, class T1, class S1,
694  class T2, class S2>
695 inline void
696 markRegionBoundaries(MultiArrayView<N, T1, S1> const & labels,
697  MultiArrayView<N, T2, S2> out,
698  NeighborhoodType neighborhood=DirectNeighborhood)
699 {
700  vigra_precondition(labels.shape() == out.shape(),
701  "markRegionBoundaries(): shape mismatch between input and output.");
702 
703  GridGraph<N, undirected_tag> graph(labels.shape(), neighborhood);
704 
705  lemon_graph::markRegionBoundaries(graph, labels, out);
706 }
707 
708 //MultiDistance which works directly on labeled data
709 
710 namespace detail
711 {
712 
713 /********************************************************/
714 /* */
715 /* boundaryDistParabola */
716 /* */
717 /********************************************************/
718 
719 template <class DestIterator, class LabelIterator>
720 void
721 boundaryDistParabola(DestIterator is, DestIterator iend,
722  LabelIterator ilabels,
723  double dmax,
724  bool array_border_is_active=false)
725 {
726  // We assume that the data in the input is distance squared and treat it as such
727  double w = iend - is;
728  if(w <= 0)
729  return;
730 
731  DestIterator id = is;
732  typedef typename LabelIterator::value_type LabelType;
733  typedef typename DestIterator::value_type DestType;
734  typedef detail::DistParabolaStackEntry<DestType> Influence;
735  typedef std::vector<Influence> Stack;
736 
737  double apex_height = array_border_is_active
738  ? 0.0
739  : dmax;
740  Stack _stack(1, Influence(apex_height, 0.0, -1.0, w));
741  LabelType current_label = *ilabels;
742  for(double begin = 0.0, current = 0.0; current <= w; ++ilabels, ++is, ++current)
743  {
744  apex_height = (current < w)
745  ? (current_label == *ilabels)
746  ? *is
747  : 0.0
748  : array_border_is_active
749  ? 0.0
750  : dmax;
751  while(true)
752  {
753  Influence & s = _stack.back();
754  double diff = current - s.center;
755  double intersection = current + (apex_height - s.apex_height - sq(diff)) / (2.0 * diff);
756 
757  if(intersection < s.left) // previous parabola has no influence
758  {
759  _stack.pop_back();
760  if(_stack.empty())
761  intersection = begin; // new parabola is valid for entire present segment
762  else
763  continue; // try new top of stack without advancing to next pixel
764  }
765  else if(intersection < s.right)
766  {
767  s.right = intersection;
768  }
769  if(intersection < w)
770  _stack.push_back(Influence(apex_height, intersection, current, w));
771  if(current < w && current_label == *ilabels)
772  break; // finished present pixel, advance to next one
773 
774  // label changed => finalize the current segment
775  typename Stack::iterator it = _stack.begin();
776  for(double c = begin; c < current; ++c, ++id)
777  {
778  while(c >= it->right)
779  ++it;
780  *id = sq(c - it->center) + it->apex_height;
781  }
782  if(current == w)
783  break; // stop when this was the last segment
784 
785  // initialize the new segment
786  begin = current;
787  current_label = *ilabels;
788  apex_height = *is;
789  Stack(1, Influence(0.0, begin-1.0, begin-1.0, w)).swap(_stack);
790  // don't advance to next pixel here, because the present pixel must also
791  // be analysed in the context of the new segment
792  }
793  }
794 }
795 
796 /********************************************************/
797 /* */
798 /* internalBoundaryMultiArrayDist */
799 /* */
800 /********************************************************/
801 
802 template <unsigned int N, class T1, class S1,
803  class T2, class S2>
804 void
805 internalBoundaryMultiArrayDist(
806  MultiArrayView<N, T1, S1> const & labels,
807  MultiArrayView<N, T2, S2> dest,
808  double dmax, bool array_border_is_active=false)
809 {
810  typedef typename MultiArrayView<N, T1, S1>::const_traverser LabelIterator;
811  typedef typename MultiArrayView<N, T2, S2>::traverser DestIterator;
812  typedef MultiArrayNavigator<LabelIterator, N> LabelNavigator;
813  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
814 
815  dest = dmax;
816  for( unsigned d = 0; d < N; ++d )
817  {
818  LabelNavigator lnav( labels.traverser_begin(), labels.shape(), d );
819  DNavigator dnav( dest.traverser_begin(), dest.shape(), d );
820 
821  for( ; dnav.hasMore(); dnav++, lnav++ )
822  {
823  boundaryDistParabola(dnav.begin(), dnav.end(),
824  lnav.begin(),
825  dmax, array_border_is_active);
826  }
827  }
828 }
829 
830 } // namespace detail
831 
832  /** \brief Specify which boundary is used for boundaryMultiDistance().
833 
834  */
835 enum BoundaryDistanceTag {
836  OuterBoundary, ///< Pixels just outside of each region
837  InterpixelBoundary, ///< Half-integer points between pixels of different labels
838  InnerBoundary ///< Pixels just inside of each region
839 };
840 
841 /********************************************************/
842 /* */
843 /* boundaryMultiDistance */
844 /* */
845 /********************************************************/
846 
847 /** \brief Euclidean distance to the implicit boundaries of a multi-dimensional label array.
848 
849 
850  <b> Declarations:</b>
851 
852  pass arbitrary-dimensional array views:
853  \code
854  namespace vigra {
855  template <unsigned int N, class T1, class S1,
856  class T2, class S2>
857  void
858  boundaryMultiDistance(MultiArrayView<N, T1, S1> const & labels,
859  MultiArrayView<N, T2, S2> dest,
860  bool array_border_is_active=false,
861  BoundaryDistanceTag boundary=InterpixelBoundary);
862  }
863  \endcode
864 
865  This function computes the distance transform of a labeled image <i>simultaneously</i>
866  for all regions. Depending on the requested type of \a boundary, three modes
867  are supported:
868  <ul>
869  <li><tt>OuterBoundary</tt>: In each region, compute the distance to the nearest pixel not
870  belonging to that regions. This is the same as if a normal distance transform
871  where applied to a binary image containing just this region.</li>
872  <li><tt>InterpixelBoundary</tt> (default): Like <tt>OuterBoundary</tt>, but shift the distance
873  to the interpixel boundary by subtractiong 1/2. This make the distences consistent
874  accross boundaries.</li>
875  <li><tt>InnerBoundary</tt>: In each region, compute the distance to the nearest pixel in the
876  region which is adjacent to the boundary. </li>
877  </ul>
878  If <tt>array_border_is_active=true</tt>, the
879  outer border of the array (i.e. the interpixel boundary between the array
880  and the infinite region) is also used. Otherwise (the default), regions
881  touching the array border are treated as if they extended to infinity.
882 
883  <b> Usage:</b>
884 
885  <b>\#include</b> <vigra/multi_distance.hxx><br/>
886  Namespace: vigra
887 
888  \code
889  Shape3 shape(width, height, depth);
890  MultiArray<3, unsigned char> source(shape);
891  MultiArray<3, UInt32> labels(shape);
892  MultiArray<3, float> dest(shape);
893  ...
894 
895  // find regions (interpixel boundaries are implied)
896  labelMultiArray(source, labels);
897 
898  // Calculate Euclidean distance to interpixel boundary for all pixels
899  boundaryMultiDistance(labels, dest);
900  \endcode
901 
902  \see vigra::distanceTransform(), vigra::separableMultiDistance()
903 */
904 doxygen_overloaded_function(template <...> void boundaryMultiDistance)
905 
906 template <unsigned int N, class T1, class S1,
907  class T2, class S2>
908 void
909 boundaryMultiDistance(MultiArrayView<N, T1, S1> const & labels,
910  MultiArrayView<N, T2, S2> dest,
911  bool array_border_is_active=false,
912  BoundaryDistanceTag boundary=InterpixelBoundary)
913 {
914  vigra_precondition(labels.shape() == dest.shape(),
915  "boundaryMultiDistance(): shape mismatch between input and output.");
916 
917  using namespace vigra::functor;
918 
919  if(boundary == InnerBoundary)
920  {
921  MultiArray<N, unsigned char> boundaries(labels.shape());
922 
923  markRegionBoundaries(labels, boundaries, IndirectNeighborhood);
924  if(array_border_is_active)
925  initMultiArrayBorder(boundaries, 1, 1);
926  separableMultiDistance(boundaries, dest, true);
927  }
928  else
929  {
930  T2 offset = 0.0;
931 
932  if(boundary == InterpixelBoundary)
933  {
934  vigra_precondition(!NumericTraits<T2>::isIntegral::value,
935  "boundaryMultiDistance(..., InterpixelBoundary): output pixel type must be float or double.");
936  offset = T2(0.5);
937  }
938  double dmax = squaredNorm(labels.shape()) + N;
939  if(dmax > double(NumericTraits<T2>::max()))
940  {
941  // need a temporary array to avoid overflows
942  typedef typename NumericTraits<T2>::RealPromote Real;
943  MultiArray<N, Real> tmpArray(labels.shape());
944  detail::internalBoundaryMultiArrayDist(labels, tmpArray,
945  dmax, array_border_is_active);
946  transformMultiArray(tmpArray, dest, sqrt(Arg1()) - Param(offset) );
947  }
948  else
949  {
950  // can work directly on the destination array
951  detail::internalBoundaryMultiArrayDist(labels, dest, dmax, array_border_is_active);
952  transformMultiArray(dest, dest, sqrt(Arg1()) - Param(offset) );
953  }
954  }
955 }
956 
957 //@}
958 
959 } //-- namespace vigra
960 
961 
962 #endif //-- VIGRA_MULTI_DISTANCE_HXX
vigra::initMultiArrayBorder
void initMultiArrayBorder(...)
Write values to the specified border values in the array.
vigra::boundaryMultiDistance
void boundaryMultiDistance(...)
Euclidean distance to the implicit boundaries of a multi-dimensional label array. ...
vigra::BoundaryDistanceTag
BoundaryDistanceTag
Specify which boundary is used for boundaryMultiDistance().
Definition: multi_distance.hxx:835
vigra::squaredNorm
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
vigra::separableMultiDistance
void separableMultiDistance(...)
Euclidean distance on multi-dimensional arrays.
vigra::InnerBoundary
Pixels just inside of each region.
Definition: multi_distance.hxx:838
vigra::OuterBoundary
Pixels just outside of each region.
Definition: multi_distance.hxx:836
vigra::sq
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
vigra::boost_graph::source
vigra::GridGraph< N, DirectedTag >::vertex_descriptor source(typename vigra::GridGraph< N, DirectedTag >::edge_descriptor const &e, vigra::GridGraph< N, DirectedTag > const &g)
Get a vertex descriptor for the start vertex of edge e in graph g (API: boost).
Definition: multi_gridgraph.hxx:2943
vigra::doxygen_overloaded_function
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::copyMultiArray
void copyMultiArray(...)
Copy a multi-dimensional array.
vigra::IndirectNeighborhood
use direct and indirect neighbors
Definition: multi_fwd.hxx:188
vigra::InterpixelBoundary
Half-integer points between pixels of different labels.
Definition: multi_distance.hxx:837
vigra::separableMultiDistSquared
void separableMultiDistSquared(...)
Euclidean distance squared on multi-dimensional arrays.
vigra::transformMultiArray
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
vigra::ceil
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
vigra::DirectNeighborhood
use only direct neighbors
Definition: multi_fwd.hxx:187
vigra::NeighborhoodType
NeighborhoodType
Choose the neighborhood system in a dimension-independent way.
Definition: multi_fwd.hxx:186
vigra::sqrt
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616

© 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.11.1 (Fri May 19 2017)