[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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) |
html generated using doxygen and Python
|