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"
57 template <
class Value>
58 struct DistParabolaStackEntry
60 double left, center, right;
63 DistParabolaStackEntry(Value
const & p,
double l,
double c,
double r)
64 : left(l), center(c), right(r), prevVal(p)
76 template <
class SrcIterator,
class SrcAccessor,
77 class DestIterator,
class DestAccessor >
78 void distParabola(SrcIterator is, SrcIterator iend, SrcAccessor sa,
79 DestIterator
id, DestAccessor da,
double sigma )
86 double sigma2 = sigma * sigma;
87 double sigma22 = 2.0 * sigma2;
89 typedef typename SrcAccessor::value_type SrcType;
90 typedef DistParabolaStackEntry<SrcType> Influence;
91 std::vector<Influence> _stack;
92 _stack.push_back(Influence(sa(is), 0.0, 0.0, w));
98 Influence & s = _stack.back();
99 double diff = current - s.center;
100 double intersection = current + (sa(is) - s.prevVal - sigma2*
sq(diff)) / (sigma22 * diff);
102 if( intersection < s.left)
107 _stack.push_back(Influence(sa(is), 0.0, current, w));
114 else if(intersection < s.right)
116 s.right = intersection;
117 _stack.push_back(Influence(sa(is), intersection, current, w));
126 typename std::vector<Influence>::iterator it = _stack.begin();
127 for(current = 0.0; current < w; ++current, ++id)
129 while( current >= it->right)
131 da.set(sigma2 *
sq(current - it->center) + it->prevVal,
id);
135 template <
class SrcIterator,
class SrcAccessor,
136 class DestIterator,
class DestAccessor>
137 inline void distParabola(triple<SrcIterator, SrcIterator, SrcAccessor> src,
138 pair<DestIterator, DestAccessor> dest,
double sigma)
140 distParabola(src.first, src.second, src.third,
141 dest.first, dest.second, sigma);
150 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
151 class DestIterator,
class DestAccessor,
class Array>
152 void internalSeparableMultiArrayDistTmp(
153 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
154 DestIterator di, DestAccessor dest, Array
const & sigmas,
bool invert)
159 enum { N = SrcShape::static_size};
162 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
165 ArrayVector<TmpType> tmp( shape[0] );
167 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
168 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
172 SNavigator snav( si, shape, 0 );
173 DNavigator dnav( di, shape, 0 );
175 using namespace vigra::functor;
177 for( ; snav.hasMore(); snav++, dnav++ )
182 transformLine( snav.begin(), snav.end(), src, tmp.begin(),
183 typename AccessorTraits<TmpType>::default_accessor(),
184 Param(NumericTraits<TmpType>::zero())-Arg1());
186 copyLine( snav.begin(), snav.end(), src, tmp.begin(),
187 typename AccessorTraits<TmpType>::default_accessor() );
189 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
190 typename AccessorTraits<TmpType>::default_const_accessor()),
191 destIter( dnav.begin(), dest ), sigmas[0] );
195 for(
int d = 1; d < N; ++d )
197 DNavigator dnav( di, shape, d );
199 tmp.resize( shape[d] );
202 for( ; dnav.hasMore(); dnav++ )
205 copyLine( dnav.begin(), dnav.end(), dest,
206 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
208 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
209 typename AccessorTraits<TmpType>::default_const_accessor()),
210 destIter( dnav.begin(), dest ), sigmas[d] );
216 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
217 class DestIterator,
class DestAccessor,
class Array>
218 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape
const & shape, SrcAccessor src,
219 DestIterator di, DestAccessor dest, Array
const & sigmas)
221 internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas,
false );
224 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
225 class DestIterator,
class DestAccessor>
226 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape
const & shape, SrcAccessor src,
227 DestIterator di, DestAccessor dest)
229 ArrayVector<double> sigmas(shape.size(), 1.0);
230 internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas,
false );
367 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
368 class DestIterator,
class DestAccessor,
class Array>
370 DestIterator d, DestAccessor dest,
bool background,
371 Array
const & pixelPitch)
373 int N = shape.size();
375 typedef typename SrcAccessor::value_type SrcType;
376 typedef typename DestAccessor::value_type DestType;
377 typedef typename NumericTraits<DestType>::RealPromote Real;
379 SrcType zero = NumericTraits<SrcType>::zero();
382 bool pixelPitchIsReal =
false;
383 for(
int k=0; k<N; ++k)
385 if(
int(pixelPitch[k]) != pixelPitch[k])
386 pixelPitchIsReal =
true;
387 dmax +=
sq(pixelPitch[k]*shape[k]);
390 using namespace vigra::functor;
392 if(dmax > NumericTraits<DestType>::toRealPromote(NumericTraits<DestType>::max())
396 Real maxDist = (Real)dmax, rzero = (Real)0.0;
397 MultiArray<SrcShape::static_size, Real> tmpArray(shape);
398 if(background ==
true)
400 tmpArray.traverser_begin(),
typename AccessorTraits<Real>::default_accessor(),
401 ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
404 tmpArray.traverser_begin(),
typename AccessorTraits<Real>::default_accessor(),
405 ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
407 detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(),
408 shape,
typename AccessorTraits<Real>::default_accessor(),
409 tmpArray.traverser_begin(),
410 typename AccessorTraits<Real>::default_accessor(), pixelPitch);
417 DestType maxDist = DestType(
std::ceil(dmax)), rzero = (DestType)0;
418 if(background ==
true)
420 ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
423 ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
425 detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest, pixelPitch);
429 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
430 class DestIterator,
class DestAccessor>
433 DestIterator d, DestAccessor dest,
bool background)
435 ArrayVector<double> pixelPitch(shape.size(), 1.0);
439 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
440 class DestIterator,
class DestAccessor,
class Array>
442 pair<DestIterator, DestAccessor>
const & dest,
bool background,
443 Array
const & pixelPitch)
446 dest.first, dest.second, background, pixelPitch );
449 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
450 class DestIterator,
class DestAccessor>
452 pair<DestIterator, DestAccessor>
const & dest,
bool background)
455 dest.first, dest.second, background );
458 template <
unsigned int N,
class T1,
class S1,
463 MultiArrayView<N, T2, S2> dest,
bool background,
464 Array
const & pixelPitch)
466 vigra_precondition(source.shape() == dest.shape(),
467 "separableMultiDistSquared(): shape mismatch between input and output.");
469 destMultiArray(dest), background, pixelPitch );
472 template <
unsigned int N,
class T1,
class S1,
476 MultiArrayView<N, T2, S2> dest,
bool background)
478 vigra_precondition(source.shape() == dest.shape(),
479 "separableMultiDistSquared(): shape mismatch between input and output.");
481 destMultiArray(dest), background );
587 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
588 class DestIterator,
class DestAccessor,
class Array>
590 DestIterator d, DestAccessor dest,
bool background,
591 Array
const & pixelPitch)
596 using namespace vigra::functor;
601 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
602 class DestIterator,
class DestAccessor>
604 DestIterator d, DestAccessor dest,
bool background)
609 using namespace vigra::functor;
614 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
615 class DestIterator,
class DestAccessor,
class Array>
617 pair<DestIterator, DestAccessor>
const & dest,
bool background,
618 Array
const & pixelPitch)
621 dest.first, dest.second, background, pixelPitch );
624 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
625 class DestIterator,
class DestAccessor>
627 pair<DestIterator, DestAccessor>
const & dest,
bool background)
630 dest.first, dest.second, background );
633 template <
unsigned int N,
class T1,
class S1,
634 class T2,
class S2,
class Array>
637 MultiArrayView<N, T2, S2> dest,
639 Array
const & pixelPitch)
641 vigra_precondition(source.shape() == dest.shape(),
642 "separableMultiDistance(): shape mismatch between input and output.");
644 destMultiArray(dest), background, pixelPitch );
647 template <
unsigned int N,
class T1,
class S1,
651 MultiArrayView<N, T2, S2> dest,
654 vigra_precondition(source.shape() == dest.shape(),
655 "separableMultiDistance(): shape mismatch between input and output.");
657 destMultiArray(dest), background );
665 #endif //-- VIGRA_MULTI_DISTANCE_HXX