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