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

vigra/multi_localminmax.hxx VIGRA

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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (20 Sep 2011)