[ 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 explicitly 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> <vigra/gaborfilter.hxx><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 redundancy. 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> <vigra/gaborfilter.hxx> 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) |
html generated using doxygen and Python
|