[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/multi_localminmax.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2010 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 00037 #ifndef VIGRA_MULTI_LOCALMINMAX_HXX 00038 #define VIGRA_MULTI_LOCALMINMAX_HXX 00039 00040 #include <vector> 00041 #include <functional> 00042 #include "multi_array.hxx" 00043 #include "localminmax.hxx" 00044 #if 0 00045 namespace vigra { 00046 00047 namespace detail { 00048 00049 // direct neighborhood 00050 template <unsigned int M> 00051 struct IsLocalExtremum2 00052 { 00053 template <class T, class Shape, class Compare> 00054 static bool exec(T * v, Shape const & stride, Compare const & compare) 00055 { 00056 return compare(*v, *(v-stride[M])) && 00057 compare(*v, *(v+stride[M])) && 00058 IsLocalExtremum2<M-1>::exec(v, stride, compare); 00059 } 00060 00061 template <class Shape, class T, class Compare> 00062 static bool execAtBorder(Shape const & point, Shape const & shape, 00063 T * v, Shape const & stride, Compare const & compare) 00064 { 00065 return (point[M] > 0 && compare(*v, *(v-stride[M]))) && 00066 (point[M] < shape[M]-1 && compare(*v, *(v+stride[M]))) && 00067 IsLocalExtremum2<M-1>::exec(point, shape, v, stride, compare); 00068 } 00069 }; 00070 00071 template <> 00072 struct IsLocalExtremum2<0> 00073 { 00074 template <class T, class Shape, class Compare> 00075 static bool exec(T * v, Shape const & stride, Compare const & compare) 00076 { 00077 return compare(*v, *(v-stride[0])) && compare(*v, *(v+stride[0])); 00078 } 00079 00080 template <class Shape, class T, class Compare> 00081 static bool execAtBorder(Shape const & point, Shape const & shape, 00082 T * v, Shape const & stride, Compare const & compare) 00083 { 00084 return (point[0] > 0 && compare(*v, *(v-stride[0]))) && 00085 (point[0] < shape[0]-1 && compare(*v, *(v+stride[0]))); 00086 } 00087 }; 00088 00089 // indirect neighborhood 00090 template <unsigned int M> 00091 struct IsLocalExtremum3 00092 { 00093 template <class T, class Shape, class Compare> 00094 static bool exec(T * v, Shape const & stride, Compare const & compare) 00095 { 00096 return exec(v, v, stride, compare, true); 00097 } 00098 00099 template <class T, class Shape, class Compare> 00100 static bool exec(T * v, T * u, Shape const & stride, 00101 Compare const & compare, bool isCenter) 00102 { 00103 return IsLocalExtremum3<M-1>::exec(v, u-stride[M], stride, compare, false) && 00104 IsLocalExtremum3<M-1>::exec(v, u, stride, compare, isCenter) && 00105 IsLocalExtremum3<M-1>::exec(v, u+stride[M], stride, compare, false); 00106 } 00107 00108 template <class Shape, class T, class Compare> 00109 static bool execAtBorder(Shape const & point, Shape const & shape, 00110 T * v, Shape const & stride, Compare const & compare) 00111 { 00112 return execAtBorder(point, shape, v, v, stride, compare, true); 00113 } 00114 00115 template <class Shape, class T, class Compare> 00116 static bool execAtBorder(Shape const & point, Shape const & shape, 00117 T * v, T * u, Shape const & stride, 00118 Compare const & compare, bool isCenter) 00119 { 00120 return (point[M] > 0 && IsLocalExtremum3<M-1>::exec(v, u-stride[M], stride, compare, false)) && 00121 IsLocalExtremum3<M-1>::exec(point, shape, v, u, stride, compare, isCenter) && 00122 (point[M] < shape[M]-1 && IsLocalExtremum3<M-1>::exec(v, u+stride[M], stride, compare, false)); 00123 } 00124 }; 00125 00126 template <> 00127 struct IsLocalExtremum3<0> 00128 { 00129 template <class T, class Shape, class Compare> 00130 static bool exec(T * v, Shape const & stride, Compare const & compare) 00131 { 00132 return compare(*v, *(v-stride[0])) && compare(*v, *(v+stride[0])); 00133 } 00134 00135 template <class T, class Shape, class Compare> 00136 static bool exec(T * v, T * u, Shape const & stride, 00137 Compare const & compare, bool isCenter) 00138 { 00139 return compare(*v, *(u-stride[0])) && 00140 (!isCenter && compare(*v, *u)) && 00141 compare(*v, *(u+stride[0])); 00142 } 00143 00144 template <class Shape, class T, class Compare> 00145 static bool execAtBorder(Shape const & point, Shape const & shape, 00146 T * v, Shape const & stride, Compare const & compare) 00147 { 00148 return (point[0] > 0 && compare(*v, *(v-stride[0]))) && 00149 (point[0] < shape[0]-1 && compare(*v, *(v+stride[0]))); 00150 } 00151 00152 template <class Shape, class T, class Compare> 00153 static bool execAtBorder(Shape const & point, Shape const & shape, 00154 T * v, T * u, Shape const & stride, 00155 Compare const & compare, bool isCenter) 00156 { 00157 return (point[M] > 0 && compare(*v, *(u-stride[0]))) && 00158 (!isCenter && compare(*v, *u)) && 00159 (point[M] < shape[M]-1 && compare(*v, *(u+stride[0]))); 00160 } 00161 }; 00162 00163 template <unsigned int N, class T1, class C1, class T2, class C2, class Compare> 00164 void 00165 localMinMax(MultiArrayView<N, T1, C1> src, 00166 MultiArrayView<N, T2, C2> dest, 00167 T2 marker, unsigned int neighborhood, 00168 T1 threshold, 00169 Compare compare, 00170 bool allowExtremaAtBorder = false) 00171 { 00172 typedef typename MultiArrayShape<N>::type Shape; 00173 typedef MultiCoordinateNavigator<N> Navigator; 00174 00175 Shape shape = src.shape(), 00176 unit = Shape(MultiArrayIndex(1)); 00177 00178 vigra_precondition(shape == dest.shape(), 00179 "localMinMax(): Shape mismatch between input and output."); 00180 vigra_precondition(neighborhood == 2*N || neighborhood == pow(3, N) - 1, 00181 "localMinMax(): Invalid neighborhood."); 00182 00183 if(allowExtremaAtBorder) 00184 { 00185 for(unsigned int d=0; d<N; ++d) 00186 { 00187 Navigator nav(shape, d); 00188 for(; nav.hasMore(); ++nav) 00189 { 00190 Shape i = nav.begin(); 00191 00192 for(; i[d] < shape[d]; i[d] += shape[d]-1) 00193 { 00194 if(!compare(src[i], threshold)) 00195 continue; 00196 00197 if(neighborhood == 2*N) 00198 { 00199 if(IsLocalExtremum2<N>::execAtBorder(i, shape, &src[i], 00200 src.stride(), compare)) 00201 dest[i] = marker; 00202 } 00203 else 00204 { 00205 if(IsLocalExtremum3<N>::execAtBorder(i, shape, &src[i], 00206 src.stride(), compare)) 00207 dest[i] = marker; 00208 } 00209 } 00210 } 00211 } 00212 } 00213 00214 src = src.subarray(unit, shape - unit); 00215 dest = dest.subarray(unit, shape - unit); 00216 shape = src.shape(); 00217 00218 Navigator nav(shape, 0); 00219 for(; nav.hasMore(); ++nav) 00220 { 00221 Shape i = nav.begin(); 00222 00223 for(; i[0] < shape[0]; ++i[0]) 00224 { 00225 if(!compare(src[i], threshold)) 00226 continue; 00227 00228 if(neighborhood == 2*N) 00229 { 00230 if(IsLocalExtremum2<N>::exec(&src[i], src.stride(), compare)) 00231 dest[i] = marker; 00232 } 00233 else 00234 { 00235 if(IsLocalExtremum3<N>::exec(&src[i], src.stride(), compare)) 00236 dest[i] = marker; 00237 } 00238 } 00239 } 00240 } 00241 00242 template <class T1, class C1, class T2, class C2, 00243 class Neighborhood, class Compare, class Equal> 00244 void 00245 extendLocalMinMax(MultiArrayView<3, T1, C1> src, 00246 MultiArrayView<3, T2, C2> dest, 00247 T2 marker, 00248 Neighborhood neighborhood, 00249 Compare compare, Equal equal, 00250 T1 threshold, 00251 bool allowExtremaAtBorder = false) 00252 { 00253 typedef typename MultiArrayView<3, T1, C1>::traverser SrcIterator; 00254 typedef typename MultiArrayView<3, T2, C2>::traverser DestIterator; 00255 typedef typename MultiArray<3, int>::traverser LabelIterator; 00256 typedef MultiArrayShape<3>::type Shape; 00257 00258 Shape shape = src.shape(); 00259 MultiArrayIndex w = shape[0], h = shape[1], d = shape[2]; 00260 00261 vigra_precondition(shape == dest.shape(), 00262 "extendLocalMinMax(): Shape mismatch between input and output."); 00263 00264 MultiArray<3, int> labels(shape); 00265 00266 int number_of_regions = labelVolume(srcMultiArrayRange(src), destMultiArray(labels), 00267 neighborhood, equal); 00268 00269 // assume that a region is a extremum until the opposite is proved 00270 ArrayVector<unsigned char> isExtremum(number_of_regions+1, (unsigned char)1); 00271 00272 SrcIterator zs = src.traverser_begin(); 00273 LabelIterator zl = labels.traverser_begin(); 00274 for(MultiArrayIndex z = 0; z != d; ++z, ++zs.dim2(), ++zl.dim2()) 00275 { 00276 SrcIterator ys(zs); 00277 LabelIterator yl(zl); 00278 00279 for(MultiArrayIndex y = 0; y != h; ++y, ++ys.dim1(), ++yl.dim1()) 00280 { 00281 SrcIterator xs(ys); 00282 LabelIterator xl(yl); 00283 00284 for(MultiArrayIndex x = 0; x != w; ++x, ++xs.dim0(), ++xl.dim0()) 00285 { 00286 int lab = *xl; 00287 T1 v = *xs; 00288 00289 if(isExtremum[lab] == 0) 00290 continue; 00291 00292 if(!compare(v, threshold)) 00293 { 00294 // mark all regions that don't exceed the threshold as non-extremum 00295 isExtremum[lab] = 0; 00296 continue; 00297 } 00298 00299 AtVolumeBorder atBorder = isAtVolumeBorder(x, y, z, w, h, d); 00300 if(atBorder == NotAtBorder) 00301 { 00302 NeighborhoodCirculator<SrcIterator, Neighborhood> cs(xs); 00303 NeighborhoodCirculator<LabelIterator, Neighborhood> cl(xl); 00304 for(i=0; i<Neighborhood::DirectionCount; ++i, ++cs, ++cl) 00305 { 00306 if(lab != *cl && compare(*cs,v)) 00307 { 00308 isExtremum[lab] = 0; 00309 break; 00310 } 00311 } 00312 } 00313 else 00314 { 00315 if(allowExtremaAtBorder) 00316 { 00317 RestrictedNeighborhoodCirculator<SrcIterator, Neighborhood> 00318 cs(xs, atBorder), scend(cs); 00319 do 00320 { 00321 if(lab != *(xl+cs.diff()) && compare(cs,v)) 00322 { 00323 isExtremum[lab] = 0; 00324 break; 00325 } 00326 } 00327 while(++cs != scend); 00328 } 00329 else 00330 { 00331 isExtremum[lab] = 0; 00332 } 00333 } 00334 } 00335 } 00336 } 00337 00338 00339 zl = labels.traverser_begin(); 00340 DestIterator zd = dest.traverser_begin(); 00341 for(MultiArrayIndex z = 0; z != d; ++z, ++zl.dim2(), ++zd.dim2()) 00342 { 00343 LabelIterator yl(zl); 00344 DestIterator yd(zd); 00345 00346 for(MultiArrayIndex y = 0; y != h; ++y, ++yl.dim1(), ++yd.dim1()) 00347 { 00348 LabelIterator xl(yl); 00349 DestIterator xd(yd); 00350 00351 for(MultiArrayIndex x = 0; x != w; ++x, ++xl.dim0(), ++xd.dim0()) 00352 { 00353 if(isExtremum[*xl]) 00354 *xd = marker; 00355 } 00356 } 00357 } 00358 } 00359 00360 } // namespace detail 00361 00362 /********************************************************/ 00363 /* */ 00364 /* localMinima */ 00365 /* */ 00366 /********************************************************/ 00367 00368 // documentation is in localminmax.hxx 00369 template <unsigned int N, class T1, class C1, class T2, class C2> 00370 void 00371 localMinima(MultiArrayView<N, T1, C1> src, 00372 MultiArrayView<N, T2, C2> dest, 00373 LocalMinmaxOptions const & options = LocalMinmaxOptions()) 00374 { 00375 T1 threshold = options.use_threshold 00376 ? std::min(NumericTraits<T1>::max(), (T1)options.thresh) 00377 : NumericTraits<T1>::max(); 00378 T2 marker = (T2)options.marker; 00379 00380 vigra_precondition(!options.allow_plateaus, 00381 "localMinima(): Option 'allowPlateaus' is not implemented for arbitrary dimensions,\n" 00382 " use extendedLocalMinima() for 2D and 3D problems."); 00383 00384 if(options.neigh == 0) 00385 options.neigh = 2*N; 00386 if(options.neigh == 1) 00387 options.neigh = pow(3, N) - 1; 00388 00389 detail::localMinMax(src, dest, marker, options.neigh, 00390 threshold, std::less<T1>(), options.allow_at_border); 00391 } 00392 00393 /********************************************************/ 00394 /* */ 00395 /* localMaxima */ 00396 /* */ 00397 /********************************************************/ 00398 00399 // documentation is in localminmax.hxx 00400 template <unsigned int N, class T1, class C1, class T2, class C2> 00401 void 00402 localMaxima(MultiArrayView<N, T1, C1> src, 00403 MultiArrayView<N, T2, C2> dest, 00404 LocalMinmaxOptions const & options = LocalMinmaxOptions()) 00405 { 00406 T1 threshold = options.use_threshold 00407 ? std::max(NumericTraits<T1>::min(), (T1)options.thresh) 00408 : NumericTraits<T1>::min(); 00409 T2 marker = (T2)options.marker; 00410 00411 vigra_precondition(!options.allow_plateaus, 00412 "localMaxima(): Option 'allowPlateaus' is not implemented for arbitrary dimensions,\n" 00413 " use extendedLocalMinima() for 2D and 3D problems."); 00414 00415 if(options.neigh == 0) 00416 options.neigh = 2*N; 00417 if(options.neigh == 1) 00418 options.neigh = pow(3, N) - 1; 00419 00420 detail::localMinMax(src, dest, marker, options.neigh, 00421 threshold, std::greater<T1>(), options.allow_at_border); 00422 } 00423 00424 /**************************************************************************/ 00425 00426 /********************************************************/ 00427 /* */ 00428 /* extendedLocalMinima */ 00429 /* */ 00430 /********************************************************/ 00431 00432 // documentation is in localminmax.hxx 00433 template <class T1, class C1, class T2, class C2, 00434 class Neighborhood, class EqualityFunctor> 00435 inline void 00436 extendedLocalMinima(MultiArrayView<3, T1, C1> src, 00437 MultiArrayView<3, T2, C2> dest, 00438 LocalMinmaxOptions const & options = LocalMinmaxOptions()) 00439 { 00440 T1 threshold = options.use_threshold 00441 ? std::min(NumericTraits<T1>::max(), (T1)options.thresh) 00442 : NumericTraits<T1>::max(); 00443 T2 marker = (T2)options.marker; 00444 00445 if(options.neigh == 0 || options.neigh == 6) 00446 { 00447 detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DSix(), 00448 threshold, std::less<T1>(), std::equal_to<T1>(), 00449 options.allow_at_border); 00450 } 00451 else if(options.neigh == 1 || options.neigh == 26) 00452 { 00453 detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DTwentySix(), 00454 threshold, std::less<T1>(), std::equal_to<T1>(), 00455 options.allow_at_border); 00456 } 00457 else 00458 vigra_precondition(false, 00459 "extendedLocalMinima(): Invalid neighborhood."); 00460 } 00461 00462 /********************************************************/ 00463 /* */ 00464 /* extendedLocalMaxima */ 00465 /* */ 00466 /********************************************************/ 00467 00468 // documentation is in localminmax.hxx 00469 template <class T1, class C1, class T2, class C2, 00470 class Neighborhood, class EqualityFunctor> 00471 inline void 00472 extendedLocalMaxima(MultiArrayView<3, T1, C1> src, 00473 MultiArrayView<3, T2, C2> dest, 00474 LocalMinmaxOptions const & options = LocalMinmaxOptions()) 00475 { 00476 T1 threshold = options.use_threshold 00477 ? std::max(NumericTraits<T1>::min(), (T1)options.thresh) 00478 : NumericTraits<T1>::min(); 00479 T2 marker = (T2)options.marker; 00480 00481 if(options.neigh == 0 || options.neigh == 6) 00482 { 00483 detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DSix(), 00484 threshold, std::greater<T1>(), std::equal_to<T1>(), 00485 options.allow_at_border); 00486 } 00487 else if(options.neigh == 1 || options.neigh == 26) 00488 { 00489 detail::extendedLocalMinMax(src, dest, marker, NeighborCode3DTwentySix(), 00490 threshold, std::greater<T1>(), std::equal_to<T1>(), 00491 options.allow_at_border); 00492 } 00493 else 00494 vigra_precondition(false, 00495 "extendedLocalMaxima(): Invalid neighborhood."); 00496 } 00497 00498 } // namespace vigra 00499 00500 #endif 00501 00502 #endif // VIGRA_MULTI_LOCALMINMAX_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|