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 )
83 double sigma2 = sigma * sigma;
84 double sigma22 = 2.0 * sigma2;
86 typedef typename SrcAccessor::value_type SrcType;
87 typedef DistParabolaStackEntry<SrcType> Influence;
88 std::vector<Influence> _stack;
89 _stack.push_back(Influence(sa(is), 0.0, 0.0, w));
95 Influence & s = _stack.back();
96 double diff = current - s.center;
97 double intersection = current + (sa(is) - s.prevVal - sigma2*
sq(diff)) / (sigma22 * diff);
99 if( intersection < s.left)
104 _stack.push_back(Influence(sa(is), 0.0, current, w));
111 else if(intersection < s.right)
113 s.right = intersection;
114 _stack.push_back(Influence(sa(is), intersection, current, w));
123 typename std::vector<Influence>::iterator it = _stack.begin();
124 for(current = 0.0; current < w; ++current, ++id)
126 while( current >= it->right)
128 da.set(sigma2 *
sq(current - it->center) + it->prevVal,
id);
132 template <
class SrcIterator,
class SrcAccessor,
133 class DestIterator,
class DestAccessor>
134 inline void distParabola(triple<SrcIterator, SrcIterator, SrcAccessor> src,
135 pair<DestIterator, DestAccessor> dest,
double sigma)
137 distParabola(src.first, src.second, src.third,
138 dest.first, dest.second, sigma);
147 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
148 class DestIterator,
class DestAccessor,
class Array>
149 void internalSeparableMultiArrayDistTmp(
150 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
151 DestIterator di, DestAccessor dest, Array
const & sigmas,
bool invert)
156 enum { N = SrcShape::static_size};
159 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
162 ArrayVector<TmpType> tmp( shape[0] );
164 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
165 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
169 SNavigator snav( si, shape, 0 );
170 DNavigator dnav( di, shape, 0 );
172 using namespace vigra::functor;
174 for( ; snav.hasMore(); snav++, dnav++ )
179 transformLine( snav.begin(), snav.end(), src, tmp.begin(),
180 typename AccessorTraits<TmpType>::default_accessor(),
181 Param(NumericTraits<TmpType>::zero())-Arg1());
183 copyLine( snav.begin(), snav.end(), src, tmp.begin(),
184 typename AccessorTraits<TmpType>::default_accessor() );
186 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
187 typename AccessorTraits<TmpType>::default_const_accessor()),
188 destIter( dnav.begin(), dest ), sigmas[0] );
192 for(
int d = 1; d < N; ++d )
194 DNavigator dnav( di, shape, d );
196 tmp.resize( shape[d] );
199 for( ; dnav.hasMore(); dnav++ )
202 copyLine( dnav.begin(), dnav.end(), dest,
203 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
205 detail::distParabola( srcIterRange(tmp.begin(), tmp.end(),
206 typename AccessorTraits<TmpType>::default_const_accessor()),
207 destIter( dnav.begin(), dest ), sigmas[d] );
213 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
214 class DestIterator,
class DestAccessor,
class Array>
215 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape
const & shape, SrcAccessor src,
216 DestIterator di, DestAccessor dest, Array
const & sigmas)
218 internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas,
false );
221 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
222 class DestIterator,
class DestAccessor>
223 inline void internalSeparableMultiArrayDistTmp( SrcIterator si, SrcShape
const & shape, SrcAccessor src,
224 DestIterator di, DestAccessor dest)
226 ArrayVector<double> sigmas(shape.size(), 1.0);
227 internalSeparableMultiArrayDistTmp( si, shape, src, di, dest, sigmas,
false );
339 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
340 class DestIterator,
class DestAccessor,
class Array>
342 DestIterator d, DestAccessor dest,
bool background,
343 Array
const & pixelPitch)
345 int N = shape.size();
347 typedef typename SrcAccessor::value_type SrcType;
348 typedef typename DestAccessor::value_type DestType;
349 typedef typename NumericTraits<DestType>::RealPromote Real;
351 SrcType zero = NumericTraits<SrcType>::zero();
354 bool pixelPitchIsReal =
false;
355 for(
int k=0; k<N; ++k)
357 if(
int(pixelPitch[k]) != pixelPitch[k])
358 pixelPitchIsReal =
true;
359 dmax +=
sq(pixelPitch[k]*shape[k]);
362 using namespace vigra::functor;
364 if(dmax > NumericTraits<DestType>::toRealPromote(NumericTraits<DestType>::max())
368 Real maxDist = (Real)dmax, rzero = (Real)0.0;
369 MultiArray<SrcShape::static_size, Real> tmpArray(shape);
370 if(background ==
true)
372 tmpArray.traverser_begin(),
typename AccessorTraits<Real>::default_accessor(),
373 ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
376 tmpArray.traverser_begin(),
typename AccessorTraits<Real>::default_accessor(),
377 ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
379 detail::internalSeparableMultiArrayDistTmp( tmpArray.traverser_begin(),
380 shape,
typename AccessorTraits<Real>::default_accessor(),
381 tmpArray.traverser_begin(),
382 typename AccessorTraits<Real>::default_accessor(), pixelPitch);
389 DestType maxDist = DestType(
std::ceil(dmax)), rzero = (DestType)0;
390 if(background ==
true)
392 ifThenElse( Arg1() == Param(zero), Param(maxDist), Param(rzero) ));
395 ifThenElse( Arg1() != Param(zero), Param(maxDist), Param(rzero) ));
397 detail::internalSeparableMultiArrayDistTmp( d, shape, dest, d, dest, pixelPitch);
401 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
402 class DestIterator,
class DestAccessor,
class Array>
404 pair<DestIterator, DestAccessor>
const & dest,
bool background,
405 Array
const & pixelPitch)
408 dest.first, dest.second, background, pixelPitch );
411 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
412 class DestIterator,
class DestAccessor>
415 DestIterator d, DestAccessor dest,
bool background)
417 ArrayVector<double> pixelPitch(shape.size(), 1.0);
421 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
422 class DestIterator,
class DestAccessor>
424 pair<DestIterator, DestAccessor>
const & dest,
bool background)
427 dest.first, dest.second, background );
509 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
510 class DestIterator,
class DestAccessor,
class Array>
512 DestIterator d, DestAccessor dest,
bool background,
513 Array
const & pixelPitch)
518 using namespace vigra::functor;
523 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
524 class DestIterator,
class DestAccessor>
526 DestIterator d, DestAccessor dest,
bool background)
531 using namespace vigra::functor;
536 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
537 class DestIterator,
class DestAccessor,
class Array>
539 pair<DestIterator, DestAccessor>
const & dest,
bool background,
540 Array
const & pixelPitch)
543 dest.first, dest.second, background, pixelPitch );
546 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
547 class DestIterator,
class DestAccessor>
549 pair<DestIterator, DestAccessor>
const & dest,
bool background)
552 dest.first, dest.second, background );
560 #endif //-- VIGRA_MULTI_DISTANCE_HXX