[ 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 
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)