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

details Fast Fourier Transform VIGRA

Classes

class  FFTWConvolvePlan< N, Real >
 
class  FFTWCorrelatePlan< N, Real >
 
class  FFTWPlan< N, Real >
 

Functions

template<... >
void applyFourierFilter (...)
 Apply a filter (defined in the frequency domain) to an image. More...
 
template<... >
void applyFourierFilterFamily (...)
 Apply an array of filters (defined in the frequency domain) to an image. More...
 
template<class T , int N>
TinyVector< T, N > fftwCorrespondingShapeR2C (TinyVector< T, N > shape)
 Find frequency domain shape for a R2C Fourier transform. More...
 
template<... >
void fourierTransform (...)
 Compute forward and inverse Fourier transforms. More...
 
template<... >
void fourierTransformInverse (...)
 Compute inverse Fourier transforms. More...
 
template<... >
void fourierTransformReal (...)
 Real Fourier transforms for even and odd boundary conditions (aka. cosine and sine transforms). More...
 
template<... >
void moveDCToCenter (...)
 Rearrange the quadrants of a Fourier image so that the origin is in the image center. More...
 
template<... >
void moveDCToUpperLeft (...)
 Rearrange the quadrants of a Fourier image so that the origin is in the image's upper left. More...
 

Detailed Description

VIGRA provides a powerful C++ API for the popular FFTW library for fast Fourier transforms. There are two versions of the API: an older one based on image iterators (and therefore restricted to 2D) and a new one based on MultiArrayView that works for arbitrary dimensions. In addition, the functions convolveFFT() and applyFourierFilter() provide an easy-to-use interface for FFT-based convolution, a major application of Fourier transforms.

Function Documentation

void vigra::moveDCToCenter (   ...)

Rearrange the quadrants of a Fourier image so that the origin is in the image center.

FFTW produces fourier images where the DC component (origin of fourier space) is located in the upper left corner of the image. The quadrants are placed like this (using a 4x4 image for example):

DC 4 3 3
4 4 3 3
1 1 2 2
1 1 2 2

After applying the function, the quadrants are at their usual places:

2 2 1 1
2 2 1 1
3 3 DC 4
3 3 4 4

This transformation can be reversed by moveDCToUpperLeft(). Note that the 2D versions of this transformation must not be executed in place - input and output images must be different. In contrast, the nD version (with MultiArrayView argument) always works in-place.

Declarations:

use MultiArrayView (this works in-place, with arbitrary dimension N):

namespace vigra {
template <unsigned int N, class T, class Stride>
void moveDCToCenter(MultiArrayView<N, T, Stride> a);
}

pass iterators explicitly (2D only, not in-place):

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class DestImageIterator, class DestAccessor>
void moveDCToCenter(SrcImageIterator src_upperleft,
SrcImageIterator src_lowerright, SrcAccessor sa,
DestImageIterator dest_upperleft, DestAccessor da);
}

use argument objects in conjunction with Argument Object Factories (2D only, not in-place):

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class DestImageIterator, class DestAccessor>
void moveDCToCenter(
triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
pair<DestImageIterator, DestAccessor> dest);
}

Usage:

#include <vigra/fftw3.hxx> (for 2D variants)
#include <vigra/multi_fft.hxx> (for nD variants)
Namespace: vigra

vigra::FFTWComplexImage spatial(width,height), fourier(width,height);
... // fill image with data
// create a plan with estimated performance optimization
fftw_plan forwardPlan = fftw_plan_dft_2d(height, width,
(fftw_complex *)spatial.begin(), (fftw_complex *)fourier.begin(),
FFTW_FORWARD, FFTW_ESTIMATE );
// calculate FFT
fftw_execute(forwardPlan);
vigra::FFTWComplexImage rearrangedFourier(width, height);
moveDCToCenter(srcImageRange(fourier), destImage(rearrangedFourier));
// delete the plan
fftw_destroy_plan(forwardPlan);
void vigra::moveDCToUpperLeft (   ...)

