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