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

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