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

vector_distance.hxx
1/************************************************************************/
2/* */
3/* Copyright 2003-2015 by Philipp Schubert and Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36#ifndef VIGRA_VECTOR_DISTANCE_HXX
37#define VIGRA_VECTOR_DISTANCE_HXX
38
39#include <vector>
40#include <set>
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#include "multi_distance.hxx"
51
52#undef VECTORIAL_DIST_DEBUG
53
54namespace vigra
55{
56
57namespace detail
58{
59
60template <class Vector, class Value>
61struct VectorialDistParabolaStackEntry
62{
63 double left, center, right;
64 Value apex_height;
65 Vector point;
66
67 VectorialDistParabolaStackEntry(const Vector& vec, Value prev, double l, double c, double r)
68 : left(l), center(c), right(r), apex_height(prev), point(vec)
69 {}
70};
71
72#ifdef VECTORIAL_DIST_DEBUG
73template <class Vector, class Value>
74std::ostream& operator<<(std::ostream&o, const VectorialDistParabolaStackEntry<Vector, Value>& e) {
75 o << "l=" << e.left << ", c=" << e.center << ", r=" << e.right << ", pV=" << e.apex_height << ", pVec=" << e.point;
76 return o;
77}
78#endif
79
80/********************************************************/
81/* */
82/* vectorialDistParabola */
83/* */
84/********************************************************/
85
86template <class VEC, class ARRAY>
87inline double
88partialSquaredMagnitude(const VEC& vec, MultiArrayIndex dim, ARRAY const & pixel_pitch)
89{
90 //computes the squared magnitude of vec
91 //considering only the first dim dimensions
92 double sqMag = 0.0;
93 for(MultiArrayIndex i=0; i<=dim; ++i)
94 {
95 sqMag += sq(pixel_pitch[i]*vec[i]);
96 }
97 return sqMag;
98}
99
100template <class SrcIterator,
101 class Array>
102void
103vectorialDistParabola(MultiArrayIndex dimension,
104 SrcIterator is, SrcIterator iend,
105 Array const & pixel_pitch )
106{
107 typedef typename SrcIterator::value_type SrcType;
108 typedef VectorialDistParabolaStackEntry<SrcType, double> Influence;
109
110 double sigma = pixel_pitch[dimension],
111 sigma2 = sq(sigma);
112 double w = iend - is; //width of the scanline
113
114 SrcIterator id = is;
115
116 std::vector<Influence> _stack; //stack of influence parabolas
117 double apex_height = partialSquaredMagnitude(*is, dimension, pixel_pitch);
118 _stack.push_back(Influence(*is, apex_height, 0.0, 0.0, w));
119 ++is;
120 double current = 1.0;
121 while(current < w)
122 {
123 apex_height = partialSquaredMagnitude(*is, dimension, pixel_pitch);
124 Influence & s = _stack.back();
125 double diff = current - s.center;
126 double intersection = current + (apex_height - s.apex_height - sq(sigma*diff)) / (2.0*sigma2 * diff);
127
128 if( intersection < s.left) // previous point has no influence
129 {
130 _stack.pop_back();
131 if(_stack.empty())
132 _stack.push_back(Influence(*is, apex_height, 0.0, current, w));
133 else
134 continue; // try new top of stack without advancing current
135 }
136 else if(intersection < s.right)
137 {
138 s.right = intersection;
139 _stack.push_back(Influence(*is, apex_height, intersection, current, w));
140 }
141 ++is;
142 ++current;
143 }
144
145 // Now we have the stack indicating which rows are influenced by (and therefore
146 // closest to) which row. We can go through the stack and calculate the
147 // distance squared for each element of the column.
148 typename std::vector<Influence>::iterator it = _stack.begin();
149 for(current = 0.0; current < w; ++current, ++id)
150 {
151 while( current >= it->right)
152 ++it;
153
154 *id = it->point;
155 (*id)[dimension] = it->center - current;
156 }
157}
158
159template <class DestIterator,
160 class LabelIterator,
161 class Array1, class Array2>
162void
163boundaryVectorDistParabola(MultiArrayIndex dimension,
164 DestIterator is, DestIterator iend,
165 LabelIterator ilabels,
166 Array1 const & pixel_pitch,
167 Array2 const & dmax,
168 bool array_border_is_active=false)
169{
170 double w = iend - is;
171 if(w <= 0)
172 return;
173
174 typedef typename LabelIterator::value_type LabelType;
175 typedef typename DestIterator::value_type DestType;
176 typedef VectorialDistParabolaStackEntry<DestType, double> Influence;
177 typedef std::vector<Influence> Stack;
178
179 DestIterator id = is;
180 DestType border_point = array_border_is_active
181 ? DestType(0)
182 : dmax;
183 double apex_height = partialSquaredMagnitude(border_point, dimension, pixel_pitch);
184 Stack _stack(1, Influence(border_point, apex_height, 0.0, -1.0, w));
185 LabelType current_label = *ilabels;
186 for(double begin = 0.0, current = 0.0; current <= w; ++ilabels, ++is, ++current)
187 {
188 DestType point = (current < w)
189 ? (current_label == *ilabels)
190 ? *is
191 : DestType(0)
192 : border_point;
193 apex_height = partialSquaredMagnitude(point, dimension, pixel_pitch);
194 while(true)
195 {
196 Influence & s = _stack.back();
197 double diff = (current - s.center)*pixel_pitch[dimension];
198 double intersection = current + (apex_height - s.apex_height - sq(diff)) / (2.0 * diff);
199
200 if(intersection < s.left) // previous parabola has no influence
201 {
202 _stack.pop_back();
203 if(_stack.empty())
204 intersection = begin; // new parabola is valid for entire present segment
205 else
206 continue; // try new top of stack without advancing to next pixel
207 }
208 else if(intersection < s.right)
209 {
210 s.right = intersection;
211 }
212 if(intersection < w)
213 _stack.push_back(Influence(point, apex_height, intersection, current, w));
214 if(current < w && current_label == *ilabels)
215 break; // finished present pixel, advance to next one
216
217 // label changed => finalize the current segment
218 typename Stack::iterator it = _stack.begin();
219 for(double c = begin; c < current; ++c, ++id)
220 {
221 while(c >= it->right)
222 ++it;
223 *id = it->point;
224 (*id)[dimension] = it->center - c;
225 }
226 if(current == w)
227 break; // stop when this was the last segment
228
229 // initialize the new segment
230 begin = current;
231 current_label = *ilabels;
232 point = *is;
233 apex_height = partialSquaredMagnitude(point, dimension, pixel_pitch);
234 Stack(1, Influence(DestType(0), 0.0, begin-1.0, begin-1.0, w)).swap(_stack);
235 // don't advance to next pixel here, because the present pixel must also
236 // be analysed in the context of the new segment
237 }
238 }
239}
240
241template <unsigned int N, class T1, class S1,
242 class T2, class S2,
243 class Array>
244void
245interpixelBoundaryVectorDistance(MultiArrayView<N, T1, S1> const & labels,
246 MultiArrayView<N, T2, S2> dest,
247 Array const & pixelPitch)
248{
249 typedef typename MultiArrayShape<N>::type Shape;
250 typedef GridGraph<N> Graph;
251 typedef typename Graph::Node Node;
252 typedef typename Graph::NodeIt graph_scanner;
253 typedef typename Graph::OutArcIt neighbor_iterator;
254
255 Graph g(labels.shape());
256 for (graph_scanner node(g); node != lemon_graph::INVALID; ++node)
257 {
258 T1 label = labels[*node];
259 double min_dist = NumericTraits<double>::max();
260 Node point = *node,
261 boundary = point + Node(dest[point]),
262 min_pos = lemon::INVALID;
263 T2 min_diff;
264
265 //go to adjacent neighbour with same label as origin pixel with smallest distance
266 if(labels.isInside(boundary))
267 {
268 for (neighbor_iterator arc(g, boundary); arc != lemon_graph::INVALID; ++arc)
269 {
270 if(label == labels[g.target(*arc)])
271 {
272 double dist = squaredNorm(pixelPitch*(g.target(*arc) - point));
273 if (dist < min_dist)
274 {
275 min_dist = dist;
276 min_pos = g.target(*arc);
277 }
278 }
279 }
280 if(min_pos == lemon::INVALID)
281 continue;
282 min_dist = NumericTraits<double>::max();
283 }
284 else
285 {
286 min_pos = clip(boundary, Shape(0), labels.shape()-Shape(1));
287 min_diff = 0.5*(boundary + min_pos) - point;
288 min_dist = squaredNorm(pixelPitch*min_diff);
289 }
290
291 //from this pixel look for the vector which points to the nearest interpixel between two label
292 for (neighbor_iterator arc(g, min_pos); arc != lemon_graph::INVALID; ++arc)
293 {
294 if(label != labels[g.target(*arc)])
295 {
296 T2 diff = 0.5*(g.target(*arc) + min_pos) - point;
297 double dist = squaredNorm(pixelPitch*diff);
298 if (dist < min_dist)
299 {
300 min_dist = dist;
301 min_diff = diff;
302 }
303 }
304 }
305 dest[point] = min_diff;
306 }
307}
308
309} // namespace detail
310
311/** \addtogroup DistanceTransform
312*/
313//@{
314
315template<bool PRED>
316struct Error_output_pixel_type_must_be_TinyVector_of_appropriate_length
317: vigra::staticAssert::AssertBool<PRED> {};
318
319/********************************************************/
320/* */
321/* separableVectorDistance */
322/* */
323/********************************************************/
324
325 /** \brief Compute the vector distance transform of a N-dimensional binary array.
326
327 <b> Declarations:</b>
328
329 \code
330 namespace vigra {
331 template <unsigned int N, class T1, class S1,
332 class T2, class S2, class Array>
333 void
334 separableVectorDistance(MultiArrayView<N, T1, S1> const & source,
335 MultiArrayView<N, T2, S2> dest,
336 bool background,
337 Array const & pixelPitch=TinyVector<double, N>(1));
338 }
339 \endcode
340
341 This function works like \ref separableMultiDistance() (see there for details),
342 but returns in each pixel the <i>vector</i> to the nearest background pixel
343 rather than the scalar distance. This enables much more powerful applications.
344
345 <b> Usage:</b>
346
347 <b>\#include</b> <vigra/vector_distance.hxx><br/>
348 Namespace: vigra
349
350 \code
351 Shape3 shape(width, height, depth);
352 MultiArray<3, unsigned char> source(shape);
353 MultiArray<3, Shape3> dest(shape);
354 ...
355
356 // For each background pixel, find the vector to the nearest foreground pixel.
357 separableVectorDistance(source, dest, true);
358 \endcode
359
360 \see vigra::separableMultiDistance(), vigra::boundaryVectorDistance()
361 */
363
364template <unsigned int N, class T1, class S1,
365 class T2, class S2, class Array>
366void
369 bool background,
370 Array const & pixelPitch)
371{
372 using namespace vigra::functor;
373 typedef typename MultiArrayView<N, T2, S2>::traverser Traverser;
374 typedef MultiArrayNavigator<Traverser, N> Navigator;
375
376 VIGRA_STATIC_ASSERT((Error_output_pixel_type_must_be_TinyVector_of_appropriate_length<N == T2::static_size>));
377 vigra_precondition(source.shape() == dest.shape(),
378 "separableVectorDistance(): shape mismatch between input and output.");
379 vigra_precondition(pixelPitch.size() == N,
380 "separableVectorDistance(): pixelPitch has wrong length.");
381
382 T2 maxDist(2*sum(source.shape()*pixelPitch)), rzero;
383 if(background == true)
384 transformMultiArray( source, dest,
385 ifThenElse( Arg1() == Param(0), Param(maxDist), Param(rzero) ));
386 else
387 transformMultiArray( source, dest,
388 ifThenElse( Arg1() != Param(0), Param(maxDist), Param(rzero) ));
389
390 for(unsigned d = 0; d < N; ++d )
391 {
392 Navigator nav( dest.traverser_begin(), dest.shape(), d);
393 for( ; nav.hasMore(); nav++ )
394 {
395 detail::vectorialDistParabola(d, nav.begin(), nav.end(), pixelPitch);
396 }
397 }
398}
399
400template <unsigned int N, class T1, class S1,
401 class T2, class S2>
402inline void
405 bool background=true)
406{
407 TinyVector<double, N> pixelPitch(1.0);
408 separableVectorDistance(source, dest, background, pixelPitch);
409}
410
411
412 /** \brief Compute the vector distance transform to the implicit boundaries of a
413 multi-dimensional label array.
414
415 <b> Declarations:</b>
416
417 \code
418 namespace vigra {
419 template <unsigned int N, class T1, class S1,
420 class T2, class S2,
421 class Array>
422 void
423 boundaryVectorDistance(MultiArrayView<N, T1, S1> const & labels,
424 MultiArrayView<N, T2, S2> dest,
425 bool array_border_is_active=false,
426 BoundaryDistanceTag boundary=OuterBoundary,
427 Array const & pixelPitch=TinyVector<double, N>(1));
428 }
429 \endcode
430
431 This function works like \ref boundaryMultiDistance() (see there for details),
432 but returns in each pixel the <i>vector</i> to the nearest boundary pixel
433 rather than the scalar distance. This enables much more powerful applications.
434 Additionally, it support a <tt>pixelPitch</tt> parameter which allows to adjust
435 the distance calculations for anisotropic grid resolution.
436
437 <b> Usage:</b>
438
439 <b>\#include</b> <vigra/vector_distance.hxx><br/>
440 Namespace: vigra
441
442 \code
443 Shape3 shape(width, height, depth);
444 MultiArray<3, UInt32> labels(shape);
445 MultiArray<3, Shape3> dest(shape);
446 ...
447
448 // For each region, find the vectors to the nearest boundary pixel, including the
449 // outer border of the array.
450 boundaryVectorDistance(labels, dest, true);
451 \endcode
452
453 \see vigra::boundaryMultiDistance(), vigra::separableVectorDistance()
454 */
456
457template <unsigned int N, class T1, class S1,
458 class T2, class S2,
459 class Array>
460void
463 bool array_border_is_active,
464 BoundaryDistanceTag boundary,
465 Array const & pixelPitch)
466{
467 VIGRA_STATIC_ASSERT((Error_output_pixel_type_must_be_TinyVector_of_appropriate_length<N == T2::static_size>));
468 vigra_precondition(labels.shape() == dest.shape(),
469 "boundaryVectorDistance(): shape mismatch between input and output.");
470 vigra_precondition(pixelPitch.size() == N,
471 "boundaryVectorDistance(): pixelPitch has wrong length.");
472
473 using namespace vigra::functor;
474
475 if(boundary == InnerBoundary)
476 {
477 MultiArray<N, unsigned char> boundaries(labels.shape());
478
479 markRegionBoundaries(labels, boundaries, IndirectNeighborhood);
480 if(array_border_is_active)
481 initMultiArrayBorder(boundaries, 1, 1);
482 separableVectorDistance(boundaries, dest, true, pixelPitch);
483 }
484 else
485 {
486 if(boundary == InterpixelBoundary)
487 {
488 vigra_precondition(!NumericTraits<T2>::isIntegral::value,
489 "boundaryVectorDistance(..., InterpixelBoundary): output pixel type must be float or double.");
490 }
491
492 typedef typename MultiArrayView<N, T1, S1>::const_traverser LabelIterator;
493 typedef typename MultiArrayView<N, T2, S2>::traverser DestIterator;
494 typedef MultiArrayNavigator<LabelIterator, N> LabelNavigator;
495 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
496
497 T2 maxDist(2*sum(labels.shape()*pixelPitch));
498 dest = maxDist;
499 for( unsigned d = 0; d < N; ++d )
500 {
501 LabelNavigator lnav( labels.traverser_begin(), labels.shape(), d );
502 DNavigator dnav( dest.traverser_begin(), dest.shape(), d );
503
504 for( ; dnav.hasMore(); dnav++, lnav++ )
505 {
506 detail::boundaryVectorDistParabola(d, dnav.begin(), dnav.end(), lnav.begin(),
507 pixelPitch, maxDist, array_border_is_active);
508 }
509 }
510
511 if(boundary == InterpixelBoundary)
512 {
513 detail::interpixelBoundaryVectorDistance(labels, dest, pixelPitch);
514 }
515 }
516}
517
518template <unsigned int N, class T1, class S1,
519 class T2, class S2>
520void
523 bool array_border_is_active=false,
525{
526 TinyVector<double, N> pixelPitch(1.0);
527 boundaryVectorDistance(labels, dest, array_border_is_active, boundary, pixelPitch);
528}
529
530//@}
531
532} //-- namespace vigra
533
534#endif //-- VIGRA_VECTOR_DISTANCE_HXX
A navigator that provides access to the 1D subranges of an n-dimensional range given by a vigra::Mult...
Definition: navigator.hxx:101
TinyVector< MultiArrayIndex, N > type
Definition: multi_shape.hxx:272
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:705
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, Tconst &, Tconst * >::type const_traverser
Definition: multi_array.hxx:769
traverser traverser_begin()
Definition: multi_array.hxx:1953
const difference_type & shape() const
Definition: multi_array.hxx:1648
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, T &, T * >::type traverser
Definition: multi_array.hxx:764
Main MultiArray class containing the memory management.
Definition: multi_array.hxx:2477
BoundaryDistanceTag
Specify which boundary is used for boundaryMultiDistance().
Definition: multi_distance.hxx:835
@ OuterBoundary
Pixels just outside of each region.
Definition: multi_distance.hxx:836
@ InnerBoundary
Pixels just inside of each region.
Definition: multi_distance.hxx:838
@ InterpixelBoundary
Half-integer points between pixels of different labels.
Definition: multi_distance.hxx:837
void separableVectorDistance(...)
Compute the vector distance transform of a N-dimensional binary array.
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
doxygen_overloaded_function(template<... > void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
@ IndirectNeighborhood
use direct and indirect neighbors
Definition: multi_fwd.hxx:188
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
TinyVector< V, SIZE > clip(TinyVector< V, SIZE > const &t, const V valLower, const V valUpper)
Clip values to an interval.
Definition: tinyvector.hxx:2307
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
std::ptrdiff_t MultiArrayIndex
Definition: multi_fwd.hxx:60
void boundaryVectorDistance(...)
Compute the vector distance transform to the implicit boundaries of a multi-dimensional label array.
void initMultiArrayBorder(...)
Write values to the specified border values in the array.

© 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