Rearrange the quadrants of a Fourier image so that the origin is in the image's upper left.

This function is the inversion of moveDCToCenter(). See there for a detailed description and usage examples.

Declarations:

use MultiArrayView (this works in-place, with arbitrary dimension N):

namespace vigra {
template <unsigned int N, class T, class Stride>
void moveDCToUpperLeft(MultiArrayView<N, T, Stride> a);
}

pass iterators explicitly (2D only, not in-place):

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class DestImageIterator, class DestAccessor>
void moveDCToUpperLeft(SrcImageIterator src_upperleft,
SrcImageIterator src_lowerright, SrcAccessor sa,
DestImageIterator dest_upperleft, DestAccessor da);
}

use argument objects in conjunction with Argument Object Factories (2D only, not in-place):

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class DestImageIterator, class DestAccessor>
void moveDCToUpperLeft(
triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
pair<DestImageIterator, DestAccessor> dest);
}
void vigra::fourierTransform (   ...)

Compute forward and inverse Fourier transforms.

The array referring to the spatial domain (i.e. the input in a forward transform, and the output in an inverse transform) may be scalar or complex. The array representing the frequency domain (i.e. output for forward transform, input for inverse transform) must always be complex.

The new implementations (those using MultiArrayView arguments) perform a normalized transform, whereas the old ones (using 2D iterators or argument objects) perform an un-normalized transform (i.e. the result of the inverse transform is scaled by the number of pixels).

In general, input and output arrays must have the same shape, with the exception of the special R2C and C2R modes defined by FFTW.

The R2C transform reduces the redundancy in the Fourier representation of a real-valued signal: Since the Fourier representation of a real signal is symmetric, about half of the Fourier coefficients can simply be dropped. By convention, this reduction is applied to the first (innermost) dimension, such that fourier.shape(0) == spatial.shape(0)/2 + 1 holds. The correct frequency domain shape can be conveniently computed by means of the function fftwCorrespondingShapeR2C().

Note that your program must always link against libfftw3. If you want to compute Fourier transforms for float or long double arrays, you must additionally link against libfftw3f and libfftw3l respectively. (Old-style functions only support double).

The Fourier transform functions internally create FFTW plans which control the algorithm details. The plans are creates with the flag FFTW_ESTIMATE, i.e. optimal settings are guessed or read from saved "wisdom" files. If you need more control over planning, you can use the class FFTWPlan.

Declarations:

use complex-valued MultiArrayView arguments (this works for arbitrary dimension N):

namespace vigra {
template <unsigned int N, class Real, class C1, class C2>
void
fourierTransform(MultiArrayView<N, FFTWComplex<Real>, C1> in,
MultiArrayView<N, FFTWComplex<Real>, C2> out);
template <unsigned int N, class Real, class C1, class C2>
void
fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
MultiArrayView<N, FFTWComplex<Real>, C2> out);
}

use real-valued MultiArrayView in the spatial domain, complex-valued MultiArrayView in the frequency domain (this works for arbitrary dimension N, and also supports the R2C and C2R transform, depending on the array shape in the frequency domain):

namespace vigra {
template <unsigned int N, class Real, class C1, class C2>
void
fourierTransform(MultiArrayView<N, Real, C1> in,
MultiArrayView<N, FFTWComplex<Real>, C2> out);
template <unsigned int N, class Real, class C1, class C2>
void
fourierTransformInverse(MultiArrayView<N, FFTWComplex<Real>, C1> in,
MultiArrayView<N, Real, C2> out);
}

pass iterators explicitly (2D only, double only):

