[ 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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)