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

details vigra/gaborfilter.hxx VIGRA

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