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