namespace vigra {
template <class SrcImageIterator, class SrcAccessor>
void fourierTransform(SrcImageIterator srcUpperLeft,
SrcImageIterator srcLowerRight, SrcAccessor src,
FFTWComplexImage::traverser destUpperLeft, FFTWComplexImage::Accessor dest);
void
fourierTransformInverse(FFTWComplexImage::const_traverser sul,
FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src,
FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
}

use argument objects in conjunction with Argument Object Factories (2D only, double only):

namespace vigra {
template <class SrcImageIterator, class SrcAccessor>
void fourierTransform(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
void
fourierTransformInverse(triple<FFTWComplexImage::const_traverser,
FFTWComplexImage::const_traverser, FFTWComplexImage::ConstAccessor> src,
pair<FFTWComplexImage::traverser, FFTWComplexImage::Accessor> dest);
}

Usage:

#include <vigra/fftw3.hxx> (old-style 2D variants)
#include <vigra/multi_fft.hxx> (new-style nD variants)
Namespace: vigra

old-style example (using FFTWComplexImage):

// compute complex Fourier transform of a real image, old-style implementation
vigra::DImage src(w, h);
vigra::FFTWComplexImage fourier(w, h);
fourierTransform(srcImageRange(src), destImage(fourier));
// compute inverse Fourier transform
// note that both source and destination image must be of type vigra::FFTWComplexImage
vigra::FFTWComplexImage inverseFourier(w, h);
fourierTransformInverse(srcImageRange(fourier), destImage(inverseFourier));

new-style examples (using MultiArray):

// compute Fourier transform of a real array, using the R2C algorithm
MultiArray<2, double> src(Shape2(w, h));
MultiArray<2, FFTWComplex<double> > fourier(fftwCorrespondingShapeR2C(src.shape()));
fourierTransform(src, fourier);
// compute inverse Fourier transform, using the C2R algorithm
MultiArray<2, double> dest(src.shape());
fourierTransformInverse(fourier, dest);
// compute Fourier transform of a real array with standard algorithm
MultiArray<2, double> src(Shape2(w, h));
MultiArray<2, FFTWComplex<double> > fourier(src.shape());
fourierTransform(src, fourier);
// compute inverse Fourier transform, using the C2R algorithm
MultiArray<2, double> dest(src.shape());
fourierTransformInverse(fourier, dest);

Complex input arrays are handled in the same way.

void vigra::fourierTransformInverse (   ...)

Compute inverse Fourier transforms.

See fourierTransform() for details.

void vigra::applyFourierFilter (   ...)

Apply a filter (defined in the frequency domain) to an image.

After transferring the image into the frequency domain, it is multiplied pixel-wise with the filter and transformed back. The result is put into the given destination image which must have the right size. The result will be normalized to compensate for the two FFTs.

If the destination image is scalar, only the real part of the result image is retained. In this case, you are responsible for choosing a filter image which ensures a zero imaginary part of the result (e.g. use a real, even symmetric filter image, or a purely imaginary, odd symmetric one).

The DC entry of the filter must be in the upper left, which is the position where FFTW expects it (see moveDCToUpperLeft()).

See also convolveFFT() for corresponding functionality on the basis of the MultiArrayView interface.

Declarations:

pass 2D array views:

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class FilterImageIterator, class FilterAccessor,
class DestImageIterator, class DestAccessor>
void applyFourierFilter(SrcImageIterator srcUpperLeft,
SrcImageIterator srcLowerRight, SrcAccessor sa,
FilterImageIterator filterUpperLeft, FilterAccessor fa,
DestImageIterator destUpperLeft, DestAccessor da);
}

pass Image Iterators and Data Accessors :

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class FilterImageIterator, class FilterAccessor,
class DestImageIterator, class DestAccessor>
void applyFourierFilter(SrcImageIterator srcUpperLeft,
SrcImageIterator srcLowerRight, SrcAccessor sa,
FilterImageIterator filterUpperLeft, FilterAccessor fa,
DestImageIterator destUpperLeft, DestAccessor da);
}

use argument objects in conjunction with Argument Object Factories :

