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

| Fast Fourier Transform |

## Namespaces

namespace | vigra::detail

## Functions

template<... >
void | applyFourierFilter (...)
| Apply a filter (defined in the frequency domain) to an image.
template<... >
void | applyFourierFilterFamily (...)
| Apply an array of filters (defined in the frequency domain) to an image.
template<... >
void | fourierTransform (...)
| Compute forward and inverse Fourier transforms.
void | fourierTransformInverse (FFTWComplexImage::const_traverser sul, FFTWComplexImage::const_traverser slr, FFTWComplexImage::ConstAccessor src, FFTWComplexImage::traverser dul, FFTWComplexImage::Accessor dest)
| Compute inverse Fourier transforms.
template<... >
void | fourierTransformReal (...)
| Real Fourier transforms for even and odd boundary conditions (aka. cosine and sine transforms).
template<... >
void | moveDCToCenter (...)
| Rearrange the quadrants of a Fourier image so that the origin is in the image center.
template<... >
void | moveDCToUpperLeft (...)
| Rearrange the quadrants of a Fourier image so that the origin is in the image's upper left.

## Detailed Description

This documentation describes the VIGRA interface to FFTW version 3. The interface to the old FFTW version 2 (file "vigra/fftw.hxx") is deprecated.

VIGRA uses the FFTW Fast Fourier Transform package to perform Fourier transformations. VIGRA provides a wrapper for FFTW's complex number type (FFTWComplex), but FFTW's functions are used verbatim. If the image is stored as a FFTWComplexImage, the simplest call to an FFT function is like this:

        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 (this can be repeated as often as needed,
        //                with fresh data written into the source array)
        fftw_execute(forwardPlan);

        // release the plan memory
        fftw_destroy_plan(forwardPlan);

        // likewise for the inverse transform
        fftw_plan backwardPlan = fftw_plan_dft_2d(height, width,
                                     (fftw_complex *)fourier.begin(), (fftw_complex *)spatial.begin(),
                                     FFTW_BACKWARD, FFTW_ESTIMATE);
        fftw_execute(backwardPlan);
        fftw_destroy_plan(backwardPlan);

        // do not forget to normalize the result according to the image size
        transformImage(srcImageRange(spatial), destImage(spatial),
                       std::bind1st(std::multiplies<FFTWComplex>(), 1.0 / width / height));

Note that in the creation of a plan, the height must be given first. Note also that `spatial.begin()` may only be passed to `fftw_plan_dft_2d` if the transform shall be applied to the entire image. When you want to restrict operation to an ROI, you can create a copy of the ROI in an image of appropriate size, or you may use the Guru interface to FFTW.

More information on using FFTW can be found here.

FFTW produces fourier images that have the DC component (the origin of the Fourier space) in the upper left corner. Often, one wants the origin in the center of the image, so that frequencies always increase towards the border of the image. This can be achieved by calling moveDCToCenter(). The inverse transformation is done by moveDCToUpperLeft().

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

## Function Documentation

template<... >
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 transformation must not be executed in place - input and output images must be different.

Declarations:

pass arguments explicitly:

        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 :

        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>
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);

template<... >
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 declarations and a usage example.

Declarations:

pass arguments explicitly:

            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 :

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

template<... >
void vigra::fourierTransform | ( |  | ... | ) |

Compute forward and inverse Fourier transforms.

In the forward direction, the input image may be scalar or complex, and the output image is always complex. In the inverse direction, both input and output must be complex.

Declarations:

pass arguments explicitly:

        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 :

        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>
Namespace: vigra

        // compute complex Fourier transform of a real image
        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);

        fourierTransform(srcImageRange(fourier), destImage(inverseFourier));

void vigra::fourierTransformInverse | ( | FFTWComplexImage::const_traverser | sul,
|  | FFTWComplexImage::const_traverser | slr,
|  | FFTWComplexImage::ConstAccessor | src,
|  | FFTWComplexImage::traverser | dul,
|  | FFTWComplexImage::Accessor | dest |
| ) |  |  |

Compute inverse Fourier transforms.

See fourierTransform() for details.

template<... >
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 on).

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

Declarations:

pass arguments explicitly:

        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.

template<... >
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 arguments explicitly:

        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);

template<... >
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 arguments explicitly:

        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));

© 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 (3 Dec 2010)