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

vector_distance.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2013 by Thorben Kroeger, Kasim Terzic, */
4 /* Christian-Dennis Rahn 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_VECTOR_DISTANCE_HXX
38 #define VIGRA_VECTOR_DISTANCE_HXX
39 
40 #include <vector>
41 #include <set>
42 #include <functional>
43 #include "array_vector.hxx"
44 #include "multi_array.hxx"
45 #include "accessor.hxx"
46 #include "numerictraits.hxx"
47 #include "navigator.hxx"
48 #include "metaprogramming.hxx"
49 #include "multi_pointoperators.hxx"
50 #include "functorexpression.hxx"
51 #include "multi_distance.hxx"
52 
53 #undef VECTORIAL_DIST_DEBUG
54 
55 namespace vigra
56 {
57 
58 namespace detail
59 {
60 
61 template <class Vector, class Value>
62 struct VectorialDistParabolaStackEntry
63 {
64  double left, center, right;
65  Value apex_height;
66  Vector point;
67 
68  VectorialDistParabolaStackEntry(const Vector& vec, Value prev, double l, double c, double r)
69  : left(l), center(c), right(r), apex_height(prev), point(vec)
70  {}
71 };
72 
73 #ifdef VECTORIAL_DIST_DEBUG
74 template <class Vector, class Value>
75 std::ostream& operator<<(std::ostream&o, const VectorialDistParabolaStackEntry<Vector, Value>& e) {
76  o << "l=" << e.left << ", c=" << e.center << ", r=" << e.right << ", pV=" << e.apex_height << ", pVec=" << e.point;
77  return o;
78 }
79 #endif
80 
81 /********************************************************/
82 /* */
83 /* vectorialDistParabola */
84 /* */
85 /********************************************************/
86 
87 template <class VEC, class ARRAY>
88 inline double
89 partialSquaredMagnitude(const VEC& vec, MultiArrayIndex dim, ARRAY const & pixel_pitch)
90 {
91  //computes the squared magnitude of vec
92  //considering only the first dim dimensions
93  double sqMag = 0.0;
94  for(MultiArrayIndex i=0; i<=dim; ++i)
95  {
96  sqMag += sq(pixel_pitch[i]*vec[i]);
97  }
98  return sqMag;
99 }
100 
101 template <class SrcIterator,
102  class Array>
103 void
104 vectorialDistParabola(MultiArrayIndex dimension,
105  SrcIterator is, SrcIterator iend,
106  Array const & pixel_pitch )
107 {
108  typedef typename SrcIterator::value_type SrcType;
109  typedef VectorialDistParabolaStackEntry<SrcType, double> Influence;
110 
111  double sigma = pixel_pitch[dimension],
112  sigma2 = sq(sigma);
113  double w = iend - is; //width of the scanline
114 
115  SrcIterator id = is;
116 
117  std::vector<Influence> _stack; //stack of influence parabolas
118  double apex_height = partialSquaredMagnitude(*is, dimension, pixel_pitch);
119  _stack.push_back(Influence(*is, apex_height, 0.0, 0.0, w));
120  ++is;
121  double current = 1.0;
122  while(current < w)
123  {
124  apex_height = partialSquaredMagnitude(*is, dimension, pixel_pitch);
125  Influence & s = _stack.back();
126  double diff = current - s.center;
127  double intersection = current + (apex_height - s.apex_height - sq(sigma*diff)) / (2.0*sigma2 * diff);
128 
129  if( intersection < s.left) // previous point has no influence
130  {
131  _stack.pop_back();
132  if(_stack.empty())
133  _stack.push_back(Influence(*is, apex_height, 0.0, current, w));
134  else
135  continue; // try new top of stack without advancing current
136  }
137  else if(intersection < s.right)
138  {
139  s.right = intersection;
140  _stack.push_back(Influence(*is, apex_height, intersection, current, w));
141  }
142  ++is;
143  ++current;
144  }
145 
146  // Now we have the stack indicating which rows are influenced by (and therefore
147  // closest to) which row. We can go through the stack and calculate the
148  // distance squared for each element of the column.
149  typename std::vector<Influence>::iterator it = _stack.begin();
150  for(current = 0.0; current < w; ++current, ++id)
151  {
152  while( current >= it->right)
153  ++it;
154 
155  *id = it->point;
156  (*id)[dimension] = it->center - current;
157  }
158 }
159 
160 template <class DestIterator,
161  class LabelIterator,
162  class Array1, class Array2>
163 void
164 boundaryVectorDistParabola(MultiArrayIndex dimension,
165  DestIterator is, DestIterator iend,
166  LabelIterator ilabels,
167  Array1 const & pixel_pitch,
168  Array2 const & dmax,
169  bool array_border_is_active=false)
170 {
171  double w = iend - is;
172  if(w <= 0)
173  return;
174 
175  typedef typename LabelIterator::value_type LabelType;
176  typedef typename DestIterator::value_type DestType;
177  typedef VectorialDistParabolaStackEntry<DestType, double> Influence;
178  typedef std::vector<Influence> Stack;
179 
180  DestIterator id = is;
181  DestType border_point = array_border_is_active
182  ? DestType(0)
183  : dmax;
184  double apex_height = partialSquaredMagnitude(border_point, dimension, pixel_pitch);
185  Stack _stack(1, Influence(border_point, apex_height, 0.0, -1.0, w));
186  LabelType current_label = *ilabels;
187  for(double begin = 0.0, current = 0.0; current <= w; ++ilabels, ++is, ++current)
188  {
189  DestType point = (current < w)
190  ? (current_label == *ilabels)
191  ? *is
192  : DestType(0)
193  : border_point;
194  apex_height = partialSquaredMagnitude(point, dimension, pixel_pitch);
195  while(true)
196  {
197  Influence & s = _stack.back();
198  double diff = (current - s.center)*pixel_pitch[dimension];
199  double intersection = current + (apex_height - s.apex_height - sq(diff)) / (2.0 * diff);
200 
201  if(intersection < s.left) // previous parabola has no influence
202  {
203  _stack.pop_back();
204  if(_stack.empty())
205  intersection = begin; // new parabola is valid for entire present segment
206  else
207  continue; // try new top of stack without advancing to next pixel
208  }
209  else if(intersection < s.right)
210  {
211  s.right = intersection;
212  }
213  if(intersection < w)
214  _stack.push_back(Influence(point, apex_height, intersection, current, w));
215  if(current < w && current_label == *ilabels)
216  break; // finished present pixel, advance to next one
217 
218  // label changed => finalize the current segment
219  typename Stack::iterator it = _stack.begin();
220  for(double c = begin; c < current; ++c, ++id)
221  {
222  while(c >= it->right)
223  ++it;
224  *id = it->point;
225  (*id)[dimension] = it->center - c;
226  }
227  if(current == w)
228  break; // stop when this was the last segment
229 
230  // initialize the new segment
231  begin = current;
232  current_label = *ilabels;
233  point = *is;
234  apex_height = partialSquaredMagnitude(point, dimension, pixel_pitch);
235  Stack(1, Influence(DestType(0), 0.0, begin-1.0, begin-1.0, w)).swap(_stack);
236  // don't advance to next pixel here, because the present pixel must also
237  // be analysed in the context of the new segment
238  }
239  }
240 }
241 
242 template <unsigned int N, class T1, class S1,
243  class T2, class S2,
244  class Array>
245 void
246 interpixelBoundaryVectorDistance(MultiArrayView<N, T1, S1> const & labels,
247  MultiArrayView<N, T2, S2> dest,
248  Array const & pixelPitch)
249 {
250  typedef typename MultiArrayShape<N>::type Shape;
251  typedef GridGraph<N> Graph;
252  typedef typename Graph::Node Node;
253  typedef typename Graph::NodeIt graph_scanner;
254  typedef typename Graph::OutArcIt neighbor_iterator;
255 
256  Graph g(labels.shape());
257  for (graph_scanner node(g); node != lemon_graph::INVALID; ++node)
258  {
259  T1 label = labels[*node];
260  double min_dist = NumericTraits<double>::max();
261  Node point = *node,
262  boundary = point + Node(dest[point]),
263  min_pos = lemon::INVALID;
264  T2 min_diff;
265 
266  //go to adjacent neighbour with same label as origin pixel with smallest distance
267  if(labels.isInside(boundary))
268  {
269  for (neighbor_iterator arc(g, boundary); arc != lemon_graph::INVALID; ++arc)
270  {
271  if(label == labels[g.target(*arc)])
272  {
273  double dist = squaredNorm(pixelPitch*(g.target(*arc) - point));
274  if (dist < min_dist)
275  {
276  min_dist = dist;
277  min_pos = g.target(*arc);
278  }
279  }
280  }
281  if(min_pos == lemon::INVALID)
282  continue;
283  min_dist = NumericTraits<double>::max();
284  }
285  else
286  {
287  min_pos = clip(boundary, Shape(0), labels.shape()-Shape(1));
288  min_diff = 0.5*(boundary + min_pos) - point;
289  min_dist = squaredNorm(pixelPitch*min_diff);
290  }
291 
292  //from this pixel look for the vector which points to the nearest interpixel between two label
293  for (neighbor_iterator arc(g, min_pos); arc != lemon_graph::INVALID; ++arc)
294  {
295  if(label != labels[g.target(*arc)])
296  {
297  T2 diff = 0.5*(g.target(*arc) + min_pos) - point;
298  double dist = squaredNorm(pixelPitch*diff);
299  if (dist < min_dist)
300  {
301  min_dist = dist;
302  min_diff = diff;
303  }
304  }
305  }
306  dest[point] = min_diff;
307  }
308 }
309 
310 } // namespace detail
311 
312 /** \addtogroup MultiArrayDistanceTransform
313 */
314 //@{
315 
316 template<bool PRED>
317 struct Error_output_pixel_type_must_be_TinyVector_of_appropriate_length
318 : vigra::staticAssert::AssertBool<PRED> {};
319 
320 /********************************************************/
321 /* */
322 /* separableVectorDistance */
323 /* */
324 /********************************************************/
325 
326  /** \brief Compute the vector distance transform of a N-dimensional binary array.
327 
328  <b> Declarations:</b>
329 
330  \code
331  namespace vigra {
332  template <unsigned int N, class T1, class S1,
333  class T2, class S2, class Array>
334  void
335  separableVectorDistance(MultiArrayView<N, T1, S1> const & source,
336  MultiArrayView<N, T2, S2> dest,
337  bool background,
338  Array const & pixelPitch=TinyVector<double, N>(1));
339  }
340  \endcode
341 
342  This function works like \ref separableMultiDistance() (see there for details),
343  but returns in each pixel the <i>vector</i> to the nearest background pixel
344  rather than the scalar distance. This enables much more powerful applications.
345 
346  <b> Usage:</b>
347 
348  <b>\#include</b> <vigra/vector_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, Shape3> dest(shape);
355  ...
356 
357  // For each background pixel, find the vector to the nearest foreground pixel.
358  separableVectorDistance(source, dest, true);
359  \endcode
360 
361  \see vigra::separableMultiDistance(), vigra::boundaryVectorDistance()
362  */
364 
365 template <unsigned int N, class T1, class S1,
366  class T2, class S2, class Array>
367 void
368 separableVectorDistance(MultiArrayView<N, T1, S1> const & source,
369  MultiArrayView<N, T2, S2> dest,
370  bool background,
371  Array const & pixelPitch)
372 {
373  using namespace vigra::functor;
374  typedef typename MultiArrayView<N, T2, S2>::traverser Traverser;
375  typedef MultiArrayNavigator<Traverser, N> Navigator;
376 
377  VIGRA_STATIC_ASSERT((Error_output_pixel_type_must_be_TinyVector_of_appropriate_length<N == T2::static_size>));
378  vigra_precondition(source.shape() == dest.shape(),
379  "separableVectorDistance(): shape mismatch between input and output.");
380  vigra_precondition(pixelPitch.size() == N,
381  "separableVectorDistance(): pixelPitch has wrong length.");
382 
383  T2 maxDist(2*sum(source.shape()*pixelPitch)), rzero;
384  if(background == true)
385  transformMultiArray( source, dest,
386  ifThenElse( Arg1() == Param(0), Param(maxDist), Param(rzero) ));
387  else
388  transformMultiArray( source, dest,
389  ifThenElse( Arg1() != Param(0), Param(maxDist), Param(rzero) ));
390 
391  for(int d = 0; d < N; ++d )
392  {
393  Navigator nav( dest.traverser_begin(), dest.shape(), d);
394  for( ; nav.hasMore(); nav++ )
395  {
396  detail::vectorialDistParabola(d, nav.begin(), nav.end(), pixelPitch);
397  }
398  }
399 }
400 
401 template <unsigned int N, class T1, class S1,
402  class T2, class S2>
403 inline void
404 separableVectorDistance(MultiArrayView<N, T1, S1> const & source,
405  MultiArrayView<N, T2, S2> dest,
406  bool background=true)
407 {
408  TinyVector<double, N> pixelPitch(1.0);
409  separableVectorDistance(source, dest, background, pixelPitch);
410 }
411 
412 
413  /** \brief Compute the vector distance transform to the implicit boundaries of a
414  multi-dimensional label array.
415 
416  <b> Declarations:</b>
417 
418  \code
419  namespace vigra {
420  template <unsigned int N, class T1, class S1,
421  class T2, class S2,
422  class Array>
423  void
424  boundaryVectorDistance(MultiArrayView<N, T1, S1> const & labels,
425  MultiArrayView<N, T2, S2> dest,
426  bool array_border_is_active=false,
427  BoundaryDistanceTag boundary=OuterBoundary,
428  Array const & pixelPitch=TinyVector<double, N>(1));
429  }
430  \endcode
431 
432  This function works like \ref boundaryMultiDistance() (see there for details),
433  but returns in each pixel the <i>vector</i> to the nearest boundary pixel
434  rather than the scalar distance. This enables much more powerful applications.
435  Additionally, it support a <tt>pixelPitch</tt> parameter which allows to adjust
436  the distance calculations for anisotropic grid resolution.
437 
438  <b> Usage:</b>
439 
440  <b>\#include</b> <vigra/vector_distance.hxx><br/>
441  Namespace: vigra
442 
443  \code
444  Shape3 shape(width, height, depth);
445  MultiArray<3, UInt32> labels(shape);
446  MultiArray<3, Shape3> dest(shape);
447  ...
448 
449  // For each region, find the vectors to the nearest boundary pixel, including the
450  // outer border of the array.
451  boundaryVectorDistance(labels, dest, true);
452  \endcode
453 
454  \see vigra::boundaryMultiDistance(), vigra::separableVectorDistance()
455  */
457 
458 template <unsigned int N, class T1, class S1,
459  class T2, class S2,
460  class Array>
461 void
462 boundaryVectorDistance(MultiArrayView<N, T1, S1> const & labels,
463  MultiArrayView<N, T2, S2> dest,
464  bool array_border_is_active,
465  BoundaryDistanceTag boundary,
466  Array const & pixelPitch)
467 {
468  VIGRA_STATIC_ASSERT((Error_output_pixel_type_must_be_TinyVector_of_appropriate_length<N == T2::static_size>));
469  vigra_precondition(labels.shape() == dest.shape(),
470  "boundaryVectorDistance(): shape mismatch between input and output.");
471  vigra_precondition(pixelPitch.size() == N,
472  "boundaryVectorDistance(): pixelPitch has wrong length.");
473 
474  using namespace vigra::functor;
475 
476  if(boundary == InnerBoundary)
477  {
478  MultiArray<N, unsigned char> boundaries(labels.shape());
479 
480  markRegionBoundaries(labels, boundaries, IndirectNeighborhood);
481  if(array_border_is_active)
482  initMultiArrayBorder(boundaries, 1, 1);
483  separableVectorDistance(boundaries, dest, true, pixelPitch);
484  }
485  else
486  {
487  if(boundary == InterpixelBoundary)
488  {
489  vigra_precondition(!NumericTraits<T2>::isIntegral::value,
490  "boundaryVectorDistance(..., InterpixelBoundary): output pixel type must be float or double.");
491  }
492 
493  typedef typename MultiArrayView<N, T1, S1>::const_traverser LabelIterator;
494  typedef typename MultiArrayView<N, T2, S2>::traverser DestIterator;
495  typedef MultiArrayNavigator<LabelIterator, N> LabelNavigator;
496  typedef MultiArrayNavigator<DestIterator, N> DNavigator;
497 
498  T2 maxDist(2*sum(labels.shape()*pixelPitch));
499  dest = maxDist;
500  for( int d = 0; d < N; ++d )
501  {
502  LabelNavigator lnav( labels.traverser_begin(), labels.shape(), d );
503  DNavigator dnav( dest.traverser_begin(), dest.shape(), d );
504 
505  for( ; dnav.hasMore(); dnav++, lnav++ )
506  {
507  detail::boundaryVectorDistParabola(d, dnav.begin(), dnav.end(), lnav.begin(),
508  pixelPitch, maxDist, array_border_is_active);
509  }
510  }
511 
512  if(boundary == InterpixelBoundary)
513  {
514  detail::interpixelBoundaryVectorDistance(labels, dest, pixelPitch);
515  }
516  }
517 }
518 
519 template <unsigned int N, class T1, class S1,
520  class T2, class S2>
521 void
522 boundaryVectorDistance(MultiArrayView<N, T1, S1> const & labels,
523  MultiArrayView<N, T2, S2> dest,
524  bool array_border_is_active=false,
526 {
527  TinyVector<double, N> pixelPitch(1.0);
528  boundaryVectorDistance(labels, dest, array_border_is_active, boundary, pixelPitch);
529 }
530 
531 //@}
532 
533 } //-- namespace vigra
534 
535 #endif //-- VIGRA_VECTOR_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.10.0 (Thu Jan 8 2015)