namespace vigra {
template <class SrcImageIterator, class SrcAccessor,
class FilterImageIterator, class FilterAccessor,
class DestImageIterator, class DestAccessor>
void applyFourierFilter(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
pair<FilterImageIterator, FilterAccessor> filter,
pair<DestImageIterator, DestAccessor> dest);
}

Usage:

#include <vigra/fftw3.hxx>
Namespace: vigra

// create a Gaussian filter in Fourier space
vigra::FImage gaussFilter(w, h), filter(w, h);
for(int y=0; y<h; ++y)
for(int x=0; x<w; ++x)
{
xx = float(x - w / 2) / w;
yy = float(y - h / 2) / h;
gaussFilter(x,y) = std::exp(-(xx*xx + yy*yy) / 2.0 * scale);
}
// applyFourierFilter() expects the filter's DC in the upper left
moveDCToUpperLeft(srcImageRange(gaussFilter), destImage(filter));
vigra::FFTWComplexImage result(w, h);
vigra::applyFourierFilter(srcImageRange(image), srcImage(filter), result);

For inspection of the result, FFTWMagnitudeAccessor might be useful. If you want to apply the same filter repeatedly, it may be more efficient to use the FFTW functions directly with FFTW plans optimized for good performance.

void vigra::applyFourierFilterFamily (   ...)

Apply an array of filters (defined in the frequency domain) to an image.

This provides the same functionality as applyFourierFilter(), but applying several filters at once allows to avoid repeated Fourier transforms of the source image.

Filters and result images must be stored in vigra::ImageArray data structures. In contrast to applyFourierFilter(), this function adjusts the size of the result images and the the length of the array.

Declarations:

pass 2D array views:

namespace vigra {
template <class SrcImageIterator, class SrcAccessor, class FilterType>
void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
SrcImageIterator srcLowerRight, SrcAccessor sa,
const ImageArray<FilterType> &filters,
ImageArray<FFTWComplexImage> &results)
}

pass Image Iterators and Data Accessors :

namespace vigra {
template <class SrcImageIterator, class SrcAccessor, class FilterType>
void applyFourierFilterFamily(SrcImageIterator srcUpperLeft,
SrcImageIterator srcLowerRight, SrcAccessor sa,
const ImageArray<FilterType> &filters,
ImageArray<FFTWComplexImage> &results)
}

use argument objects in conjunction with Argument Object Factories :

namespace vigra {
template <class SrcImageIterator, class SrcAccessor, class FilterType>
void applyFourierFilterFamily(triple<SrcImageIterator, SrcImageIterator, SrcAccessor> src,
const ImageArray<FilterType> &filters,
ImageArray<FFTWComplexImage> &results)
}

Usage:

#include <vigra/fftw3.hxx>
Namespace: vigra

// assuming the presence of a real-valued image named "image" to
// be filtered in this example
vigra::ImageArray<vigra::FImage> filters(16, image.size());
for (int i=0; i<filters.size(); i++)
// create some meaningful filters here
createMyFilterOfScale(i, destImage(filters[i]));
vigra::ImageArray<vigra::FFTWComplexImage> results();
vigra::applyFourierFilterFamily(srcImageRange(image), filters, results);
void vigra::fourierTransformReal (   ...)

Real Fourier transforms for even and odd boundary conditions (aka. cosine and sine transforms).

If the image is real and has even symmetry, its Fourier transform is also real and has even symmetry. The Fourier transform of a real image with odd symmetry is imaginary and has odd symmetry. In either case, only about a quarter of the pixels need to be stored because the rest can be calculated from the symmetry properties. This is especially useful, if the original image is implicitly assumed to have reflective or anti-reflective boundary conditions. Then the "negative" pixel locations are defined as

even (reflective boundary conditions): f[-x] = f[x] (x = 1,...,N-1)
odd (anti-reflective boundary conditions): f[-1] = 0
f[-x] = -f[x-2] (x = 2,...,N-1)

