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

details vigra/gaborfilter.hxx VIGRA

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