37 #ifndef VIGRA_VECTOR_DISTANCE_HXX
38 #define VIGRA_VECTOR_DISTANCE_HXX
43 #include "array_vector.hxx"
44 #include "multi_array.hxx"
45 #include "accessor.hxx"
46 #include "numerictraits.hxx"
47 #include "navigator.hxx"
48 #include "metaprogramming.hxx"
49 #include "multi_pointoperators.hxx"
50 #include "functorexpression.hxx"
51 #include "multi_distance.hxx"
53 #undef VECTORIAL_DIST_DEBUG
61 template <
class Vector,
class Value>
62 struct VectorialDistParabolaStackEntry
64 double left, center, right;
68 VectorialDistParabolaStackEntry(
const Vector& vec, Value prev,
double l,
double c,
double r)
69 : left(l), center(c), right(r), apex_height(prev), point(vec)
73 #ifdef VECTORIAL_DIST_DEBUG
74 template <
class Vector,
class Value>
75 std::ostream& operator<<(std::ostream&o, const VectorialDistParabolaStackEntry<Vector, Value>& e) {
76 o <<
"l=" << e.left <<
", c=" << e.center <<
", r=" << e.right <<
", pV=" << e.apex_height <<
", pVec=" << e.point;
87 template <
class VEC,
class ARRAY>
89 partialSquaredMagnitude(
const VEC& vec,
MultiArrayIndex dim, ARRAY
const & pixel_pitch)
96 sqMag +=
sq(pixel_pitch[i]*vec[i]);
101 template <
class SrcIterator,
105 SrcIterator is, SrcIterator iend,
106 Array
const & pixel_pitch )
108 typedef typename SrcIterator::value_type SrcType;
109 typedef VectorialDistParabolaStackEntry<SrcType, double> Influence;
111 double sigma = pixel_pitch[dimension],
113 double w = iend - is;
117 std::vector<Influence> _stack;
118 double apex_height = partialSquaredMagnitude(*is, dimension, pixel_pitch);
119 _stack.push_back(Influence(*is, apex_height, 0.0, 0.0, w));
121 double current = 1.0;
124 apex_height = partialSquaredMagnitude(*is, dimension, pixel_pitch);
125 Influence & s = _stack.back();
126 double diff = current - s.center;
127 double intersection = current + (apex_height - s.apex_height -
sq(sigma*diff)) / (2.0*sigma2 * diff);
129 if( intersection < s.left)
133 _stack.push_back(Influence(*is, apex_height, 0.0, current, w));
137 else if(intersection < s.right)
139 s.right = intersection;
140 _stack.push_back(Influence(*is, apex_height, intersection, current, w));
149 typename std::vector<Influence>::iterator it = _stack.begin();
150 for(current = 0.0; current < w; ++current, ++id)
152 while( current >= it->right)
156 (*id)[dimension] = it->center - current;
160 template <
class DestIterator,
162 class Array1,
class Array2>
165 DestIterator is, DestIterator iend,
166 LabelIterator ilabels,
167 Array1
const & pixel_pitch,
169 bool array_border_is_active=
false)
171 double w = iend - is;
175 typedef typename LabelIterator::value_type LabelType;
176 typedef typename DestIterator::value_type DestType;
177 typedef VectorialDistParabolaStackEntry<DestType, double> Influence;
178 typedef std::vector<Influence> Stack;
180 DestIterator
id = is;
181 DestType border_point = array_border_is_active
184 double apex_height = partialSquaredMagnitude(border_point, dimension, pixel_pitch);
185 Stack _stack(1, Influence(border_point, apex_height, 0.0, -1.0, w));
186 LabelType current_label = *ilabels;
187 for(
double begin = 0.0, current = 0.0; current <= w; ++ilabels, ++is, ++current)
189 DestType point = (current < w)
190 ? (current_label == *ilabels)
194 apex_height = partialSquaredMagnitude(point, dimension, pixel_pitch);
197 Influence & s = _stack.back();
198 double diff = (current - s.center)*pixel_pitch[dimension];
199 double intersection = current + (apex_height - s.apex_height -
sq(diff)) / (2.0 * diff);
201 if(intersection < s.left)
205 intersection = begin;
209 else if(intersection < s.right)
211 s.right = intersection;
214 _stack.push_back(Influence(point, apex_height, intersection, current, w));
215 if(current < w && current_label == *ilabels)
219 typename Stack::iterator it = _stack.begin();
220 for(
double c = begin; c < current; ++c, ++id)
222 while(c >= it->right)
225 (*id)[dimension] = it->center - c;
232 current_label = *ilabels;
234 apex_height = partialSquaredMagnitude(point, dimension, pixel_pitch);
235 Stack(1, Influence(DestType(0), 0.0, begin-1.0, begin-1.0, w)).swap(_stack);
242 template <
unsigned int N,
class T1,
class S1,
246 interpixelBoundaryVectorDistance(MultiArrayView<N, T1, S1>
const & labels,
247 MultiArrayView<N, T2, S2> dest,
248 Array
const & pixelPitch)
250 typedef typename MultiArrayShape<N>::type Shape;
251 typedef GridGraph<N> Graph;
252 typedef typename Graph::Node Node;
253 typedef typename Graph::NodeIt graph_scanner;
254 typedef typename Graph::OutArcIt neighbor_iterator;
256 Graph g(labels.shape());
257 for (graph_scanner node(g); node != lemon_graph::INVALID; ++node)
259 T1 label = labels[*node];
260 double min_dist = NumericTraits<double>::max();
262 boundary = point + Node(dest[point]),
263 min_pos = lemon::INVALID;
267 if(labels.isInside(boundary))
269 for (neighbor_iterator arc(g, boundary); arc != lemon_graph::INVALID; ++arc)
271 if(label == labels[g.target(*arc)])
273 double dist =
squaredNorm(pixelPitch*(g.target(*arc) - point));
277 min_pos = g.target(*arc);
281 if(min_pos == lemon::INVALID)
283 min_dist = NumericTraits<double>::max();
287 min_pos =
clip(boundary, Shape(0), labels.shape()-Shape(1));
288 min_diff = 0.5*(boundary + min_pos) - point;
293 for (neighbor_iterator arc(g, min_pos); arc != lemon_graph::INVALID; ++arc)
295 if(label != labels[g.target(*arc)])
297 T2 diff = 0.5*(g.target(*arc) + min_pos) - point;
306 dest[point] = min_diff;
317 struct Error_output_pixel_type_must_be_TinyVector_of_appropriate_length
318 : vigra::staticAssert::AssertBool<PRED> {};
365 template <
unsigned int N,
class T1,
class S1,
366 class T2,
class S2,
class Array>
369 MultiArrayView<N, T2, S2> dest,
371 Array
const & pixelPitch)
373 using namespace vigra::functor;
374 typedef typename MultiArrayView<N, T2, S2>::traverser Traverser;
375 typedef MultiArrayNavigator<Traverser, N> Navigator;
377 VIGRA_STATIC_ASSERT((Error_output_pixel_type_must_be_TinyVector_of_appropriate_length<N == T2::static_size>));
378 vigra_precondition(source.shape() == dest.shape(),
379 "separableVectorDistance(): shape mismatch between input and output.");
380 vigra_precondition(pixelPitch.size() == N,
381 "separableVectorDistance(): pixelPitch has wrong length.");
383 T2 maxDist(2*
sum(source.shape()*pixelPitch)), rzero;
384 if(background ==
true)
386 ifThenElse( Arg1() == Param(0), Param(maxDist), Param(rzero) ));
389 ifThenElse( Arg1() != Param(0), Param(maxDist), Param(rzero) ));
391 for(
int d = 0; d < N; ++d )
393 Navigator nav( dest.traverser_begin(), dest.shape(), d);
394 for( ; nav.hasMore(); nav++ )
396 detail::vectorialDistParabola(d, nav.begin(), nav.end(), pixelPitch);
401 template <
unsigned int N,
class T1,
class S1,
405 MultiArrayView<N, T2, S2> dest,
406 bool background=
true)
408 TinyVector<double, N> pixelPitch(1.0);
458 template <
unsigned int N,
class T1,
class S1,
463 MultiArrayView<N, T2, S2> dest,
464 bool array_border_is_active,
466 Array
const & pixelPitch)
468 VIGRA_STATIC_ASSERT((Error_output_pixel_type_must_be_TinyVector_of_appropriate_length<N == T2::static_size>));
469 vigra_precondition(labels.shape() == dest.shape(),
470 "boundaryVectorDistance(): shape mismatch between input and output.");
471 vigra_precondition(pixelPitch.size() == N,
472 "boundaryVectorDistance(): pixelPitch has wrong length.");
474 using namespace vigra::functor;
478 MultiArray<N, unsigned char> boundaries(labels.shape());
481 if(array_border_is_active)
489 vigra_precondition(!NumericTraits<T2>::isIntegral::value,
490 "boundaryVectorDistance(..., InterpixelBoundary): output pixel type must be float or double.");
493 typedef typename MultiArrayView<N, T1, S1>::const_traverser LabelIterator;
494 typedef typename MultiArrayView<N, T2, S2>::traverser DestIterator;
495 typedef MultiArrayNavigator<LabelIterator, N> LabelNavigator;
496 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
498 T2 maxDist(2*
sum(labels.shape()*pixelPitch));
500 for(
int d = 0; d < N; ++d )
502 LabelNavigator lnav( labels.traverser_begin(), labels.shape(), d );
503 DNavigator dnav( dest.traverser_begin(), dest.shape(), d );
505 for( ; dnav.hasMore(); dnav++, lnav++ )
507 detail::boundaryVectorDistParabola(d, dnav.begin(), dnav.end(), lnav.begin(),
508 pixelPitch, maxDist, array_border_is_active);
514 detail::interpixelBoundaryVectorDistance(labels, dest, pixelPitch);
519 template <
unsigned int N,
class T1,
class S1,
523 MultiArrayView<N, T2, S2> dest,
524 bool array_border_is_active=
false,
527 TinyVector<double, N> pixelPitch(1.0);
535 #endif //-- VIGRA_VECTOR_DISTANCE_HXX