end similar at the other boundary (see the FFTW documentation for details). This has the advantage that more efficient Fourier transforms that use only real numbers can be implemented. These are also known as cosine and sine transforms respectively.

If you use the odd transform it is important to note that in the Fourier domain, the DC component is always zero and is therefore dropped from the data structure. This means that index 0 in an odd symmetric Fourier domain image refers to the first harmonic. This is especially important if an image is first cosine transformed (even symmetry), then in the Fourier domain multiplied with an odd symmetric filter (e.g. a first derivative) and finally transformed back to the spatial domain with a sine transform (odd symmetric). For this to work properly the image must be shifted left or up by one pixel (depending on whether the x- or y-axis is odd symmetric) before the inverse transform can be applied. (see example below).

The real Fourier transform functions are named fourierTransformReal?? where the questions marks stand for either E or O indicating whether the x- and y-axis is to be transformed using even or odd symmetry. The same functions can be used for both the forward and inverse transforms, only the normalization changes. For signal processing, the following normalization factors are most appropriate:

forward inverse
------------------------------------------------------------
X even, Y even 1.0 4.0 * (w-1) * (h-1)
X even, Y odd -1.0 -4.0 * (w-1) * (h+1)
X odd, Y even -1.0 -4.0 * (w+1) * (h-1)
X odd, Y odd 1.0 4.0 * (w+1) * (h+1)

where w and h denote the image width and height.

Declarations:

pass 2D array views:

namespace vigra {
template <class SrcTraverser, class SrcAccessor,
class DestTraverser, class DestAccessor>
void
fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
DestTraverser dul, DestAccessor dest, fftw_real norm);
fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
}

pass Image Iterators and Data Accessors :

namespace vigra {
template <class SrcTraverser, class SrcAccessor,
class DestTraverser, class DestAccessor>
void
fourierTransformRealEE(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
DestTraverser dul, DestAccessor dest, fftw_real norm);
fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
}

use argument objects in conjunction with Argument Object Factories :

namespace vigra {
template <class SrcTraverser, class SrcAccessor,
class DestTraverser, class DestAccessor>
void
fourierTransformRealEE(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
pair<DestTraverser, DestAccessor> dest, fftw_real norm);
fourierTransformRealEO, fourierTransformRealOE, fourierTransformRealOO likewise
}

Usage:

#include <vigra/fftw3.hxx>
Namespace: vigra

vigra::FImage spatial(width,height), fourier(width,height);
... // fill image with data
// forward cosine transform == reflective boundary conditions
fourierTransformRealEE(srcImageRange(spatial), destImage(fourier), (fftw_real)1.0);
// multiply with a first derivative of Gaussian in x-direction
for(int y = 0; y < height; ++y)
{
for(int x = 1; x < width; ++x)
{
double dx = x * M_PI / (width - 1);
double dy = y * M_PI / (height - 1);
fourier(x-1, y) = fourier(x, y) * dx * std::exp(-(dx*dx + dy*dy) * scale*scale / 2.0);
}
fourier(width-1, y) = 0.0;
}
// inverse transform -- odd symmetry in x-direction, even in y,
// due to symmetry of the filter
fourierTransformRealOE(srcImageRange(fourier), destImage(spatial),
(fftw_real)-4.0 * (width+1) * (height-1));
TinyVector<T, N> vigra::fftwCorrespondingShapeR2C ( TinyVector< T, N >  shape)

Find frequency domain shape for a R2C Fourier transform.

When a real valued array is transformed to the frequency domain, about half of the Fourier coefficients are redundant. The transform can be optimized as a R2C transform that doesn't compute and store the redundant coefficients. This function computes the appropriate frequency domain shape for a given shape in the spatial domain. It simply replaces shape[0] with shape[0] / 2 + 1.

#include <vigra/multi_fft.hxx>
Namespace: vigra

© 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.11.1 (Fri May 19 2017)