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

multi_distance.hxx
1/************************************************************************/
2/* */
3/* Copyright 2003-2007 by Kasim Terzic, Christian-Dennis Rahn, */
4/* 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
54namespace vigra
55{
56
57namespace detail
58{
59
60template <class Value>
61struct 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
79template <class SrcIterator, class SrcAccessor,
80 class DestIterator, class DestAccessor >
81void 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
138template <class SrcIterator, class SrcAccessor,
139 class DestIterator, class DestAccessor>
140inline 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
153template <class SrcIterator, class SrcShape, class SrcAccessor,
154 class DestIterator, class DestAccessor, class Array>
155void 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(),
187 Param(NumericTraits<TmpType>::zero())-Arg1());
188 else
189 copyLine( snav.begin(), snav.end(), src, tmp.begin(),
191
192 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
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,
210
211 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
213 destIter( dnav.begin(), dest ), sigmas[d] );
214 }
215 }
216 if(invert) transformMultiArray( di, shape, dest, di, dest, -Arg1());
217}
218
219template <class SrcIterator, class SrcShape, class SrcAccessor,
220 class DestIterator, class DestAccessor, class Array>
221inline 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
227template <class SrcIterator, class SrcShape, class SrcAccessor,
228 class DestIterator, class DestAccessor>
229inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape const & shape, SrcAccessor src,
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*/
363doxygen_overloaded_function(template <...> void separableMultiDistSquared)
364
365template <class SrcIterator, class SrcShape, class SrcAccessor,
366 class DestIterator, class DestAccessor, class Array>
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;
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(),
407 tmpArray.traverser_begin(),
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
427template <class SrcIterator, class SrcShape, class SrcAccessor,
428 class DestIterator, class DestAccessor>
429inline
432{
433 ArrayVector<double> pixelPitch(shape.size(), 1.0);
435}
436
437template <class SrcIterator, class SrcShape, class SrcAccessor,
438 class DestIterator, class DestAccessor, class Array>
441 Array const & pixelPitch)
442{
443 separableMultiDistSquared( source.first, source.second, source.third,
444 dest.first, dest.second, background, pixelPitch );
445}
446
447template <class SrcIterator, class SrcShape, class SrcAccessor,
448 class DestIterator, class DestAccessor>
451{
452 separableMultiDistSquared( source.first, source.second, source.third,
453 dest.first, dest.second, background );
454}
455
456template <unsigned int N, class T1, class S1,
457 class T2, class S2,
458 class Array>
459inline void
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
470template <unsigned int N, class T1, class S1,
471 class T2, class S2>
472inline void
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*/
583doxygen_overloaded_function(template <...> void separableMultiDistance)
584
585template <class SrcIterator, class SrcShape, class SrcAccessor,
586 class DestIterator, class DestAccessor, class Array>
589 Array const & pixelPitch)
590{
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
599template <class SrcIterator, class SrcShape, class SrcAccessor,
600 class DestIterator, class DestAccessor>
603{
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
612template <class SrcIterator, class SrcShape, class SrcAccessor,
613 class DestIterator, class DestAccessor, class Array>
616 Array const & pixelPitch)
617{
618 separableMultiDistance( source.first, source.second, source.third,
619 dest.first, dest.second, background, pixelPitch );
620}
621
622template <class SrcIterator, class SrcShape, class SrcAccessor,
623 class DestIterator, class DestAccessor>
626{
627 separableMultiDistance( source.first, source.second, source.third,
628 dest.first, dest.second, background );
629}
630
631template <unsigned int N, class T1, class S1,
632 class T2, class S2, class Array>
633inline void
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
645template <unsigned int N, class T1, class S1,
646 class T2, class S2>
647inline void
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
661namespace lemon_graph {
662
663template <class Graph, class T1Map, class T2Map>
664void
665markRegionBoundaries(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
691doxygen_overloaded_function(template <...> unsigned int markRegionBoundaries)
692
693template <unsigned int N, class T1, class S1,
694 class T2, class S2>
695inline void
696markRegionBoundaries(MultiArrayView<N, T1, S1> const & labels,
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
710namespace detail
711{
712
713/********************************************************/
714/* */
715/* boundaryDistParabola */
716/* */
717/********************************************************/
718
719template <class DestIterator, class LabelIterator>
720void
721boundaryDistParabola(DestIterator is, DestIterator iend,
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
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
802template <unsigned int N, class T1, class S1,
803 class T2, class S2>
804void
805internalBoundaryMultiArrayDist(
806 MultiArrayView<N, T1, S1> const & labels,
808 double dmax, bool array_border_is_active=false)
809{
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(),
826 }
827 }
828}
829
830} // namespace detail
831
832 /** \brief Specify which boundary is used for boundaryMultiDistance().
833
834 */
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
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*/
904doxygen_overloaded_function(template <...> void boundaryMultiDistance)
905
906template <unsigned int N, class T1, class S1,
907 class T2, class S2>
908void
911 bool array_border_is_active=false,
913{
914 vigra_precondition(labels.shape() == dest.shape(),
915 "boundaryMultiDistance(): shape mismatch between input and output.");
916
917 using namespace vigra::functor;
918
920 {
922
923 markRegionBoundaries(labels, boundaries, IndirectNeighborhood);
927 }
928 else
929 {
930 T2 offset = 0.0;
931
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,
947 }
948 else
949 {
950 // can work directly on the destination array
951 detail::internalBoundaryMultiArrayDist(labels, dest, dmax, array_border_is_active);
953 }
954 }
955}
956
957//@}
958
959} //-- namespace vigra
960
961
962#endif //-- VIGRA_MULTI_DISTANCE_HXX
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:1955
const difference_type & shape() const
Definition multi_array.hxx:1650
vigra::detail::MultiIteratorChooser< StrideTag >::template Traverser< actual_dimension, T, T &, T * >::type traverser
Definition multi_array.hxx:764
Class for a single RGB value.
Definition rgbvalue.hxx:128
size_type size() const
Definition tinyvector.hxx:913
iterator end()
Definition tinyvector.hxx:864
iterator begin()
Definition tinyvector.hxx:861
void separableMultiDistSquared(...)
Euclidean distance squared on multi-dimensional arrays.
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 separableMultiDistance(...)
Euclidean distance on multi-dimensional arrays.
void boundaryMultiDistance(...)
Euclidean distance to the implicit boundaries of a multi-dimensional label array.
NeighborhoodType
Choose the neighborhood system in a dimension-independent way.
Definition multi_fwd.hxx:186
@ IndirectNeighborhood
use direct and indirect neighbors
Definition multi_fwd.hxx:188
@ DirectNeighborhood
use only direct neighbors
Definition multi_fwd.hxx:187
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
void copyMultiArray(...)
Copy a multi-dimensional 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.12.2