[ 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#include "multi_shape.hxx"
48
49#include <cmath>
50#include <functional>
51#include <memory>
52#include <vector>
53
54 namespace vigra {
55
56/** \addtogroup GaborFilter Gabor Filter
57 Functions to create or apply gabor filter (latter based on FFTW).
58*/
59//@{
60
61/********************************************************/
62/* */
63/* createGaborFilter */
64/* */
65/********************************************************/
66
67/** \brief Create a gabor filter in frequency space.
68
69 The orientation is given in radians, the other parameters are the
70 center frequency (for example 0.375 or smaller) and the two
71 angular and radial sigmas of the gabor filter. (See \ref
72 angularGaborSigma() for an explanation of possible values.)
73
74 The energy of the filter is explicitly normalized to 1.0.
75
76 <b> Declarations:</b>
77
78 pass 2D array views:
79 \code
80 namespace vigra {
81 template <class T, class S>
82 void
83 createGaborFilter(MultiArrayView<2, T, S> dest,
84 double orientation, double centerFrequency,
85 double angularSigma, double radialSigma);
86 }
87 \endcode
88
89 \deprecatedAPI{createGaborFilter}
90 pass \ref ImageIterators and \ref DataAccessors :
91 \code
92 namespace vigra {
93 template <class DestImageIterator, class DestAccessor>
94 void createGaborFilter(DestImageIterator destUpperLeft,
95 DestImageIterator destLowerRight,
96 DestAccessor da,
97 double orientation, double centerFrequency,
98 double angularSigma, double radialSigma)
99 }
100 \endcode
101 use argument objects in conjunction with \ref ArgumentObjectFactories :
102 \code
103 namespace vigra {
104 template <class DestImageIterator, class DestAccessor>
105 void createGaborFilter(triple<DestImageIterator,
106 DestImageIterator,
107 DestAccessor> dest,
108 double orientation, double centerFrequency,
109 double angularSigma, double radialSigma)
110 }
111 \endcode
112 \deprecatedEnd
113
114 <b> Usage:</b>
115
116 <b>\#include</b> <vigra/gaborfilter.hxx><br/>
117 Namespace: vigra
118
119 \code
120 MultiArray<2, float> gabor(w,h);
121
122 createGaborFilter(gabor, orient, freq,
123 angularGaborSigma(directionCount, freq),
124 radialGaborSigma(freq));
125 \endcode
126
127 \deprecatedUsage{createGaborFilter}
128 \code
129 vigra::FImage gabor(w,h);
130
131 vigra::createGaborFilter(destImageRange(gabor), orient, freq,
132 angularGaborSigma(directionCount, freq),
133 radialGaborSigma(freq));
134 \endcode
135 \deprecatedEnd
136*/
137doxygen_overloaded_function(template <...> void createGaborFilter)
138
139 template <class DestImageIterator, class DestAccessor>
140 void createGaborFilter(DestImageIterator destUpperLeft,
141 DestImageIterator destLowerRight, DestAccessor da,
142 double orientation, double centerFrequency,
143 double angularSigma, double radialSigma)
144{
145 int w = int(destLowerRight.x - destUpperLeft.x);
146 int h = int(destLowerRight.y - destUpperLeft.y);
147
148 double squaredSum = 0.0;
149 double cosTheta= VIGRA_CSTD::cos(orientation);
150 double sinTheta= VIGRA_CSTD::sin(orientation);
151
152 double radialSigma2 = radialSigma*radialSigma;
153 double angularSigma2 = angularSigma*angularSigma;
154
155 double wscale = w % 1 ?
156 1.0f / (w-1) :
157 1.0f / w;
158 double hscale = h % 1 ?
159 1.0f / (h-1) :
160 1.0f / h;
161
162 int dcX= (w+1)/2, dcY= (h+1)/2;
163
164 double u, v;
165 for ( int y=0; y<h; y++, destUpperLeft.y++ )
166 {
167 typename DestImageIterator::row_iterator dix = destUpperLeft.rowIterator();
168
169 v = hscale * ((h - (y - dcY))%h - dcY);
170 for ( int x=0; x<w; x++, dix++ )
171 {
172 u= wscale*((x - dcX + w)%w - dcX);
173
174 double uu = cosTheta*u + sinTheta*v - centerFrequency;
175 double vv = -sinTheta*u + cosTheta*v;
176 double gabor;
177
178 gabor = VIGRA_CSTD::exp(-0.5*(uu*uu / radialSigma2 + vv*vv / angularSigma2));
179 squaredSum += gabor * gabor;
180 da.set( gabor, dix );
181 }
182 }
183 destUpperLeft.y -= h;
184
185 // clear out DC value and remove it from the squared sum
186 double dcValue = da(destUpperLeft);
187 squaredSum -= dcValue * dcValue;
188 da.set( 0.0, destUpperLeft );
189
190 // normalize energy to one
191 double factor = VIGRA_CSTD::sqrt(squaredSum);
192 for ( int y=0; y<h; y++, destUpperLeft.y++ )
193 {
194 typename DestImageIterator::row_iterator dix = destUpperLeft.rowIterator();
195
196 for ( int x=0; x<w; x++, dix++ )
197 {
198 da.set( da(dix) / factor, dix );
199 }
200 }
201}
202
203 template <class DestImageIterator, class DestAccessor>
204 inline void
205 createGaborFilter(triple<DestImageIterator, DestImageIterator, DestAccessor> dest,
206 double orientation, double centerFrequency,
207 double angularSigma, double radialSigma)
208{
209 createGaborFilter(dest.first, dest.second, dest.third,
210 orientation, centerFrequency,
211 angularSigma, radialSigma);
212}
213
214 template <class T, class S>
215 inline void
216 createGaborFilter(MultiArrayView<2, T, S> dest,
217 double orientation, double centerFrequency,
218 double angularSigma, double radialSigma)
219{
220 createGaborFilter(destImageRange(dest),
221 orientation, centerFrequency,
222 angularSigma, radialSigma);
223}
224
225/********************************************************/
226/* */
227/* radialGaborSigma */
228/* */
229/********************************************************/
230
231/** \brief Calculate sensible radial sigma for given parameters.
232
233 For a brief introduction what is meant with "sensible" sigmas, see
234 \ref angularGaborSigma().
235
236 <b> Declaration:</b>
237
238 \code
239 namespace vigra {
240 double radialGaborSigma(double centerFrequency)
241 }
242 \endcode
243 */
244
245 inline double radialGaborSigma(double centerFrequency)
246{
247 double sfactor = 3.0 * VIGRA_CSTD::sqrt(VIGRA_CSTD::log(4.0));
248 return centerFrequency / sfactor;
249}
250
251/********************************************************/
252/* */
253/* angularGaborSigma */
254/* */
255/********************************************************/
256
257/** \brief Calculate sensible angular sigma for given parameters.
258
259 "Sensible" means: If you use a range of gabor filters for feature
260 detection, you are interested in minimal redundancy. This is hard
261 to define but one possible try is to arrange the filters in
262 frequency space, so that the half-peak-magnitude ellipses touch
263 each other.
264
265 To do so, you must know the number of directions (first parameter
266 for the angular sigma function) and the center frequency of the
267 filter you want to calculate the sigmas for.
268
269 The exact formulas are:
270 \code
271 sigma_radial= 1/sqrt(ln(4)) * centerFrequency/3
272 \endcode
273
274 \code
275 sigma_angular= 1/sqrt(ln(4)) * tan(pi/(directions*2))
276 * sqrt(8/9) * centerFrequency
277 \endcode
278
279 <b> Declaration:</b>
280
281 \code
282 namespace vigra {
283 double angularGaborSigma(int directionCount, double centerFrequency)
284 }
285 \endcode
286 */
287
288 inline double angularGaborSigma(int directionCount, double centerFrequency)
289{
290 return VIGRA_CSTD::tan(M_PI/directionCount/2.0) * centerFrequency
291 * VIGRA_CSTD::sqrt(8.0 / (9 * VIGRA_CSTD::log(4.0)));
292}
293
294/********************************************************/
295/* */
296/* GaborFilterFamily */
297/* */
298/********************************************************/
299
300/** \brief Family of gabor filters of different scale and direction.
301
302 A GaborFilterFamily can be used to quickly create a whole family
303 of gabor filters in frequency space. Especially useful in
304 conjunction with \ref applyFourierFilterFamily, since it's derived
305 from \ref ImageArray.
306
307 The filter parameters are chosen to make the center frequencies
308 decrease in octaves with increasing scale indices, and to make the
309 half-peak-magnitude ellipses touch each other to somewhat reduce
310 redundancy in the filter answers. This is done by using \ref
311 angularGaborSigma() and \ref radialGaborSigma(), you'll find more
312 information there.
313
314 The template parameter ImageType should be a scalar image type suitable for filling in
315
316 <b>\#include</b> <vigra/gaborfilter.hxx><br/>
317 Namespace: vigra
318*/
319 template <class ImageType,
320 class Alloc = typename std::allocator_traits<typename ImageType::allocator_type>::template rebind_alloc<ImageType> >
321 class GaborFilterFamily
322: public ImageArray<ImageType, Alloc>
323{
324 typedef ImageArray<ImageType, Alloc> ParentClass;
325 int scaleCount_, directionCount_;
326 double maxCenterFrequency_;
327
328 protected:
329 void initFilters()
330 {
331 for(int direction= 0; direction<directionCount_; direction++)
332 for(int scale= 0; scale<scaleCount_; scale++)
333 {
334 double angle = direction * M_PI / directionCount();
335 double centerFrequency =
336 maxCenterFrequency_ / VIGRA_CSTD::pow(2.0, (double)scale);
337 createGaborFilter(destImageRange(this->images_[filterIndex(direction, scale)]),
338 angle, centerFrequency,
339 angularGaborSigma(directionCount(), centerFrequency),
340 radialGaborSigma(centerFrequency));
341 }
342 }
343
344 public:
345 enum { stdFilterSize= 128, stdDirectionCount= 6, stdScaleCount= 4 };
346
347 /** Constructs a family of gabor filters in frequency
348 space. The filters will be calculated on construction, so
349 it makes sense to provide good parameters right now
350 although they can be changed later, too. If you leave them
351 out, the defaults are a \ref directionCount of 6, a \ref
352 scaleCount of 4 and a \ref maxCenterFrequency of
353 3/8(=0.375).
354 */
355 GaborFilterFamily(const Diff2D & size,
356 int directionCount = stdDirectionCount, int scaleCount = stdScaleCount,
357 double maxCenterFrequency = 3.0/8.0,
358 Alloc const & alloc = Alloc())
359 : ParentClass(directionCount*scaleCount, size, alloc),
360 scaleCount_(scaleCount),
361 directionCount_(directionCount),
362 maxCenterFrequency_(maxCenterFrequency)
363 {
364 initFilters();
365 }
366
367 /** Convenience variant of the above constructor taking width
368 and height separately. Also, this one serves as default
369 constructor constructing 128x128 pixel filters.
370 */
371 GaborFilterFamily(int width= stdFilterSize, int height= -1,
372 int directionCount = stdDirectionCount, int scaleCount = stdScaleCount,
373 double maxCenterFrequency = 3.0/8.0,
374 Alloc const & alloc = Alloc())
375 : ParentClass(directionCount*scaleCount,
376 Size2D(width, height > 0 ? height : width), alloc),
377 scaleCount_(scaleCount),
378 directionCount_(directionCount),
379 maxCenterFrequency_(maxCenterFrequency)
380 {
381 initFilters();
382 }
383
384 /** Return the index of the filter with the given direction and
385 scale in this ImageArray. direction must in the range
386 0..directionCount()-1 and scale in the range
387 0..rangeCount()-1. This is useful for example if you used
388 \ref applyFourierFilterFamily() and got a resulting
389 ImageArray which still has the same order of images, but no
390 \ref getFilter() method anymore.
391 */
392 int filterIndex(int direction, int scale) const
393 {
394 return scale*directionCount()+direction;
395 }
396
397 /** Return the filter with the given direction and
398 scale. direction must in the range 0..directionCount()-1
399 and scale in the range 0..rangeCount()-1.
400 <tt>filters.getFilter(direction, scale)</tt> is the same as
401 <tt>filters[filterIndex(direction, scale)]</tt>.
402 */
403 ImageType const & getFilter(int direction, int scale) const
404 {
405 return this->images_[filterIndex(direction, scale)];
406 }
407
408 /** Resize all filters (causing their recalculation).
409 */
410 virtual void resizeImages(const Diff2D &newSize)
411 {
412 ParentClass::resizeImages(newSize);
413 initFilters();
414 }
415
416 /** Query the number of filter scales available.
417 */
418 int scaleCount() const
419 { return scaleCount_; }
420
421 /** Query the number of filter directions available.
422 */
423 int directionCount() const
424 { return directionCount_; }
425
426 /** Change the number of directions / scales. This causes the
427 recalculation of all filters.
428 */
429 void setDirectionScaleCounts(int directionCount, int scaleCount)
430 {
431 this->resize(directionCount * scaleCount);
432 scaleCount_ = scaleCount;
433 directionCount_ = directionCount;
434 initFilters();
435 }
436
437 /** Return the center frequency of the filter(s) with
438 scale==0. Filters with scale>0 will have a center frequency
439 reduced in octaves:
440 <tt>centerFrequency= maxCenterFrequency / 2.0^scale</tt>
441 */
442 double maxCenterFrequency()
443 { return maxCenterFrequency_; }
444
445 /** Change the center frequency of the filter(s) with
446 scale==0. See \ref maxCenterFrequency().
447 */
448 void setMaxCenterFrequency(double maxCenterFrequency)
449 {
450 maxCenterFrequency_ = maxCenterFrequency;
451 initFilters();
452 }
453};
454
455//@}
456
457} // namespace vigra
458
459#endif // VIGRA_GABORFILTER_HXX
vigra::Diff2D
Two dimensional difference vector.
Definition diff2d.hxx:186
vigra::GaborFilterFamily::scaleCount
int scaleCount() const
Definition gaborfilter.hxx:418
vigra::GaborFilterFamily::setMaxCenterFrequency
void setMaxCenterFrequency(double maxCenterFrequency)
Definition gaborfilter.hxx:448
vigra::GaborFilterFamily::maxCenterFrequency
double maxCenterFrequency()
Definition gaborfilter.hxx:442
vigra::GaborFilterFamily::GaborFilterFamily
GaborFilterFamily(const Diff2D &size, int directionCount=stdDirectionCount, int scaleCount=stdScaleCount, double maxCenterFrequency=3.0/8.0, Alloc const &alloc=Alloc())
Definition gaborfilter.hxx:355
vigra::GaborFilterFamily::resizeImages
virtual void resizeImages(const Diff2D &newSize)
Definition gaborfilter.hxx:410
vigra::GaborFilterFamily::getFilter
ImageType const & getFilter(int direction, int scale) const
Definition gaborfilter.hxx:403
vigra::GaborFilterFamily::filterIndex
int filterIndex(int direction, int scale) const
Definition gaborfilter.hxx:392
vigra::GaborFilterFamily::setDirectionScaleCounts
void setDirectionScaleCounts(int directionCount, int scaleCount)
Definition gaborfilter.hxx:429
vigra::GaborFilterFamily::GaborFilterFamily
GaborFilterFamily(int width=stdFilterSize, int height=-1, int directionCount=stdDirectionCount, int scaleCount=stdScaleCount, double maxCenterFrequency=3.0/8.0, Alloc const &alloc=Alloc())
Definition gaborfilter.hxx:371
vigra::GaborFilterFamily::directionCount
int directionCount() const
Definition gaborfilter.hxx:423
vigra::ImageArray< ImageType, Alloc >::resizeImages
virtual void resizeImages(const Diff2D &newSize)
Definition imagecontainer.hxx:419
vigra::ImageArray::size
size_type size() const
Definition imagecontainer.hxx:229
vigra::ImageArray::resize
void resize(size_type newSize)
Definition imagecontainer.hxx:311
vigra::ImageArray::ImageArray
ImageArray(unsigned int numImages, const Diff2D &imageSize, Alloc const &alloc=Alloc())
Definition imagecontainer.hxx:100
vigra::MultiArrayView
Base class for, and view to, MultiArray.
Definition multi_array.hxx:705
vigra::Size2D
Two dimensional size object.
Definition diff2d.hxx:483
vigra
vigra::createGaborFilter
void createGaborFilter(...)
Create a gabor filter in frequency space.
vigra::angularGaborSigma
double angularGaborSigma(int directionCount, double centerFrequency)
Calculate sensible angular sigma for given parameters.
Definition gaborfilter.hxx:288
vigra::radialGaborSigma
double radialGaborSigma(double centerFrequency)
Calculate sensible radial sigma for given parameters.
Definition gaborfilter.hxx:245
utilities.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.12.4