public class AlgorithmFFT2
extends AlgorithmBase

Processing images by filtering in the frequency domain is a 3 step process: 1.) Performing a forward fast fourier transform to convert a spatial image into its frequency image. 2.) Enhancing some frequency components of the image and attenuating other frequency components of the image by mulitplying by a lowpass, highpass, bandpass, or bandstop filter. Frequency filters may be constructed with 1 of 3 methods - finite impulse response filters constructed with Hamming windows, Gaussian filters, and Butterworth filters. However, for the Gaussian filters only lowpass and highpass filters are available. 3.) Performing an inverse fast fourier transform to convert from the frequency domain back into the spatial domain. This software module performs all 3 steps of this process, but only performs one step after a dialog OK. The AlgorithmFrequencyFilter module performs all 3 steps after a dialog OK.

The module also creates a Gabor filter, which is essentially a tilted Gaussian with 2 unequal axes at an offset (freqU, freqV) from the origin. A Gabor filter only responds to a texture having both a particular frequency and a particular orientation. Note that a filter and its mirror image reflected across the u and v frequency axes produce identical frequency responses.

The core algorithm of this module, the fast fourier transform algorithm found in exec(), requires that all the dimensions of an N-dimensional dataset be powers of 2. To be able to use this algorithm on datasets with arbitrary dimensions, the data is zero padded to powers of 2 before applying the forward fast fourier transform and stripped down to the original dimensions after applying the inverse fast fourier transform. If unequalDim is equal to true, the fourier pictures are allowed to have unequal dimensions so as to save memory. If unequalDim is equal to false, the pictures have equal dimensions so as to maximize symmetry.

The typical full sequence is as follows: 1.) Data from a real spatial image is exported into a double array realData. 2.) An equally sized double array called imagData is created and filled with zeros. 3.) realData and imagData are enlarged to the same length in every dimension. If finite impulse repsonse filters are constructed with Hamming windows and no cropping, the new dimension size is equal to the minimum power of two number that equals or exceeds the maximum original dimension size + kDim - 1, where kDim is the diameter of a circular or spherical convolution kernel. If finite impulse response filters with Hamming windows and cropping or infinite impulse response Gaussian or Butterworth filters are used, the new dimension size is equal to the minimum power of two number that equals or exceeds the maximum original dimension size. The data is padded with zeros at the end of each dimension. 4.) exec() is invoked to run the fast fourier transform algorithm. 5.) The center() algorithm is invoked to reorder the data for display. The Fourier transform of real data is conjugate symmetric; the real parts are even and the imaginary parts are odd. In other words, the magnitude is even and the phase is odd. Since images display the magnitude or log10(1 + magnitude) a symmetry reflected across the center of the image is observed. 6.) The complex data is imported to the display image with the image minimum and maximum being calculated. If logMagDisplay is set equal to false, the minimum and maximum of sqrt(realData(i)*realData(i) + imagData(i)*imagData(i)) will be used as the image minimum and maximum in other modules. If logMagDisplay is set equal to true, the log10(1 + sqrt(realData(i)*realData(i) + imagData(i)*imagData(i))) will be used as the image minimum and maximum in other modules. Note that while this module stores the complex data in the source or destination image, the display oriented modules in the MIPAV package will use the magnitude or log magnitude as indicated above in the screen displays since complex numbers cannot be directly displayed. 7.) The command setOriginalExtents is invoked so that the original dimensions of the source image are available along with the display image. The original dimensions are needed so that the original source dimensions can be restored after the inverse fast fourier transform. If finite impulse response filters with windows are used: 8a.) An ideal filter kernel is constructed in the spatial domain. 9a.) A Hamming window kernel is constructed in the spatial domain. 10a.) The ideal kernel and the Hamming window kernel are multiplied together. 11a.) The kernel is zero padded up to the same dimensions that the image data was. 12a.) exec() is run to obtain the FFT of the kernel. 13a.) The center() algorithm is run to reorder the data. If Gaussian or Butterworth filters are used: For Gaussian and Butterworth filters the transfer functions affect the real and imaginary parts of the FFT of the image in exactly the same manner. These filters are zero-phase- shift filters because they do not alter the phase of the transform. Since these filters are conjugate symmetric, the inverse FFTs of these filters are purely real. Since this filtering is equivalent to convolving to real 2D data sets, the result must be purely real. 8b.) The real part of the FFT is set equal to the appropriate filter magnitude and the imaginary part is set equal to zero. This filter has the same dimensions as the padded image data. 14) The data FFT is set equal to the product of the data FFT and the filter FFT. 15.) The inverse FFT process is invoked. The complex data is exported into the 2 double arrays realData and imagData. There should be no need for zero padding at this point since the dimensions should already be all powers of 2 from before. 16.) The center() algorithm is invoked to restore the data to its original ordering. 17.) exec() is invoked to run the inverse fast fourier transform algorithm. 18.) The realData now holds the correct response if Gaussian or Butterworth filtering was used. If FIR filtering with windows was used then realData holds a version of the correct response shifted by (kDim - 1)/2 toward the end of each dimension. The imagData should contain only roundoff error. 19.) For FIR filters shift the data back by (kDim - 1)/2 toward the start of each dimension. 20.) Stripping is performed to return the image to its original dimensions. 21.) The realData() is imported into the the new spatial image.

Methods included calculate the magnitude and phase as shown below:

magnitude = ( (realData)^2 + (imagData)^2 )^(1/2); phase = arctan(imagData/realData);