[ 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 MultiArrayDistanceTransform Euclidean distance transform for multi-dimensional arrays.
239 
240  These functions perform variants of the Euclidean distance transform on
241  arbitrary dimensional arrays.
242 */
243 //@{
244 
245 /********************************************************/
246 /* */
247 /* separableMultiDistSquared */
248 /* */
249 /********************************************************/
250 
251 /** \brief Euclidean distance squared on multi-dimensional arrays.
252 
253  The algorithm is taken from Donald Bailey: "An Efficient Euclidean Distance Transform",
254  Proc. IWCIA'04, Springer LNCS 3322, 2004.
255 
256  <b> Declarations:</b>
257 
258  pass arbitrary-dimensional array views:
259  \code
260  namespace vigra {
261  // explicitly specify pixel pitch for each coordinate
262  template <unsigned int N, class T1, class S1,
263  class T2, class S2,
264  class Array>
265  void
266  separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
267  MultiArrayView<N, T2, S2> dest,
268  bool background,
269  Array const & pixelPitch);
270 
271  // use default pixel pitch = 1.0 for each coordinate
272  template <unsigned int N, class T1, class S1,
273  class T2, class S2>
274  void
275  separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
276  MultiArrayView<N, T2, S2> dest,
277  bool background);
278  }
279  \endcode
280 
281  \deprecatedAPI{separableMultiDistSquared}
282  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
283  \code
284  namespace vigra {
285  // explicitly specify pixel pitch for each coordinate
286  template <class SrcIterator, class SrcShape, class SrcAccessor,
287  class DestIterator, class DestAccessor, class Array>
288  void
289  separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
290  DestIterator d, DestAccessor dest,
291  bool background,
292  Array const & pixelPitch);
293 
294  // use default pixel pitch = 1.0 for each coordinate
295  template <class SrcIterator, class SrcShape, class SrcAccessor,
296  class DestIterator, class DestAccessor>
297  void
298  separableMultiDistSquared(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
299  DestIterator diter, DestAccessor dest,
300  bool background);
301 
302  }
303  \endcode
304  use argument objects in conjunction with \ref ArgumentObjectFactories :
305  \code
306  namespace vigra {
307  // explicitly specify pixel pitch for each coordinate
308  template <class SrcIterator, class SrcShape, class SrcAccessor,
309  class DestIterator, class DestAccessor, class Array>
310  void
311  separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
312  pair<DestIterator, DestAccessor> const & dest,
313  bool background,
314  Array const & pixelPitch);
315 
316  // use default pixel pitch = 1.0 for each coordinate
317  template <class SrcIterator, class SrcShape, class SrcAccessor,
318  class DestIterator, class DestAccessor>
319  void
320  separableMultiDistSquared(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
321  pair<DestIterator, DestAccessor> const & dest,
322  bool background);
323 
324  }
325  \endcode
326  \deprecatedEnd
327 
328  This function performs a squared Euclidean squared distance transform on the given
329  multi-dimensional array. Both source and destination
330  arrays are represented by iterators, shape objects and accessors.
331  The destination array is required to already have the correct size.
332 
333  This function expects a mask as its source, where background pixels are
334  marked as zero, and non-background pixels as non-zero. If the parameter
335  <i>background</i> is true, then the squared distance of all background
336  pixels to the nearest object is calculated. Otherwise, the distance of all
337  object pixels to the nearest background pixel is calculated.
338 
339  Optionally, one can pass an array that specifies the pixel pitch in each direction.
340  This is necessary when the data have non-uniform resolution (as is common in confocal
341  microscopy, for example).
342 
343  This function may work in-place, which means that <tt>siter == diter</tt> is allowed.
344  A full-sized internal array is only allocated if working on the destination
345  array directly would cause overflow errors (i.e. if
346  <tt> NumericTraits<typename DestAccessor::value_type>::max() < N * M*M</tt>, where M is the
347  size of the largest dimension of the array.
348 
349  <b> Usage:</b>
350 
351  <b>\#include</b> <vigra/multi_distance.hxx><br/>
352  Namespace: vigra
353 
354  \code
355  Shape3 shape(width, height, depth);
356  MultiArray<3, unsigned char> source(shape);
357  MultiArray<3, unsigned int> dest(shape);
358  ...
359 
360  // Calculate Euclidean distance squared for all background pixels
361  separableMultiDistSquared(source, dest, true);
362  \endcode
363 
364  \see vigra::distanceTransform(), vigra::separableMultiDistance()
365 */
367 
368 template <class SrcIterator, class SrcShape, class SrcAccessor,
369  class DestIterator, class DestAccessor, class Array>
370 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
371  DestIterator d, DestAccessor dest, bool background,
372  Array const & pixelPitch)
373 {
374  int N = shape.size();
375 
376  typedef typename SrcAccessor::value_type SrcType;
377  typedef typename DestAccessor::value_type DestType;
378  typedef typename NumericTraits<DestType>::RealPromote Real;
379 
380  SrcType zero = NumericTraits<SrcType>::zero();
381 
382  double dmax = 0.0;
383  bool pixelPitchIsReal = false;
384  for( int k=0; k<N; ++k)
385  {
386  if(int(pixelPitch[k]) != pixelPitch[k])
387  pixelPitchIsReal = true;
388  dmax += sq(pixelPitch[k]*shape[k]);
389  }
390 
391  using namespace vigra::functor;
392 
393  if(dmax > NumericTraits<DestType>::toRealPromote(NumericTraits<DestType>::max())
394  || pixelPitchIsReal) // need a temporary array to avoid overflows
395  {
396  // Threshold the values so all objects have infinity value in the beginning
397  Real maxDist = (Real)dmax, rzero = (Real)0.0;
398  MultiArray<SrcShape::static_size, Real> tmpArray(shape);
399  if(background == true)
400  transformMultiArray( s, shape, src,
401  tmpArray.traverser_begin(), typename AccessorTraits<Real>::default_accessor(),
402  ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
403  else
404  transformMultiArray( s, shape, src,
405  tmpArray.traverser_begin(), typename AccessorTraits<Real>::default_accessor(),
406  ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
407 
408  detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(),
409  shape, typename AccessorTraits<Real>::default_accessor(),
410  tmpArray.traverser_begin(),
411  typename AccessorTraits<Real>::default_accessor(), pixelPitch);
412 
413  copyMultiArray(srcMultiArrayRange(tmpArray), destIter(d, dest));
414  }
415  else // work directly on the destination array
416  {
417  // Threshold the values so all objects have infinity value in the beginning
418  DestType maxDist = DestType(std::ceil(dmax)), rzero = (DestType)0;
419  if(background == true)
420  transformMultiArray( s, shape, src, d, dest,
421  ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
422  else
423  transformMultiArray( s, shape, src, d, dest,
424  ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
425 
426  detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest, pixelPitch);
427  }
428 }
429 
430 template <class SrcIterator, class SrcShape, class SrcAccessor,
431  class DestIterator, class DestAccessor>
432 inline
433 void separableMultiDistSquared( SrcIterator s, SrcShape const & shape, SrcAccessor src,
434  DestIterator d, DestAccessor dest, bool background)
435 {
436  ArrayVector<double> pixelPitch(shape.size(), 1.0);
437  separableMultiDistSquared( s, shape, src, d, dest, background, pixelPitch );
438 }
439 
440 template <class SrcIterator, class SrcShape, class SrcAccessor,
441  class DestIterator, class DestAccessor, class Array>
442 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
443  pair<DestIterator, DestAccessor> const & dest, bool background,
444  Array const & pixelPitch)
445 {
446  separableMultiDistSquared( source.first, source.second, source.third,
447  dest.first, dest.second, background, pixelPitch );
448 }
449 
450 template <class SrcIterator, class SrcShape, class SrcAccessor,
451  class DestIterator, class DestAccessor>
452 inline void separableMultiDistSquared( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
453  pair<DestIterator, DestAccessor> const & dest, bool background)
454 {
455  separableMultiDistSquared( source.first, source.second, source.third,
456  dest.first, dest.second, background );
457 }
458 
459 template <unsigned int N, class T1, class S1,
460  class T2, class S2,
461  class Array>
462 inline void
463 separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
464  MultiArrayView<N, T2, S2> dest, bool background,
465  Array const & pixelPitch)
466 {
467  vigra_precondition(source.shape() == dest.shape(),
468  "separableMultiDistSquared(): shape mismatch between input and output.");
469  separableMultiDistSquared( srcMultiArrayRange(source),
470  destMultiArray(dest), background, pixelPitch );
471 }
472 
473 template <unsigned int N, class T1, class S1,
474  class T2, class S2>
475 inline void
476 separableMultiDistSquared(MultiArrayView<N, T1, S1> const & source,
477  MultiArrayView<N, T2, S2> dest, bool background)
478 {
479  vigra_precondition(source.shape() == dest.shape(),
480  "separableMultiDistSquared(): shape mismatch between input and output.");
481  separableMultiDistSquared( srcMultiArrayRange(source),
482  destMultiArray(dest), background );
483 }
484 
485 /********************************************************/
486 /* */
487 /* separableMultiDistance */
488 /* */
489 /********************************************************/
490 
491 /** \brief Euclidean distance on multi-dimensional arrays.
492 
493  <b> Declarations:</b>
494 
495  pass arbitrary-dimensional array views:
496  \code
497  namespace vigra {
498  // explicitly specify pixel pitch for each coordinate
499  template <unsigned int N, class T1, class S1,
500  class T2, class S2, class Array>
501  void
502  separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
503  MultiArrayView<N, T2, S2> dest,
504  bool background,
505  Array const & pixelPitch);
506 
507  // use default pixel pitch = 1.0 for each coordinate
508  template <unsigned int N, class T1, class S1,
509  class T2, class S2>
510  void
511  separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
512  MultiArrayView<N, T2, S2> dest,
513  bool background);
514  }
515  \endcode
516 
517  \deprecatedAPI{separableMultiDistance}
518  pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
519  \code
520  namespace vigra {
521  // explicitly specify pixel pitch for each coordinate
522  template <class SrcIterator, class SrcShape, class SrcAccessor,
523  class DestIterator, class DestAccessor, class Array>
524  void
525  separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
526  DestIterator d, DestAccessor dest,
527  bool background,
528  Array const & pixelPitch);
529 
530  // use default pixel pitch = 1.0 for each coordinate
531  template <class SrcIterator, class SrcShape, class SrcAccessor,
532  class DestIterator, class DestAccessor>
533  void
534  separableMultiDistance(SrcIterator siter, SrcShape const & shape, SrcAccessor src,
535  DestIterator diter, DestAccessor dest,
536  bool background);
537 
538  }
539  \endcode
540  use argument objects in conjunction with \ref ArgumentObjectFactories :
541  \code
542  namespace vigra {
543  // explicitly specify pixel pitch for each coordinate
544  template <class SrcIterator, class SrcShape, class SrcAccessor,
545  class DestIterator, class DestAccessor, class Array>
546  void
547  separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
548  pair<DestIterator, DestAccessor> const & dest,
549  bool background,
550  Array const & pixelPitch);
551 
552  // use default pixel pitch = 1.0 for each coordinate
553  template <class SrcIterator, class SrcShape, class SrcAccessor,
554  class DestIterator, class DestAccessor>
555  void
556  separableMultiDistance(triple<SrcIterator, SrcShape, SrcAccessor> const & source,
557  pair<DestIterator, DestAccessor> const & dest,
558  bool background);
559 
560  }
561  \endcode
562  \deprecatedEnd
563 
564  This function performs a Euclidean distance transform on the given
565  multi-dimensional array. It simply calls \ref separableMultiDistSquared()
566  and takes the pixel-wise square root of the result. See \ref separableMultiDistSquared()
567  for more documentation.
568 
569  <b> Usage:</b>
570 
571  <b>\#include</b> <vigra/multi_distance.hxx><br/>
572  Namespace: vigra
573 
574  \code
575  Shape3 shape(width, height, depth);
576  MultiArray<3, unsigned char> source(shape);
577  MultiArray<3, float> dest(shape);
578  ...
579 
580  // Calculate Euclidean distance for all background pixels
581  separableMultiDistance(source, dest, true);
582  \endcode
583 
584  \see vigra::distanceTransform(), vigra::separableMultiDistSquared()
585 */
587 
588 template <class SrcIterator, class SrcShape, class SrcAccessor,
589  class DestIterator, class DestAccessor, class Array>
590 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
591  DestIterator d, DestAccessor dest, bool background,
592  Array const & pixelPitch)
593 {
594  separableMultiDistSquared( s, shape, src, d, dest, background, pixelPitch);
595 
596  // Finally, calculate the square root of the distances
597  using namespace vigra::functor;
598 
599  transformMultiArray( d, shape, dest, d, dest, sqrt(Arg1()) );
600 }
601 
602 template <class SrcIterator, class SrcShape, class SrcAccessor,
603  class DestIterator, class DestAccessor>
604 void separableMultiDistance( SrcIterator s, SrcShape const & shape, SrcAccessor src,
605  DestIterator d, DestAccessor dest, bool background)
606 {
607  separableMultiDistSquared( s, shape, src, d, dest, background);
608 
609  // Finally, calculate the square root of the distances
610  using namespace vigra::functor;
611 
612  transformMultiArray( d, shape, dest, d, dest, sqrt(Arg1()) );
613 }
614 
615 template <class SrcIterator, class SrcShape, class SrcAccessor,
616  class DestIterator, class DestAccessor, class Array>
617 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
618  pair<DestIterator, DestAccessor> const & dest, bool background,
619  Array const & pixelPitch)
620 {
621  separableMultiDistance( source.first, source.second, source.third,
622  dest.first, dest.second, background, pixelPitch );
623 }
624 
625 template <class SrcIterator, class SrcShape, class SrcAccessor,
626  class DestIterator, class DestAccessor>
627 inline void separableMultiDistance( triple<SrcIterator, SrcShape, SrcAccessor> const & source,
628  pair<DestIterator, DestAccessor> const & dest, bool background)
629 {
630  separableMultiDistance( source.first, source.second, source.third,
631  dest.first, dest.second, background );
632 }
633 
634 template <unsigned int N, class T1, class S1,
635  class T2, class S2, class Array>
636 inline void
637 separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
638  MultiArrayView<N, T2, S2> dest,
639  bool background,
640  Array const & pixelPitch)
641 {
642  vigra_precondition(source.shape() == dest.shape(),
643  "separableMultiDistance(): shape mismatch between input and output.");
644  separableMultiDistance( srcMultiArrayRange(source),
645  destMultiArray(dest), background, pixelPitch );
646 }
647 
648 template <unsigned int N, class T1, class S1,
649  class T2, class S2>
650 inline void
651 separableMultiDistance(MultiArrayView<N, T1, S1> const & source,
652  MultiArrayView<N, T2, S2> dest,
653  bool background)
654 {
655  vigra_precondition(source.shape() == dest.shape(),
656  "separableMultiDistance(): shape mismatch between input and output.");
657  separableMultiDistance( srcMultiArrayRange(source),
658  destMultiArray(dest), background );
659 }
660 
661 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BoundaryDistanceTransform %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662 
663 //rewrite labeled data and work with separableMultiDist
664 namespace lemon_graph {
665 
666 template <class Graph, class T1Map, class T2Map>
667 void
668 markRegionBoundaries(Graph const & g,
669  T1Map const & labels,
670  T2Map & out)
671 {
672  typedef typename Graph::NodeIt graph_scanner;
673  typedef typename Graph::OutBackArcIt neighbor_iterator;
674 
675  //find faces
676  for (graph_scanner node(g); node != INVALID; ++node)
677  {
678  typename T1Map::value_type center = labels[*node];
679 
680  for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
681  {
682  // set adjacent nodes with different labels to 1
683  if(center != labels[g.target(*arc)])
684  {
685  out[*node] = 1;
686  out[g.target(*arc)] = 1;
687  }
688  }
689  }
690 }
691 
692 } //-- namspace lemon_graph
693 
694 doxygen_overloaded_function(template <...> unsigned int markRegionBoundaries)
695 
696 template <unsigned int N, class T1, class S1,
697  class T2, class S2>
698 inline void
699 markRegionBoundaries(MultiArrayView<N, T1, S1> const & labels,
700  MultiArrayView<N, T2, S2> out,
702 {
703  vigra_precondition(labels.shape() == out.shape(),
704  "markRegionBoundaries(): shape mismatch between input and output.");
705 
706  GridGraph<N, undirected_tag> graph(labels.shape(), neighborhood);
707 
708  lemon_graph::markRegionBoundaries(graph, labels, out);
709 }
710 
711 //MultiDistance which works directly on labeled data
712 
713 namespace detail
714 {
715 
716 /********************************************************/
717 /* */
718 /* boundaryDistParabola */
719 /* */
720 /********************************************************/
721 
722 template <class DestIterator, class LabelIterator>
723 void
724 boundaryDistParabola(DestIterator is, DestIterator iend,
725  LabelIterator ilabels,
726  double dmax,
727  bool array_border_is_active=false)
728 {
729  // We assume that the data in the input is distance squared and treat it as such
730  double w = iend - is;
731  if(w <= 0)
732  return;
733 
734  DestIterator id = is;
735  typedef typename LabelIterator::value_type LabelType;
736  typedef typename DestIterator::value_type DestType;
737  typedef detail::DistParabolaStackEntry<DestType> Influence;
738  typedef std::vector<Influence> Stack;
739 
740  double apex_height = array_border_is_active
741  ? 0.0
742  : dmax;
743  Stack _stack(1, Influence(apex_height, 0.0, -1.0, w));
744  LabelType current_label = *ilabels;
745  for(double begin = 0.0, current = 0.0; current <= w; ++ilabels, ++is, ++current)
746  {
747  apex_height = (current < w)
748  ? (current_label == *ilabels)
749  ? *is
750  : 0.0
751  : array_border_is_active
752  ? 0.0
753  : dmax;
754  while(true)
755  {
756  Influence & s = _stack.back();
757  double diff = current - s.center;
758  double intersection = current + (apex_height - s.apex_height - sq(diff)) / (2.0 * diff);
759 
760  if(intersection < s.left) // previous parabola has no influence
761  {
762  _stack.pop_back();
763  if(_stack.empty())
764  intersection = begin; // new parabola is valid for entire present segment
765  else
766  continue; // try new top of stack without advancing to next pixel
767  }
768  else if(intersection < s.right)
769  {
770  s.right = intersection;
771  }
772  if(intersection < w)
773  _stack.push_back(Influence(apex_height, intersection, current, w));
774  if(current < w && current_label == *ilabels)
775  break; // finished present pixel, advance to next one
776 
777  // label changed => finalize the current segment
778  typename Stack::iterator it = _stack.begin();
779  for(double c = begin; c < current; ++c, ++id)
780  {
781  while(c >= it->right)
782  ++it;
783  *id = sq(c - it->center) + it->apex_height;
784  }
785  if(current == w)
786  break; // stop when this was the last segment
787 
788  // initialize the new segment
789  begin = current;
790  current_label = *ilabels;
791  apex_height = *is;
792  Stack(1, Influence(0.0, begin-1.0, begin-1.0, w)).swap(_stack);
793  // don't advance to next pixel here, because the present pixel must also
794  // be analysed in the context of the new segment
795  }
796  }
797 }
798 
799 /********************************************************/
800 /* */
801 /* internalBoundaryMultiArrayDist */
802 /* */
803 /********************************************************/
804 
805 template <unsigned int N, class T1, class S1,
806  class T2, class S2>
807 void
808 internalBoundaryMultiArrayDist(
809  MultiArrayView<N, T1, S1> const & labels,
810  MultiArrayView<N, T2, S2> dest,
811  double dmax, bool array_border_is_active=false)
812 {
813  typedef typename MultiArrayView<N, T1, S1>::const_traverser LabelIterator;
814  typedef typename MultiArrayView<N, T2, S2>::traverser DestIterator;
815  typedef MultiArrayNavigator<LabelIterator, N> LabelNavigator;
816  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
817 
818  dest = dmax;
819  for( int d = 0; d < N; ++d )
820  {
821  LabelNavigator lnav( labels.traverser_begin(), labels.shape(), d );
822  DNavigator dnav( dest.traverser_begin(), dest.shape(), d );
823 
824  for( ; dnav.hasMore(); dnav++, lnav++ )
825  {
826  boundaryDistParabola(dnav.begin(), dnav.end(),
827  lnav.begin(),
828  dmax, array_border_is_active);
829  }
830  }
831 }
832 
833 } // namespace detail
834 
835  /** \brief Specify which boundary is used for boundaryMultiDistance().
836 
837  */
839  OuterBoundary, ///< Pixels just outside of each region
840  InterpixelBoundary, ///< Half-integer points between pixels of different labels
841  InnerBoundary ///< Pixels just inside of each region
842 };
843 
844 /********************************************************/
845 /* */
846 /* boundaryMultiDistance */
847 /* */
848 /********************************************************/
849 
850 /** \brief Euclidean distance to the implicit boundaries of a multi-dimensional label array.
851 
852 
853  <b> Declarations:</b>
854 
855  pass arbitrary-dimensional array views:
856  \code
857  namespace vigra {
858  template <unsigned int N, class T1, class S1,
859  class T2, class S2>
860  void
861  boundaryMultiDistance(MultiArrayView<N, T1, S1> const & labels,
862  MultiArrayView<N, T2, S2> dest,
863  bool array_border_is_active=false,
864  BoundaryDistanceTag boundary=InterpixelBoundary);
865  }
866  \endcode
867 
868  This function computes the distance transform of a labeled image <i>simultaneously</i>
869  for all regions. Depending on the requested type of \a boundary, three modes
870  are supported:
871  <ul>
872  <li><tt>OuterBoundary</tt>: In each region, compute the distance to the nearest pixel not
873  belonging to that regions. This is the same as if a normal distance transform
874  where applied to a binary image containing just this region.</li>
875  <li><tt>InterpixelBoundary</tt> (default): Like <tt>OuterBoundary</tt>, but shift the distance
876  to the interpixel boundary by subtractiong 1/2. This make the distences consistent
877  accross boundaries.</li>
878  <li><tt>InnerBoundary</tt>: In each region, compute the distance to the nearest pixel in the
879  region which is adjacent to the boundary. </li>
880  </ul>
881  If <tt>array_border_is_active=true</tt>, the
882  outer border of the array (i.e. the interpixel boundary between the array
883  and the infinite region) is also used. Otherwise (the default), regions
884  touching the array border are treated as if they extended to infinity.
885 
886  <b> Usage:</b>
887 
888  <b>\#include</b> <vigra/multi_distance.hxx><br/>
889  Namespace: vigra
890 
891  \code
892  Shape3 shape(width, height, depth);
893  MultiArray<3, unsigned char> source(shape);
894  MultiArray<3, UInt32> labels(shape);
895  MultiArray<3, float> dest(shape);
896  ...
897 
898  // find regions (interpixel boundaries are implied)
899  labelMultiArray(source, labels);
900 
901  // Calculate Euclidean distance to interpixel boundary for all pixels
902  boundaryMultiDistance(labels, dest);
903  \endcode
904 
905  \see vigra::distanceTransform(), vigra::separableMultiDistance()
906 */
908 
909 template <unsigned int N, class T1, class S1,
910  class T2, class S2>
911 void
914  bool array_border_is_active=false,
916 {
917  vigra_precondition(labels.shape() == dest.shape(),
918  "boundaryMultiDistance(): shape mismatch between input and output.");
919 
920  using namespace vigra::functor;
921 
922  if(boundary == InnerBoundary)
923  {
924  MultiArray<N, unsigned char> boundaries(labels.shape());
925 
926  markRegionBoundaries(labels, boundaries, IndirectNeighborhood);
927  if(array_border_is_active)
928  initMultiArrayBorder(boundaries, 1, 1);
929  separableMultiDistance(boundaries, dest, true);
930  }
931  else
932  {
933  T2 offset = 0.0;
934 
935  if(boundary == InterpixelBoundary)
936  {
937  vigra_precondition(!NumericTraits<T2>::isIntegral::value,
938  "boundaryMultiDistance(..., InterpixelBoundary): output pixel type must be float or double.");
939  offset = T2(0.5);
940  }
941  double dmax = squaredNorm(labels.shape()) + N;
942  if(dmax > double(NumericTraits<T2>::max()))
943  {
944  // need a temporary array to avoid overflows
945  typedef typename NumericTraits<T2>::RealPromote Real;
946  MultiArray<N, Real> tmpArray(labels.shape());
947  detail::internalBoundaryMultiArrayDist(labels, tmpArray,
948  dmax, array_border_is_active);
949  transformMultiArray(tmpArray, dest, sqrt(Arg1()) - Param(offset) );
950  }
951  else
952  {
953  // can work directly on the destination array
954  detail::internalBoundaryMultiArrayDist(labels, dest, dmax, array_border_is_active);
955  transformMultiArray(dest, dest, sqrt(Arg1()) - Param(offset) );
956  }
957  }
958 }
959 
960 //@}
961 
962 } //-- namespace vigra
963 
964 
965 #endif //-- VIGRA_MULTI_DISTANCE_HXX
void boundaryMultiDistance(...)
Euclidean distance to the implicit boundaries of a multi-dimensional label array. ...
void initMultiArrayBorder(...)
Write values to the specified border values in the array.
BoundaryDistanceTag
Specify which boundary is used for boundaryMultiDistance().
Definition: multi_distance.hxx:838
const difference_type & shape() const
Definition: multi_array.hxx:1596
Pixels just inside of each region.
Definition: multi_distance.hxx:841
void separableMultiDistSquared(...)
Euclidean distance squared on multi-dimensional arrays.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2422
Definition: accessor.hxx:43
Pixels just outside of each region.
Definition: multi_distance.hxx:839
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, T const &, T const * >::type const_traverser
Definition: multi_array.hxx:717
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, T &, T * >::type traverser
Definition: multi_array.hxx:712
void separableMultiDistance(...)
Euclidean distance on multi-dimensional arrays.
Definition: inspectimage.hxx:255
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:365
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
void copyMultiArray(...)
Copy a multi-dimensional array.
use direct and indirect neighbors
Definition: multi_fwd.hxx:188
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:652
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
Half-integer points between pixels of different labels.
Definition: multi_distance.hxx:840
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
use only direct neighbors
Definition: multi_fwd.hxx:187
NeighborhoodType
Choose the neighborhood system in a dimension-independent way.
Definition: multi_fwd.hxx:186
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.10.0