37 #ifndef VIGRA_MULTI_LOCALMINMAX_HXX
38 #define VIGRA_MULTI_LOCALMINMAX_HXX
42 #include "multi_array.hxx"
43 #include "localminmax.hxx"
50 template <
unsigned int M>
51 struct IsLocalExtremum2
53 template <
class T,
class Shape,
class Compare>
54 static bool exec(T * v, Shape
const & stride, Compare
const & compare)
56 return compare(*v, *(v-stride[M])) &&
57 compare(*v, *(v+stride[M])) &&
58 IsLocalExtremum2<M-1>::exec(v, stride, compare);
61 template <
class Shape,
class T,
class Compare>
62 static bool execAtBorder(Shape
const & point, Shape
const & shape,
63 T * v, Shape
const & stride, Compare
const & compare)
65 return (point[M] > 0 && compare(*v, *(v-stride[M]))) &&
66 (point[M] < shape[M]-1 && compare(*v, *(v+stride[M]))) &&
67 IsLocalExtremum2<M-1>::exec(point, shape, v, stride, compare);
72 struct IsLocalExtremum2<0>
74 template <
class T,
class Shape,
class Compare>
75 static bool exec(T * v, Shape
const & stride, Compare
const & compare)
77 return compare(*v, *(v-stride[0])) && compare(*v, *(v+stride[0]));
80 template <
class Shape,
class T,
class Compare>
81 static bool execAtBorder(Shape
const & point, Shape
const & shape,
82 T * v, Shape
const & stride, Compare
const & compare)
84 return (point[0] > 0 && compare(*v, *(v-stride[0]))) &&
85 (point[0] < shape[0]-1 && compare(*v, *(v+stride[0])));
90 template <
unsigned int M>
91 struct IsLocalExtremum3
93 template <
class T,
class Shape,
class Compare>
94 static bool exec(T * v, Shape
const & stride, Compare
const & compare)
96 return exec(v, v, stride, compare,
true);
99 template <
class T,
class Shape,
class Compare>
100 static bool exec(T * v, T * u, Shape
const & stride,
101 Compare
const & compare,
bool isCenter)
103 return IsLocalExtremum3<M-1>::exec(v, u-stride[M], stride, compare,
false) &&
104 IsLocalExtremum3<M-1>::exec(v, u, stride, compare, isCenter) &&
105 IsLocalExtremum3<M-1>::exec(v, u+stride[M], stride, compare,
false);
108 template <
class Shape,
class T,
class Compare>
109 static bool execAtBorder(Shape
const & point, Shape
const & shape,
110 T * v, Shape
const & stride, Compare
const & compare)
112 return execAtBorder(point, shape, v, v, stride, compare,
true);
115 template <
class Shape,
class T,
class Compare>
116 static bool execAtBorder(Shape
const & point, Shape
const & shape,
117 T * v, T * u, Shape
const & stride,
118 Compare
const & compare,
bool isCenter)
120 return (point[M] > 0 && IsLocalExtremum3<M-1>::exec(v, u-stride[M], stride, compare,
false)) &&
121 IsLocalExtremum3<M-1>::exec(point, shape, v, u, stride, compare, isCenter) &&
122 (point[M] < shape[M]-1 && IsLocalExtremum3<M-1>::exec(v, u+stride[M], stride, compare,
false));
127 struct IsLocalExtremum3<0>
129 template <
class T,
class Shape,
class Compare>
130 static bool exec(T * v, Shape
const & stride, Compare
const & compare)
132 return compare(*v, *(v-stride[0])) && compare(*v, *(v+stride[0]));
135 template <
class T,
class Shape,
class Compare>
136 static bool exec(T * v, T * u, Shape
const & stride,
137 Compare
const & compare,
bool isCenter)
139 return compare(*v, *(u-stride[0])) &&
140 (!isCenter && compare(*v, *u)) &&
141 compare(*v, *(u+stride[0]));
144 template <
class Shape,
class T,
class Compare>
145 static bool execAtBorder(Shape
const & point, Shape
const & shape,
146 T * v, Shape
const & stride, Compare
const & compare)
148 return (point[0] > 0 && compare(*v, *(v-stride[0]))) &&
149 (point[0] < shape[0]-1 && compare(*v, *(v+stride[0])));
152 template <
class Shape,
class T,
class Compare>
153 static bool execAtBorder(Shape
const & point, Shape
const & shape,
154 T * v, T * u, Shape
const & stride,
155 Compare
const & compare,
bool isCenter)
157 return (point[M] > 0 && compare(*v, *(u-stride[0]))) &&
158 (!isCenter && compare(*v, *u)) &&
159 (point[M] < shape[M]-1 && compare(*v, *(u+stride[0])));
163 template <
unsigned int N,
class T1,
class C1,
class T2,
class C2,
class Compare>
165 localMinMax(MultiArrayView<N, T1, C1> src,
166 MultiArrayView<N, T2, C2> dest,
167 T2 marker,
unsigned int neighborhood,
170 bool allowExtremaAtBorder =
false)
172 typedef typename MultiArrayShape<N>::type Shape;
173 typedef MultiCoordinateNavigator<N> Navigator;
175 Shape shape = src.shape(),
178 vigra_precondition(shape == dest.shape(),
179 "localMinMax(): Shape mismatch between input and output.");
180 vigra_precondition(neighborhood == 2*N || neighborhood == pow(3, N) - 1,
181 "localMinMax(): Invalid neighborhood.");
183 if(allowExtremaAtBorder)
185 for(
unsigned int d=0; d<N; ++d)
187 Navigator nav(shape, d);
188 for(; nav.hasMore(); ++nav)
190 Shape i = nav.begin();
192 for(; i[d] < shape[d]; i[d] += shape[d]-1)
194 if(!compare(src[i], threshold))
197 if(neighborhood == 2*N)
199 if(IsLocalExtremum2<N>::execAtBorder(i, shape, &src[i],
200 src.stride(), compare))
205 if(IsLocalExtremum3<N>::execAtBorder(i, shape, &src[i],
206 src.stride(), compare))
214 src = src.subarray(unit, shape - unit);
215 dest = dest.subarray(unit, shape - unit);
218 Navigator nav(shape, 0);
219 for(; nav.hasMore(); ++nav)
221 Shape i = nav.begin();
223 for(; i[0] < shape[0]; ++i[0])
225 if(!compare(src[i], threshold))
228 if(neighborhood == 2*N)
230 if(IsLocalExtremum2<N>::exec(&src[i], src.stride(), compare))
235 if(IsLocalExtremum3<N>::exec(&src[i], src.stride(), compare))
242 template <
class T1,
class C1,
class T2,
class C2,
243 class Neighborhood,
class Compare,
class Equal>
245 extendLocalMinMax(MultiArrayView<3, T1, C1> src,
246 MultiArrayView<3, T2, C2> dest,
248 Neighborhood neighborhood,
249 Compare compare, Equal equal,
251 bool allowExtremaAtBorder =
false)
253 typedef typename MultiArrayView<3, T1, C1>::traverser SrcIterator;
254 typedef typename MultiArrayView<3, T2, C2>::traverser DestIterator;
255 typedef typename MultiArray<3, int>::traverser LabelIterator;
256 typedef MultiArrayShape<3>::type Shape;
258 Shape shape = src.shape();
261 vigra_precondition(shape == dest.shape(),
262 "extendLocalMinMax(): Shape mismatch between input and output.");
264 MultiArray<3, int> labels(shape);
266 int number_of_regions =
labelVolume(srcMultiArrayRange(src), destMultiArray(labels),
267 neighborhood, equal);
270 ArrayVector<unsigned char> isExtremum(number_of_regions+1, (
unsigned char)1);
272 SrcIterator zs = src.traverser_begin();
273 LabelIterator zl = labels.traverser_begin();
277 LabelIterator yl(zl);
282 LabelIterator xl(yl);
289 if(isExtremum[lab] == 0)
292 if(!compare(v, threshold))
302 NeighborhoodCirculator<SrcIterator, Neighborhood> cs(xs);
303 NeighborhoodCirculator<LabelIterator, Neighborhood> cl(xl);
304 for(i=0; i<Neighborhood::DirectionCount; ++i, ++cs, ++cl)
306 if(lab != *cl && compare(*cs,v))
315 if(allowExtremaAtBorder)
317 RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood>
318 cs(xs, atBorder), scend(cs);
321 if(lab != *(xl+cs.diff()) && compare(cs,v))
327 while(++cs != scend);
339 zl = labels.traverser_begin();
340 DestIterator zd = dest.traverser_begin();
343 LabelIterator yl(zl);
348 LabelIterator xl(yl);
369 template <
unsigned int N,
class T1,
class C1,
class T2,
class C2>
372 MultiArrayView<N, T2, C2> dest,
373 LocalMinmaxOptions
const & options = LocalMinmaxOptions())
375 T1 threshold = options.use_threshold
376 ? std::min(NumericTraits<T1>::max(), (T1)options.thresh)
377 : NumericTraits<T1>::max();
378 T2 marker = (T2)options.marker;
380 vigra_precondition(!options.allow_plateaus,
381 "localMinima(): Option 'allowPlateaus' is not implemented for arbitrary dimensions,\n"
382 " use extendedLocalMinima() for 2D and 3D problems.");
384 if(options.neigh == 0)
386 if(options.neigh == 1)
387 options.neigh = pow(3, N) - 1;
389 detail::localMinMax(src, dest, marker, options.neigh,
390 threshold, std::less<T1>(), options.allow_at_border);
400 template <
unsigned int N,
class T1,
class C1,
class T2,
class C2>
403 MultiArrayView<N, T2, C2> dest,
404 LocalMinmaxOptions
const & options = LocalMinmaxOptions())
406 T1 threshold = options.use_threshold
407 ? std::max(NumericTraits<T1>::min(), (T1)options.thresh)
408 : NumericTraits<T1>::min();
409 T2 marker = (T2)options.marker;
411 vigra_precondition(!options.allow_plateaus,
412 "localMaxima(): Option 'allowPlateaus' is not implemented for arbitrary dimensions,\n"
413 " use extendedLocalMinima() for 2D and 3D problems.");
415 if(options.neigh == 0)
417 if(options.neigh == 1)
418 options.neigh = pow(3, N) - 1;
420 detail::localMinMax(src, dest, marker, options.neigh,
421 threshold, std::greater<T1>(), options.allow_at_border);
433 template <
class T1,
class C1,
class T2,
class C2,
434 class Neighborhood,
class EqualityFunctor>
437 MultiArrayView<3, T2, C2> dest,
438 LocalMinmaxOptions
const & options = LocalMinmaxOptions())
440 T1 threshold = options.use_threshold
441 ? std::min(NumericTraits<T1>::max(), (T1)options.thresh)
442 : NumericTraits<T1>::max();
443 T2 marker = (T2)options.marker;
445 if(options.neigh == 0 || options.neigh == 6)
448 threshold, std::less<T1>(), std::equal_to<T1>(),
449 options.allow_at_border);
451 else if(options.neigh == 1 || options.neigh == 26)
454 threshold, std::less<T1>(), std::equal_to<T1>(),
455 options.allow_at_border);
458 vigra_precondition(
false,
459 "extendedLocalMinima(): Invalid neighborhood.");
469 template <
class T1,
class C1,
class T2,
class C2,
470 class Neighborhood,
class EqualityFunctor>
473 MultiArrayView<3, T2, C2> dest,
474 LocalMinmaxOptions
const & options = LocalMinmaxOptions())
476 T1 threshold = options.use_threshold
477 ? std::max(NumericTraits<T1>::min(), (T1)options.thresh)
478 : NumericTraits<T1>::min();
479 T2 marker = (T2)options.marker;
481 if(options.neigh == 0 || options.neigh == 6)
484 threshold, std::greater<T1>(), std::equal_to<T1>(),
485 options.allow_at_border);
487 else if(options.neigh == 1 || options.neigh == 26)
490 threshold, std::greater<T1>(), std::equal_to<T1>(),
491 options.allow_at_border);
494 vigra_precondition(
false,
495 "extendedLocalMaxima(): Invalid neighborhood.");
502 #endif // VIGRA_MULTI_LOCALMINMAX_HXX
void extendedLocalMinima(...)
Find local minimal regions in an image or volume.
unsigned int labelVolume(...)
Find the connected components of a segmented volume.
AtImageBorder AtVolumeBorder
Encode whether a voxel is near the volume border.
Definition: voxelneighborhood.hxx:72
Neighborhood3DTwentySix::NeighborCode3D NeighborCode3DTwentySix
Definition: voxelneighborhood.hxx:1602
void localMinima(...)
Find local minima in an image or multi-dimensional array.
AtVolumeBorder isAtVolumeBorder(int x, int y, int z, int width, int height, int depth)
Find out whether a voxel is at the volume border.
Definition: voxelneighborhood.hxx:82
std::ptrdiff_t MultiArrayIndex
Definition: multi_iterator.hxx:348
void localMaxima(...)
Find local maxima in an image or multi-dimensional array.
Neighborhood3DSix::NeighborCode3D NeighborCode3DSix
Definition: voxelneighborhood.hxx:468
void extendedLocalMaxima(...)
Find local maximal regions in an image or volume.
Definition: pixelneighborhood.hxx:70