37 #ifndef VIGRA_MULTI_DISTANCE_HXX 38 #define VIGRA_MULTI_DISTANCE_HXX 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" 51 #include "multi_gridgraph.hxx" 52 #include "union_find.hxx" 60 template <
class Value>
61 struct DistParabolaStackEntry
63 double left, center, right;
66 DistParabolaStackEntry(Value
const & p,
double l,
double c,
double r)
67 : left(l), center(c), right(r), apex_height(p)
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 )
89 double sigma2 = sigma * sigma;
90 double sigma22 = 2.0 * sigma2;
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));
99 for(;current < w; ++is, ++current)
105 Influence & s = _stack.back();
106 double diff = current - s.center;
107 intersection = current + (sa(is) - s.apex_height - sigma2*
sq(diff)) / (sigma22 * diff);
109 if( intersection < s.left)
117 else if(intersection < s.right)
119 s.right = intersection;
123 _stack.push_back(Influence(sa(is), intersection, current, w));
129 typename std::vector<Influence>::iterator it = _stack.begin();
130 for(current = 0.0; current < w; ++current, ++id)
132 while( current >= it->right)
134 da.set(sigma2 *
sq(current - it->center) + it->apex_height,
id);
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)
143 distParabola(src.first, src.second, src.third,
144 dest.first, dest.second, sigma);
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)
162 enum { N = SrcShape::static_size};
165 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
168 ArrayVector<TmpType> tmp( shape[0] );
170 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
171 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
175 SNavigator snav( si, shape, 0 );
176 DNavigator dnav( di, shape, 0 );
180 for( ; snav.hasMore(); snav++, dnav++ )
185 transformLine( snav.begin(), snav.end(), src, tmp.begin(),
186 typename AccessorTraits<TmpType>::default_accessor(),
187 Param(NumericTraits<TmpType>::zero())-Arg1());
189 copyLine( snav.begin(), snav.end(), src, tmp.begin(),
190 typename AccessorTraits<TmpType>::default_accessor() );
192 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
193 typename AccessorTraits<TmpType>::default_const_accessor()),
194 destIter( dnav.begin(), dest ), sigmas[0] );
198 for(
int d = 1; d < N; ++d )
200 DNavigator dnav( di, shape, d );
202 tmp.resize( shape[d] );
205 for( ; dnav.hasMore(); dnav++ )
208 copyLine( dnav.begin(), dnav.end(), dest,
209 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
211 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
212 typename AccessorTraits<TmpType>::default_const_accessor()),
213 destIter( dnav.begin(), dest ), sigmas[d] );
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)
224 internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas,
false );
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)
232 ArrayVector<double> sigmas(shape.size(), 1.0);
233 internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas,
false );
368 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
369 class DestIterator,
class DestAccessor,
class Array>
371 DestIterator d, DestAccessor dest,
bool background,
372 Array
const & pixelPitch)
374 int N = shape.size();
376 typedef typename SrcAccessor::value_type SrcType;
377 typedef typename DestAccessor::value_type DestType;
378 typedef typename NumericTraits<DestType>::RealPromote Real;
380 SrcType zero = NumericTraits<SrcType>::zero();
383 bool pixelPitchIsReal =
false;
384 for(
int k=0; k<N; ++k)
386 if(
int(pixelPitch[k]) != pixelPitch[k])
387 pixelPitchIsReal =
true;
388 dmax +=
sq(pixelPitch[k]*shape[k]);
393 if(dmax > NumericTraits<DestType>::toRealPromote(NumericTraits<DestType>::max())
397 Real maxDist = (Real)dmax, rzero = (Real)0.0;
398 MultiArray<SrcShape::static_size, Real> tmpArray(shape);
399 if(background ==
true)
401 tmpArray.traverser_begin(),
typename AccessorTraits<Real>::default_accessor(),
402 ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
405 tmpArray.traverser_begin(),
typename AccessorTraits<Real>::default_accessor(),
406 ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
408 detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(),
409 shape,
typename AccessorTraits<Real>::default_accessor(),
410 tmpArray.traverser_begin(),
411 typename AccessorTraits<Real>::default_accessor(), pixelPitch);
418 DestType maxDist = DestType(
std::ceil(dmax)), rzero = (DestType)0;
419 if(background ==
true)
421 ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
424 ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
426 detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest, pixelPitch);
430 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
431 class DestIterator,
class DestAccessor>
434 DestIterator d, DestAccessor dest,
bool background)
436 ArrayVector<double> pixelPitch(shape.size(), 1.0);
440 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
441 class DestIterator,
class DestAccessor,
class Array>
443 pair<DestIterator, DestAccessor>
const & dest,
bool background,
444 Array
const & pixelPitch)
447 dest.first, dest.second, background, pixelPitch );
450 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
451 class DestIterator,
class DestAccessor>
453 pair<DestIterator, DestAccessor>
const & dest,
bool background)
456 dest.first, dest.second, background );
459 template <
unsigned int N,
class T1,
class S1,
464 MultiArrayView<N, T2, S2> dest,
bool background,
465 Array
const & pixelPitch)
467 vigra_precondition(source.shape() == dest.shape(),
468 "separableMultiDistSquared(): shape mismatch between input and output.");
470 destMultiArray(dest), background, pixelPitch );
473 template <
unsigned int N,
class T1,
class S1,
477 MultiArrayView<N, T2, S2> dest,
bool background)
479 vigra_precondition(source.shape() == dest.shape(),
480 "separableMultiDistSquared(): shape mismatch between input and output.");
482 destMultiArray(dest), background );
588 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
589 class DestIterator,
class DestAccessor,
class Array>
591 DestIterator d, DestAccessor dest,
bool background,
592 Array
const & pixelPitch)
602 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
603 class DestIterator,
class DestAccessor>
605 DestIterator d, DestAccessor dest,
bool background)
615 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
616 class DestIterator,
class DestAccessor,
class Array>
618 pair<DestIterator, DestAccessor>
const & dest,
bool background,
619 Array
const & pixelPitch)
622 dest.first, dest.second, background, pixelPitch );
625 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
626 class DestIterator,
class DestAccessor>
628 pair<DestIterator, DestAccessor>
const & dest,
bool background)
631 dest.first, dest.second, background );
634 template <
unsigned int N,
class T1,
class S1,
635 class T2,
class S2,
class Array>
638 MultiArrayView<N, T2, S2> dest,
640 Array
const & pixelPitch)
642 vigra_precondition(source.shape() == dest.shape(),
643 "separableMultiDistance(): shape mismatch between input and output.");
645 destMultiArray(dest), background, pixelPitch );
648 template <
unsigned int N,
class T1,
class S1,
652 MultiArrayView<N, T2, S2> dest,
655 vigra_precondition(source.shape() == dest.shape(),
656 "separableMultiDistance(): shape mismatch between input and output.");
658 destMultiArray(dest), background );
664 namespace lemon_graph {
666 template <
class Graph,
class T1Map,
class T2Map>
668 markRegionBoundaries(Graph
const & g,
669 T1Map
const & labels,
672 typedef typename Graph::NodeIt graph_scanner;
673 typedef typename Graph::OutBackArcIt neighbor_iterator;
676 for (graph_scanner node(g); node != INVALID; ++node)
678 typename T1Map::value_type center = labels[*node];
680 for (neighbor_iterator arc(g, node); arc != INVALID; ++arc)
683 if(center != labels[g.target(*arc)])
686 out[g.target(*arc)] = 1;
696 template <
unsigned int N,
class T1,
class S1,
699 markRegionBoundaries(MultiArrayView<N, T1, S1>
const & labels,
700 MultiArrayView<N, T2, S2> out,
703 vigra_precondition(labels.shape() == out.shape(),
704 "markRegionBoundaries(): shape mismatch between input and output.");
706 GridGraph<N, undirected_tag> graph(labels.shape(), neighborhood);
708 lemon_graph::markRegionBoundaries(graph, labels, out);
722 template <
class DestIterator,
class LabelIterator>
724 boundaryDistParabola(DestIterator is, DestIterator iend,
725 LabelIterator ilabels,
727 bool array_border_is_active=
false)
730 double w = iend - is;
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;
740 double apex_height = array_border_is_active
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)
747 apex_height = (current < w)
748 ? (current_label == *ilabels)
751 : array_border_is_active
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);
760 if(intersection < s.left)
764 intersection = begin;
768 else if(intersection < s.right)
770 s.right = intersection;
773 _stack.push_back(Influence(apex_height, intersection, current, w));
774 if(current < w && current_label == *ilabels)
778 typename Stack::iterator it = _stack.begin();
779 for(
double c = begin; c < current; ++c, ++id)
781 while(c >= it->right)
783 *
id =
sq(c - it->center) + it->apex_height;
790 current_label = *ilabels;
792 Stack(1, Influence(0.0, begin-1.0, begin-1.0, w)).swap(_stack);
805 template <
unsigned int N,
class T1,
class S1,
808 internalBoundaryMultiArrayDist(
809 MultiArrayView<N, T1, S1>
const & labels,
810 MultiArrayView<N, T2, S2> dest,
811 double dmax,
bool array_border_is_active=
false)
815 typedef MultiArrayNavigator<LabelIterator, N> LabelNavigator;
816 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
819 for(
int d = 0; d < N; ++d )
821 LabelNavigator lnav( labels.traverser_begin(), labels.shape(), d );
822 DNavigator dnav( dest.traverser_begin(), dest.shape(), d );
824 for( ; dnav.hasMore(); dnav++, lnav++ )
826 boundaryDistParabola(dnav.begin(), dnav.end(),
828 dmax, array_border_is_active);
909 template <
unsigned int N,
class T1,
class S1,
914 bool array_border_is_active=
false,
917 vigra_precondition(labels.
shape() == dest.
shape(),
918 "boundaryMultiDistance(): shape mismatch between input and output.");
927 if(array_border_is_active)
937 vigra_precondition(!NumericTraits<T2>::isIntegral::value,
938 "boundaryMultiDistance(..., InterpixelBoundary): output pixel type must be float or double.");
942 if(dmax >
double(NumericTraits<T2>::max()))
945 typedef typename NumericTraits<T2>::RealPromote Real;
947 detail::internalBoundaryMultiArrayDist(labels, tmpArray,
948 dmax, array_border_is_active);
954 detail::internalBoundaryMultiArrayDist(labels, dest, dmax, array_border_is_active);
965 #endif //-- VIGRA_MULTI_DISTANCE_HXX
void initMultiArrayBorder(...)
Write values to the specified border values in the array.
const difference_type & shape() const
Definition: multi_array.hxx:1596
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
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
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.
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