[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

gaborfilter.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2002-2004 by Ullrich Koethe and Hans Meine */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_GABORFILTER_HXX
38 #define VIGRA_GABORFILTER_HXX
39 
40 #include "imagecontainer.hxx"
41 #include "config.hxx"
42 #include "stdimage.hxx"
43 #include "copyimage.hxx"
44 #include "transformimage.hxx"
45 #include "combineimages.hxx"
46 #include "utilities.hxx"
47 
48 #include <functional>
49 #include <vector>
50 #include <cmath>
51 
52 namespace vigra {
53 
54 /** \addtogroup GaborFilter Gabor Filter
55  Functions to create or apply gabor filter (latter based on FFTW).
56 */
57 //@{
58 
59 /********************************************************/
60 /* */
61 /* createGaborFilter */
62 /* */
63 /********************************************************/
64 
65 /** \brief Create a gabor filter in frequency space.
66 
67  The orientation is given in radians, the other parameters are the
68  center frequency (for example 0.375 or smaller) and the two
69  angular and radial sigmas of the gabor filter. (See \ref
70  angularGaborSigma() for an explanation of possible values.)
71 
72  The energy of the filter is explicitely normalized to 1.0.
73 
74  <b> Declarations:</b>
75 
76  pass arguments explicitly:
77  \code
78  namespace vigra {
79  template <class DestImageIterator, class DestAccessor>
80  void createGaborFilter(DestImageIterator destUpperLeft,
81  DestImageIterator destLowerRight,
82  DestAccessor da,
83  double orientation, double centerFrequency,
84  double angularSigma, double radialSigma)
85  }
86  \endcode
87 
88  use argument objects in conjunction with \ref ArgumentObjectFactories :
89  \code
90  namespace vigra {
91  template <class DestImageIterator, class DestAccessor>
92  void createGaborFilter(triple<DestImageIterator,
93  DestImageIterator,
94  DestAccessor> dest,
95  double orientation, double centerFrequency,
96  double angularSigma, double radialSigma)
97  }
98  \endcode
99 
100  <b> Usage:</b>
101 
102  <b>\#include</b> <<a href="gaborfilter_8hxx-source.html">vigra/gaborfilter.hxx</a>><br>
103 
104  Namespace: vigra
105 
106  \code
107  vigra::FImage gabor(w,h);
108 
109  vigra::createGaborFilter(destImageRange(gabor), orient, freq,
110  angularGaborSigma(directionCount, freq)
111  radialGaborSigma(freq));
112  \endcode
113 */
115 
116 template <class DestImageIterator, class DestAccessor>
117 void createGaborFilter(DestImageIterator destUpperLeft,
118  DestImageIterator destLowerRight, DestAccessor da,
119  double orientation, double centerFrequency,
120  double angularSigma, double radialSigma)
121 {
122  int w = int(destLowerRight.x - destUpperLeft.x);
123  int h = int(destLowerRight.y - destUpperLeft.y);
124 
125  double squaredSum = 0.0;
126  double cosTheta= VIGRA_CSTD::cos(orientation);
127  double sinTheta= VIGRA_CSTD::sin(orientation);
128 
129  double radialSigma2 = radialSigma*radialSigma;
130  double angularSigma2 = angularSigma*angularSigma;
131 
132  double wscale = w % 1 ?
133  1.0f / (w-1) :
134  1.0f / w;
135  double hscale = h % 1 ?
136  1.0f / (h-1) :
137  1.0f / h;
138 
139  int dcX= (w+1)/2, dcY= (h+1)/2;
140 
141  double u, v;
142  for ( int y=0; y<h; y++, destUpperLeft.y++ )
143  {
144  typename DestImageIterator::row_iterator dix = destUpperLeft.rowIterator();
145 
146  v = hscale * ((h - (y - dcY))%h - dcY);
147  for ( int x=0; x<w; x++, dix++ )
148  {
149  u= wscale*((x - dcX + w)%w - dcX);
150 
151  double uu = cosTheta*u + sinTheta*v - centerFrequency;
152  double vv = -sinTheta*u + cosTheta*v;
153  double gabor;
154 
155  gabor = VIGRA_CSTD::exp(-0.5*(uu*uu / radialSigma2 + vv*vv / angularSigma2));
156  squaredSum += gabor * gabor;
157  da.set( gabor, dix );
158  }
159  }
160  destUpperLeft.y -= h;
161 
162  // clear out DC value and remove it from the squared sum
163  double dcValue = da(destUpperLeft);
164  squaredSum -= dcValue * dcValue;
165  da.set( 0.0, destUpperLeft );
166 
167  // normalize energy to one
168  double factor = VIGRA_CSTD::sqrt(squaredSum);
169  for ( int y=0; y<h; y++, destUpperLeft.y++ )
170  {
171  typename DestImageIterator::row_iterator dix = destUpperLeft.rowIterator();
172 
173  for ( int x=0; x<w; x++, dix++ )
174  {
175  da.set( da(dix) / factor, dix );
176  }
177  }
178 }
179 
180 template <class DestImageIterator, class DestAccessor>
181 inline
182 void createGaborFilter(triple<DestImageIterator, DestImageIterator,
183  DestAccessor> dest,
184  double orientation, double centerFrequency,
185  double angularSigma, double radialSigma)
186 {
187  createGaborFilter(dest.first, dest.second, dest.third,
188  orientation, centerFrequency,
189  angularSigma, radialSigma);
190 }
191 
192 /********************************************************/
193 /* */
194 /* radialGaborSigma */
195 /* */
196 /********************************************************/
197 
198 /** \brief Calculate sensible radial sigma for given parameters.
199 
200  For a brief introduction what is meant with "sensible" sigmas, see
201  \ref angularGaborSigma().
202 
203  <b> Declaration:</b>
204 
205  \code
206  namespace vigra {
207  double radialGaborSigma(double centerFrequency)
208  }
209  \endcode
210  */
211 
212 inline double radialGaborSigma(double centerFrequency)
213 {
214  static double sfactor = 3.0 * VIGRA_CSTD::sqrt(VIGRA_CSTD::log(4.0));
215  return centerFrequency / sfactor;
216 }
217 
218 /********************************************************/
219 /* */
220 /* angularGaborSigma */
221 /* */
222 /********************************************************/
223 
224 /** \brief Calculate sensible angular sigma for given parameters.
225 
226  "Sensible" means: If you use a range of gabor filters for feature
227  detection, you are interested in minimal redundance. This is hard
228  to define but one possible try is to arrange the filters in
229  frequency space, so that the half-peak-magnitude ellipses touch
230  each other.
231 
232  To do so, you must know the number of directions (first parameter
233  for the angular sigma function) and the center frequency of the
234  filter you want to calculate the sigmas for.
235 
236  The exact formulas are:
237  \code
238  sigma_radial= 1/sqrt(ln(4)) * centerFrequency/3
239  \endcode
240 
241  \code
242  sigma_angular= 1/sqrt(ln(4)) * tan(pi/(directions*2))
243  * sqrt(8/9) * centerFrequency
244  \endcode
245 
246  <b> Declaration:</b>
247 
248  \code
249  namespace vigra {
250  double angularGaborSigma(int directionCount, double centerFrequency)
251  }
252  \endcode
253  */
254 
255 inline double angularGaborSigma(int directionCount, double centerFrequency)
256 {
257  return VIGRA_CSTD::tan(M_PI/directionCount/2.0) * centerFrequency
258  * VIGRA_CSTD::sqrt(8.0 / (9 * VIGRA_CSTD::log(4.0)));
259 }
260 
261 /********************************************************/
262 /* */
263 /* GaborFilterFamily */
264 /* */
265 /********************************************************/
266 
267 /** \brief Family of gabor filters of different scale and direction.
268 
269  A GaborFilterFamily can be used to quickly create a whole family
270  of gabor filters in frequency space. Especially useful in
271  conjunction with \ref applyFourierFilterFamily, since it's derived
272  from \ref ImageArray.
273 
274  The filter parameters are chosen to make the center frequencies
275  decrease in octaves with increasing scale indices, and to make the
276  half-peak-magnitude ellipses touch each other to somewhat reduce
277  redundancy in the filter answers. This is done by using \ref
278  angularGaborSigma() and \ref radialGaborSigma(), you'll find more
279  information there.
280 
281  The template parameter ImageType should be a scalar image type suitable for filling in
282 
283  <b>\#include</b> <<a href="gaborfilter_8hxx-source.html">vigra/gaborfilter.hxx</a>>
284 
285  Namespace: vigra
286 */
287 template <class ImageType,
288  class Alloc = typename ImageType::allocator_type::template rebind<ImageType>::other >
290 : public ImageArray<ImageType, Alloc>
291 {
293  int scaleCount_, directionCount_;
294  double maxCenterFrequency_;
295 
296 protected:
297  void initFilters()
298  {
299  for(int direction= 0; direction<directionCount_; direction++)
300  for(int scale= 0; scale<scaleCount_; scale++)
301  {
302  double angle = direction * M_PI / directionCount();
303  double centerFrequency =
304  maxCenterFrequency_ / VIGRA_CSTD::pow(2.0, (double)scale);
305  createGaborFilter(destImageRange(this->images_[filterIndex(direction, scale)]),
306  angle, centerFrequency,
307  angularGaborSigma(directionCount(), centerFrequency),
308  radialGaborSigma(centerFrequency));
309  }
310  }
311 
312 public:
313  enum { stdFilterSize= 128, stdDirectionCount= 6, stdScaleCount= 4 };
314 
315  /** Constructs a family of gabor filters in frequency
316  space. The filters will be calculated on construction, so
317  it makes sense to provide good parameters right now
318  although they can be changed later, too. If you leave them
319  out, the defaults are a \ref directionCount of 6, a \ref
320  scaleCount of 4 and a \ref maxCenterFrequency of
321  3/8(=0.375).
322  */
324  int directionCount = stdDirectionCount, int scaleCount = stdScaleCount,
325  double maxCenterFrequency = 3.0/8.0,
326  Alloc const & alloc = Alloc())
327  : ParentClass(directionCount*scaleCount, size, alloc),
328  scaleCount_(scaleCount),
329  directionCount_(directionCount),
330  maxCenterFrequency_(maxCenterFrequency)
331  {
332  initFilters();
333  }
334 
335  /** Convenience variant of the above constructor taking width
336  and height separately. Also, this one serves as default
337  constructor constructing 128x128 pixel filters.
338  */
339  GaborFilterFamily(int width= stdFilterSize, int height= -1,
340  int directionCount = stdDirectionCount, int scaleCount = stdScaleCount,
341  double maxCenterFrequency = 3.0/8.0,
342  Alloc const & alloc = Alloc())
344  Size2D(width, height > 0 ? height : width), alloc),
345  scaleCount_(scaleCount),
346  directionCount_(directionCount),
347  maxCenterFrequency_(maxCenterFrequency)
348  {
349  initFilters();
350  }
351 
352  /** Return the index of the filter with the given direction and
353  scale in this ImageArray. direction must in the range
354  0..directionCount()-1 and scale in the range
355  0..rangeCount()-1. This is useful for example if you used
356  \ref applyFourierFilterFamily() and got a resulting
357  ImageArray which still has the same order of images, but no
358  \ref getFilter() method anymore.
359  */
360  int filterIndex(int direction, int scale) const
361  {
362  return scale*directionCount()+direction;
363  }
364 
365  /** Return the filter with the given direction and
366  scale. direction must in the range 0..directionCount()-1
367  and scale in the range 0..rangeCount()-1.
368  <tt>filters.getFilter(direction, scale)</tt> is the same as
369  <tt>filters[filterIndex(direction, scale)]</tt>.
370  */
371  ImageType const & getFilter(int direction, int scale) const
372  {
373  return this->images_[filterIndex(direction, scale)];
374  }
375 
376  /** Resize all filters (causing their recalculation).
377  */
378  virtual void resizeImages(const Diff2D &newSize)
379  {
380  ParentClass::resizeImages(newSize);
381  initFilters();
382  }
383 
384  /** Query the number of filter scales available.
385  */
386  int scaleCount() const
387  { return scaleCount_; }
388 
389  /** Query the number of filter directions available.
390  */
391  int directionCount() const
392  { return directionCount_; }
393 
394  /** Change the number of directions / scales. This causes the
395  recalculation of all filters.
396  */
398  {
399  this->resize(directionCount * scaleCount);
400  scaleCount_ = scaleCount;
401  directionCount_ = directionCount;
402  initFilters();
403  }
404 
405  /** Return the center frequency of the filter(s) with
406  scale==0. Filters with scale>0 will have a center frequency
407  reduced in octaves:
408  <tt>centerFrequency= maxCenterFrequency / 2.0^scale</tt>
409  */
411  { return maxCenterFrequency_; }
412 
413  /** Change the center frequency of the filter(s) with
414  scale==0. See \ref maxCenterFrequency().
415  */
417  {
418  maxCenterFrequency_ = maxCenterFrequency;
419  initFilters();
420  }
421 };
422 
423 //@}
424 
425 } // namespace vigra
426 
427 #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.7.1 (Tue Jul 10 2012)