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

Namespaces | Functions
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

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

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.

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);
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.0 (Thu Aug 25 2011)