[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

multi_localminmax.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2010 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_MULTI_LOCALMINMAX_HXX
38 #define VIGRA_MULTI_LOCALMINMAX_HXX
39 
40 #include <vector>
41 #include <functional>
42 #include "multi_array.hxx"
43 #include "localminmax.hxx"
44 #if 0
45 namespace vigra {
46 
47 namespace detail {
48 
49 // direct neighborhood
50 template <unsigned int M>
51 struct IsLocalExtremum2
52 {
53  template <class T, class Shape, class Compare>
54  static bool exec(T * v, Shape const & stride, Compare const & compare)
55  {
56  return compare(*v, *(v-stride[M])) &&
57  compare(*v, *(v+stride[M])) &&
58  IsLocalExtremum2<M-1>::exec(v, stride, compare);
59  }
60 
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)
64  {
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);
68  }
69 };
70 
71 template <>
72 struct IsLocalExtremum2<0>
73 {
74  template <class T, class Shape, class Compare>
75  static bool exec(T * v, Shape const & stride, Compare const & compare)
76  {
77  return compare(*v, *(v-stride[0])) && compare(*v, *(v+stride[0]));
78  }
79 
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)
83  {
84  return (point[0] > 0 && compare(*v, *(v-stride[0]))) &&
85  (point[0] < shape[0]-1 && compare(*v, *(v+stride[0])));
86  }
87 };
88 
89 // indirect neighborhood
90 template <unsigned int M>
91 struct IsLocalExtremum3
92 {
93  template <class T, class Shape, class Compare>
94  static bool exec(T * v, Shape const & stride, Compare const & compare)
95  {
96  return exec(v, v, stride, compare, true);
97  }
98 
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)
102  {
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);
106  }
107 
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)
111  {
112  return execAtBorder(point, shape, v, v, stride, compare, true);
113  }
114 
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)
119  {
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));
123  }
124 };
125 
126 template <>
127 struct IsLocalExtremum3<0>
128 {
129  template <class T, class Shape, class Compare>
130  static bool exec(T * v, Shape const & stride, Compare const & compare)
131  {
132  return compare(*v, *(v-stride[0])) && compare(*v, *(v+stride[0]));
133  }
134 
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)
138  {
139  return compare(*v, *(u-stride[0])) &&
140  (!isCenter && compare(*v, *u)) &&
141  compare(*v, *(u+stride[0]));
142  }
143 
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)
147  {
148  return (point[0] > 0 && compare(*v, *(v-stride[0]))) &&
149  (point[0] < shape[0]-1 && compare(*v, *(v+stride[0])));
150  }
151 
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)
156  {
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])));
160  }
161 };
162 
163 template <unsigned int N, class T1, class C1, class T2, class C2, class Compare>
164 void
165 localMinMax(MultiArrayView<N, T1, C1> src,
166  MultiArrayView<N, T2, C2> dest,
167  T2 marker, unsigned int neighborhood,
168  T1 threshold,
169  Compare compare,
170  bool allowExtremaAtBorder = false)
171 {
172  typedef typename MultiArrayShape<N>::type Shape;
173  typedef MultiCoordinateNavigator<N> Navigator;
174 
175  Shape shape = src.shape(),
176  unit = Shape(MultiArrayIndex(1));
177 
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.");
182 
183  if(allowExtremaAtBorder)
184  {
185  for(unsigned int d=0; d<N; ++d)
186  {
187  Navigator nav(shape, d);
188  for(; nav.hasMore(); ++nav)
189  {
190  Shape i = nav.begin();
191 
192  for(; i[d] < shape[d]; i[d] += shape[d]-1)
193  {
194  if(!compare(src[i], threshold))
195  continue;
196 
197  if(neighborhood == 2*N)
198  {
199  if(IsLocalExtremum2<N>::execAtBorder(i, shape, &src[i],
200  src.stride(), compare))
201  dest[i] = marker;
202  }
203  else
204  {
205  if(IsLocalExtremum3<N>::execAtBorder(i, shape, &src[i],
206  src.stride(), compare))
207  dest[i] = marker;
208  }
209  }
210  }
211  }
212  }
213 
214  src = src.subarray(unit, shape - unit);
215  dest = dest.subarray(unit, shape - unit);
216  shape = src.shape();
217 
218  Navigator nav(shape, 0);
219  for(; nav.hasMore(); ++nav)
220  {
221  Shape i = nav.begin();
222 
223  for(; i[0] < shape[0]; ++i[0])
224  {
225  if(!compare(src[i], threshold))
226  continue;
227 
228  if(neighborhood == 2*N)
229  {
230  if(IsLocalExtremum2<N>::exec(&src[i], src.stride(), compare))
231  dest[i] = marker;
232  }
233  else
234  {
235  if(IsLocalExtremum3<N>::exec(&src[i], src.stride(), compare))
236  dest[i] = marker;
237  }
238  }
239  }
240 }
241 
242 template <class T1, class C1, class T2, class C2,
243  class Neighborhood, class Compare, class Equal>
244 void
245 extendLocalMinMax(MultiArrayView<3, T1, C1> src,
246  MultiArrayView<3, T2, C2> dest,
247  T2 marker,
248  Neighborhood neighborhood,
249  Compare compare, Equal equal,
250  T1 threshold,
251  bool allowExtremaAtBorder = false)
252 {
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;
257 
258  Shape shape = src.shape();
259  MultiArrayIndex w = shape[0], h = shape[1], d = shape[2];
260 
261  vigra_precondition(shape == dest.shape(),
262  "extendLocalMinMax(): Shape mismatch between input and output.");
263 
264  MultiArray<3, int> labels(shape);
265 
266  int number_of_regions = labelVolume(srcMultiArrayRange(src), destMultiArray(labels),
267  neighborhood, equal);
268 
269  // assume that a region is a extremum until the opposite is proved
270  ArrayVector<unsigned char> isExtremum(number_of_regions+1, (unsigned char)1);
271 
272  SrcIterator zs = src.traverser_begin();
273  LabelIterator zl = labels.traverser_begin();
274  for(MultiArrayIndex z = 0; z != d; ++z, ++zs.dim2(), ++zl.dim2())
275  {
276  SrcIterator ys(zs);
277  LabelIterator yl(zl);
278 
279  for(MultiArrayIndex y = 0; y != h; ++y, ++ys.dim1(), ++yl.dim1())
280  {
281  SrcIterator xs(ys);
282  LabelIterator xl(yl);
283 
284  for(MultiArrayIndex x = 0; x != w; ++x, ++xs.dim0(), ++xl.dim0())
285  {
286  int lab = *xl;
287  T1 v = *xs;
288 
289  if(isExtremum[lab] == 0)
290  continue;
291 
292  if(!compare(v, threshold))
293  {
294  // mark all regions that don't exceed the threshold as non-extremum
295  isExtremum[lab] = 0;
296  continue;
297  }
298 
299  AtVolumeBorder atBorder = isAtVolumeBorder(x, y, z, w, h, d);
300  if(atBorder == NotAtBorder)
301  {
302  NeighborhoodCirculator<SrcIterator, Neighborhood> cs(xs);
303  NeighborhoodCirculator<LabelIterator, Neighborhood> cl(xl);
304  for(i=0; i<Neighborhood::DirectionCount; ++i, ++cs, ++cl)
305  {
306  if(lab != *cl && compare(*cs,v))
307  {
308  isExtremum[lab] = 0;
309  break;
310  }
311  }
312  }
313  else
314  {
315  if(allowExtremaAtBorder)
316  {
317  RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood>
318  cs(xs, atBorder), scend(cs);
319  do
320  {
321  if(lab != *(xl+cs.diff()) && compare(cs,v))
322  {
323  isExtremum[lab] = 0;
324  break;
325  }
326  }
327  while(++cs != scend);
328  }
329  else
330  {
331  isExtremum[lab] = 0;
332  }
333  }
334  }
335  }
336  }
337 
338 
339  zl = labels.traverser_begin();
340  DestIterator zd = dest.traverser_begin();
341  for(MultiArrayIndex z = 0; z != d; ++z, ++zl.dim2(), ++zd.dim2())
342  {
343  LabelIterator yl(zl);
344  DestIterator yd(zd);
345 
346  for(MultiArrayIndex y = 0; y != h; ++y, ++yl.dim1(), ++yd.dim1())
347  {
348  LabelIterator xl(yl);
349  DestIterator xd(yd);
350 
351  for(MultiArrayIndex x = 0; x != w; ++x, ++xl.dim0(), ++xd.dim0())
352  {
353  if(isExtremum[*xl])
354  *xd = marker;
355  }
356  }
357  }
358 }
359 
360 } // namespace detail
361 
362 /********************************************************/
363 /* */
364 /* localMinima */
365 /* */
366 /********************************************************/
367 
368 // documentation is in localminmax.hxx
369 template <unsigned int N, class T1, class C1, class T2, class C2>
370 void
371 localMinima(MultiArrayView<N, T1, C1> src,
372  MultiArrayView<N, T2, C2> dest,
373  LocalMinmaxOptions const & options = LocalMinmaxOptions())
374 {
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;
379 
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.");
383 
384  if(options.neigh == 0)
385  options.neigh = 2*N;
386  if(options.neigh == 1)
387  options.neigh = pow(3, N) - 1;
388 
389  detail::localMinMax(src, dest, marker, options.neigh,
390  threshold, std::less<T1>(), options.allow_at_border);
391 }
392 
393 /********************************************************/
394 /* */
395 /* localMaxima */
396 /* */
397 /********************************************************/
398 
399 // documentation is in localminmax.hxx
400 template <unsigned int N, class T1, class C1, class T2, class C2>
401 void
402 localMaxima(MultiArrayView<N, T1, C1> src,
403  MultiArrayView<N, T2, C2> dest,
404  LocalMinmaxOptions const & options = LocalMinmaxOptions())
405 {
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;
410 
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.");
414 
415  if(options.neigh == 0)
416  options.neigh = 2*N;
417  if(options.neigh == 1)
418  options.neigh = pow(3, N) - 1;
419 
420  detail::localMinMax(src, dest, marker, options.neigh,
421  threshold, std::greater<T1>(), options.allow_at_border);
422 }
423 
424 /**************************************************************************/
425 
426 /********************************************************/
427 /* */
428 /* extendedLocalMinima */
429 /* */
430 /********************************************************/
431 
432 // documentation is in localminmax.hxx
433 template <class T1, class C1, class T2, class C2,
434  class Neighborhood, class EqualityFunctor>
435 inline void
436 extendedLocalMinima(MultiArrayView<3, T1, C1> src,
437  MultiArrayView<3, T2, C2> dest,
438  LocalMinmaxOptions const & options = LocalMinmaxOptions())
439 {
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;
444 
445  if(options.neigh == 0 || options.neigh == 6)
446  {
447  detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DSix(),
448  threshold, std::less<T1>(), std::equal_to<T1>(),
449  options.allow_at_border);
450  }
451  else if(options.neigh == 1 || options.neigh == 26)
452  {
453  detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DTwentySix(),
454  threshold, std::less<T1>(), std::equal_to<T1>(),
455  options.allow_at_border);
456  }
457  else
458  vigra_precondition(false,
459  "extendedLocalMinima(): Invalid neighborhood.");
460 }
461 
462 /********************************************************/
463 /* */
464 /* extendedLocalMaxima */
465 /* */
466 /********************************************************/
467 
468 // documentation is in localminmax.hxx
469 template <class T1, class C1, class T2, class C2,
470  class Neighborhood, class EqualityFunctor>
471 inline void
472 extendedLocalMaxima(MultiArrayView<3, T1, C1> src,
473  MultiArrayView<3, T2, C2> dest,
474  LocalMinmaxOptions const & options = LocalMinmaxOptions())
475 {
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;
480 
481  if(options.neigh == 0 || options.neigh == 6)
482  {
483  detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DSix(),
484  threshold, std::greater<T1>(), std::equal_to<T1>(),
485  options.allow_at_border);
486  }
487  else if(options.neigh == 1 || options.neigh == 26)
488  {
489  detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DTwentySix(),
490  threshold, std::greater<T1>(), std::equal_to<T1>(),
491  options.allow_at_border);
492  }
493  else
494  vigra_precondition(false,
495  "extendedLocalMaxima(): Invalid neighborhood.");
496 }
497 
498 } // namespace vigra
499 
500 #endif
501 
502 #endif // VIGRA_MULTI_LOCALMINMAX_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.9.0 (Tue Oct 22 2013)