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

| vigra/gaborfilter.hxx |
--- | --- | ---

    00001 #ifndef VIGRA_GABORFILTER_HXX
    00002 #define VIGRA_GABORFILTER_HXX
    00003
    00004 #include "vigra/imagecontainer.hxx"
    00005 #include "vigra/config.hxx"
    00006 #include "vigra/stdimage.hxx"
    00007 #include "vigra/copyimage.hxx"
    00008 #include "vigra/transformimage.hxx"
    00009 #include "vigra/combineimages.hxx"
    00010 #include "vigra/utilities.hxx"
    00011
    00012 #include <functional>
    00013 #include <vector>
    00014 #include <cmath>
    00015
    00016 namespace vigra {
    00017
    00018 /** \addtogroup GaborFilter Gabor Filter
    00019     Functions to create or apply gabor filter (latter based on FFTW).
    00020 */
    00021 //@{
    00022
    00023 /********************************************************/
    00024 /*                                                      */
    00025 /*                   createGaborFilter                  */
    00026 /*                                                      */
    00027 /********************************************************/
    00028
    00029 /** \brief Create a gabor filter in frequency space.
    00030
    00031     The orientation is given in radians, the other parameters are the
    00032     center frequency (for example 0.375 or smaller) and the two
    00033     angular and radial sigmas of the gabor filter. (See \ref
    00034     angularGaborSigma() for an explanation of possible values.)
    00035
    00036     The energy of the filter is explicitely normalized to 1.0.
    00037
    00038     <b> Declarations:</b>
    00039
    00040     pass arguments explicitly:
    00041     \code
    00042     namespace vigra {
    00043         template <class DestImageIterator, class DestAccessor>
    00044         void createGaborFilter(DestImageIterator destUpperLeft,
    00045                                DestImageIterator destLowerRight,
    00046                                DestAccessor da,
    00047                                double orientation, double centerFrequency,
    00048                                double angularSigma, double radialSigma)
    00049     }
    00050     \endcode
    00051
    00052     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00053     \code
    00054     namespace vigra {
    00055         template <class DestImageIterator, class DestAccessor>
    00056         inline
    00057         void createGaborFilter(triple<DestImageIterator,
    00058                                       DestImageIterator,
    00059                                       DestAccessor> dest,
    00060                                double orientation, double centerFrequency,
    00061                                double angularSigma, double radialSigma)
    00062     }
    00063     \endcode
    00064
    00065     <b> Usage:</b>
    00066
    00067     <b>\#include</b> "<a href="gaborfilter_8hxx-source.html">vigra/gaborfilter.hxx</a>"<br>
    00068
    00069     Namespace: vigra
    00070
    00071     \code
    00072     vigra::FImage gabor(w,h);
    00073
    00074     vigra::createGaborFilter(destImageRange(gabor), orient, freq,
    00075                              angularGaborSigma(directionCount, freq)
    00076                              radialGaborSigma(freq));
    00077     \endcode
    00078 */
    00079
    00080 template <class DestImageIterator, class DestAccessor>
    00081 void createGaborFilter(DestImageIterator destUpperLeft,
    00082                        DestImageIterator destLowerRight, DestAccessor da,
    00083                        double orientation, double centerFrequency,
    00084                        double angularSigma, double radialSigma)
    00085 {
    00086     int w= destLowerRight.x - destUpperLeft.x;
    00087     int h= destLowerRight.y - destUpperLeft.y;
    00088
    00089     double squaredSum = 0.0;
    00090     double cosTheta= VIGRA_CSTD::cos(orientation);
    00091     double sinTheta= VIGRA_CSTD::sin(orientation);
    00092
    00093     double radialSigma2 = radialSigma*radialSigma;
    00094     double angularSigma2 = angularSigma*angularSigma;
    00095
    00096     double wscale = w % 1 ?
    00097                     1.0f / (w-1) :
    00098                     1.0f / w;
    00099     double hscale = h % 1 ?
    00100                     1.0f / (h-1) :
    00101                     1.0f / h;
    00102
    00103     int dcX= (w+1)/2, dcY= (h+1)/2;
    00104
    00105     double u, v;
    00106     for ( int y=0; y<h; y++, destUpperLeft.y++ )
    00107     {
    00108         typename DestImageIterator::row_iterator dix = destUpperLeft.rowIterator();
    00109
    00110         v = hscale * ((h - (y - dcY))%h - dcY);
    00111         for ( int x=0; x<w; x++, dix++ )
    00112         {
    00113             u= wscale*((x - dcX + w)%w - dcX);
    00114
    00115             double uu = cosTheta*u + sinTheta*v - centerFrequency;
    00116             double vv = -sinTheta*u + cosTheta*v;
    00117             double gabor;
    00118
    00119             gabor = VIGRA_CSTD::exp(-0.5*(uu*uu / radialSigma2 + vv*vv / angularSigma2));
    00120             squaredSum += gabor * gabor;
    00121             da.set( gabor, dix );
    00122         }
    00123     }
    00124     destUpperLeft.y -= h;
    00125
    00126     // clear out DC value and remove it from the squared sum
    00127     double dcValue = da(destUpperLeft);
    00128     squaredSum -= dcValue * dcValue;
    00129     da.set( 0.0, destUpperLeft );
    00130
    00131     // normalize energy to one
    00132     double factor = VIGRA_CSTD::sqrt(squaredSum);
    00133     for ( int y=0; y<h; y++, destUpperLeft.y++ )
    00134     {
    00135         typename DestImageIterator::row_iterator dix = destUpperLeft.rowIterator();
    00136
    00137         for ( int x=0; x<w; x++, dix++ )
    00138         {
    00139             da.set( da(dix) / factor, dix );
    00140         }
    00141     }
    00142 }
    00143
    00144 template <class DestImageIterator, class DestAccessor>
    00145 inline
    00146 void createGaborFilter(triple<DestImageIterator, DestImageIterator,
    00147                               DestAccessor> dest,
    00148                               double orientation, double centerFrequency,
    00149                               double angularSigma, double radialSigma)
    00150 {
    00151     createGaborFilter(dest.first, dest.second, dest.third,
    00152                       orientation, centerFrequency,
    00153                       angularSigma, radialSigma);
    00154 }
    00155
    00156 /********************************************************/
    00157 /*                                                      */
    00158 /*                   radialGaborSigma                   */
    00159 /*                                                      */
    00160 /********************************************************/
    00161
    00162 /** \brief Calculate sensible radial sigma for given parameters.
    00163
    00164     For a brief introduction what is meant with "sensible" sigmas, see
    00165     \ref angularGaborSigma().
    00166
    00167     <b> Declaration:</b>
    00168
    00169     \code
    00170     namespace vigra {
    00171         inline
    00172         double radialGaborSigma(double centerFrequency)
    00173     }
    00174     \endcode
    00175  */
    00176
    00177 inline double radialGaborSigma(double centerFrequency)
    00178 {
    00179     static double sfactor = 3.0 * VIGRA_CSTD::sqrt(VIGRA_CSTD::log(4.0));
    00180     return centerFrequency / sfactor;
    00181 }
    00182
    00183 /********************************************************/
    00184 /*                                                      */
    00185 /*                   angularGaborSigma                  */
    00186 /*                                                      */
    00187 /********************************************************/
    00188
    00189 /** \brief Calculate sensible angular sigma for given parameters.
    00190
    00191     "Sensible" means: If you use a range of gabor filters for feature
    00192     detection, you are interested in minimal redundance. This is hard
    00193     to define but one possible try is to arrange the filters in
    00194     frequency space, so that the half-peak-magnitude ellipses touch
    00195     each other.
    00196
    00197     To do so, you must know the number of directions (first parameter
    00198     for the angular sigma function) and the center frequency of the
    00199     filter you want to calculate the sigmas for.
    00200
    00201     The exact formulas are:
    00202     \code
    00203     sigma_radial= 1/sqrt(ln(4)) * centerFrequency/3
    00204     \endcode
    00205
    00206     \code
    00207     sigma_angular= 1/sqrt(ln(4)) * tan(pi/(directions*2))
    00208                    * sqrt(8/9) * centerFrequency
    00209     \endcode
    00210
    00211     <b> Declaration:</b>
    00212
    00213     \code
    00214     namespace vigra {
    00215         inline
    00216         double angularGaborSigma(int directionCount, double centerFrequency)
    00217     }
    00218     \endcode
    00219  */
    00220
    00221 inline double angularGaborSigma(int directionCount, double centerFrequency)
    00222 {
    00223     return VIGRA_CSTD::tan(M_PI/directionCount/2.0) * centerFrequency
    00224         * VIGRA_CSTD::sqrt(8.0 / (9 * VIGRA_CSTD::log(4.0)));
    00225 }
    00226
    00227 /********************************************************/
    00228 /*                                                      */
    00229 /*                   GaborFilterFamily                  */
    00230 /*                                                      */
    00231 /********************************************************/
    00232
    00233 /** \brief Family of gabor filters of different scale and direction.
    00234
    00235     A GaborFilterFamily can be used to quickly create a whole family
    00236     of gabor filters in frequency space. Especially useful in
    00237     conjunction with \ref applyFourierFilterFamily, since it's derived
    00238     from \ref ImageArray.
    00239
    00240     The filter parameters are chosen to make the center frequencies
    00241     decrease in octaves with increasing scale indices, and to make the
    00242     half-peak-magnitude ellipses touch each other to somewhat reduce
    00243     redundancy in the filter answers. This is done by using \ref
    00244     angularGaborSigma() and \ref radialGaborSigma(), you'll find more
    00245     information there.
    00246
    00247     The template parameter ImageType should be a scalar image type suitable for filling in
    00248
    00249     <b>\#include</b> "<a href="gaborfilter_8hxx-source.html">vigra/gaborfilter.hxx</a>"
    00250
    00251     Namespace: vigra
    00252 */
    00253 template <class ImageType,
    00254       class Alloc = typename ImageType::allocator_type::template rebind<ImageType>::other >
    00255 class GaborFilterFamily
    00256 : public ImageArray<ImageType, Alloc>
    00257 {
    00258     typedef ImageArray<ImageType, Alloc> ParentClass;
    00259     int scaleCount_, directionCount_;
    00260     double maxCenterFrequency_;
    00261
    00262 protected:
    00263     void initFilters()
    00264     {
    00265         for(int direction= 0; direction<directionCount_; direction++)
    00266             for(int scale= 0; scale<scaleCount_; scale++)
    00267             {
    00268                 double angle = direction * M_PI / directionCount();
    00269                 double centerFrequency =
    00270                     maxCenterFrequency_ / VIGRA_CSTD::pow(2.0, (double)scale);
    00271                 createGaborFilter(destImageRange(this->images_[filterIndex(direction, scale)]),
    00272                                   angle, centerFrequency,
    00273                                   angularGaborSigma(directionCount(), centerFrequency),
    00274                                   radialGaborSigma(centerFrequency));
    00275             }
    00276     }
    00277
    00278 public:
    00279     enum { stdFilterSize= 128, stdDirectionCount= 6, stdScaleCount= 4 };
    00280
    00281         /** Constructs a family of gabor filters in frequency
    00282             space. The filters will be calculated on construction, so
    00283             it makes sense to provide good parameters right now
    00284             although they can be changed later, too. If you leave them
    00285             out, the defaults are a \ref directionCount of 6, a \ref
    00286             scaleCount of 4 and a \ref maxCenterFrequency of
    00287             3/8(=0.375).
    00288         */
    00289     GaborFilterFamily(const Diff2D & size,
    00290                       int directionCount = stdDirectionCount, int scaleCount = stdScaleCount,
    00291                       double maxCenterFrequency = 3.0/8.0,
    00292                       Alloc const & alloc = Alloc())
    00293         : ParentClass(directionCount*scaleCount, size, alloc),
    00294           scaleCount_(scaleCount),
    00295           directionCount_(directionCount),
    00296           maxCenterFrequency_(maxCenterFrequency)
    00297     {
    00298         initFilters();
    00299     }
    00300
    00301         /** Convenience variant of the above constructor taking width
    00302             and height separately. Also, this one serves as default
    00303             constructor constructing 128x128 pixel filters.
    00304          */
    00305     GaborFilterFamily(int width= stdFilterSize, int height= -1,
    00306                       int directionCount = stdDirectionCount, int scaleCount = stdScaleCount,
    00307                       double maxCenterFrequency = 3.0/8.0,
    00308                       Alloc const & alloc = Alloc())
    00309         : ParentClass(directionCount*scaleCount,
    00310                       Size2D(width, height > 0 ? height : width), alloc),
    00311           scaleCount_(scaleCount),
    00312           directionCount_(directionCount),
    00313           maxCenterFrequency_(maxCenterFrequency)
    00314     {
    00315         initFilters();
    00316     }
    00317
    00318         /** Return the index of the filter with the given direction and
    00319             scale in this ImageArray. direction must in the range
    00320             0..directionCount()-1 and scale in the range
    00321             0..rangeCount()-1. This is useful for example if you used
    00322             \ref applyFourierFilterFamily() and got a resulting
    00323             ImageArray which still has the same order of images, but no
    00324             \ref getFilter() method anymore.
    00325          */
    00326     int filterIndex(int direction, int scale) const
    00327     {
    00328         return scale*directionCount()+direction;
    00329     }
    00330
    00331         /** Return the filter with the given direction and
    00332             scale. direction must in the range 0..directionCount()-1
    00333             and scale in the range 0..rangeCount()-1.
    00334             <tt>filters.getFilter(direction, scale)</tt> is the same as
    00335             <tt>filters[filterIndex(direction, scale)]</tt>.
    00336          */
    00337     ImageType const & getFilter(int direction, int scale) const
    00338     {
    00339         return this->images_[filterIndex(direction, scale)];
    00340     }
    00341
    00342         /** Resize all filters (causing their recalculation).
    00343          */
    00344     virtual void resizeImages(const Diff2D &newSize)
    00345     {
    00346         ParentClass::resizeImages(newSize);
    00347         initFilters();
    00348     }
    00349
    00350         /** Query the number of filter scales available.
    00351          */
    00352     int scaleCount() const
    00353         { return scaleCount_; }
    00354
    00355         /** Query the number of filter directions available.
    00356          */
    00357     int directionCount() const
    00358         { return directionCount_; }
    00359
    00360         /** Change the number of directions / scales. This causes the
    00361             recalculation of all filters.
    00362          */
    00363     void setDirectionScaleCounts(int directionCount, int scaleCount)
    00364     {
    00365         this->resize(directionCount * scaleCount);
    00366         scaleCount_ = scaleCount;
    00367         directionCount_ = directionCount;
    00368         initFilters();
    00369     }
    00370
    00371         /** Return the center frequency of the filter(s) with
    00372             scale==0. Filters with scale>0 will have a center frequency
    00373             reduced in octaves:
    00374             <tt>centerFrequency= maxCenterFrequency / 2.0^scale</tt>
    00375         */
    00376     double maxCenterFrequency()
    00377         { return maxCenterFrequency_; }
    00378
    00379         /** Change the center frequency of the filter(s) with
    00380             scale==0. See \ref maxCenterFrequency().
    00381          */
    00382     void setMaxCenterFrequency(double maxCenterFrequency)
    00383     {
    00384         maxCenterFrequency_ = maxCenterFrequency;
    00385         initFilters();
    00386     }
    00387 };
    00388
    00389 //@}
    00390
    00391 } // namespace vigra
    00392
    00393 #endif // VIGRA_GABORFILTER_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)