-
Notifications
You must be signed in to change notification settings - Fork 1
FourierTransformations
Fourier stuff is distributed in 3 files:
- xmipp_fftw.h for those things related to fftw (external library for Fourier transforms, and preferred implementation)
- xmipp_fft.h for those things related to frequency conversions, and the old implementation of Fourier transform.
- fourier_filter.h for those things related to filters.
In Xmipp we have the concept of a FourierTransformer. It is an object that keeps a pointer to the image in real space (the image can be real-valued or complex valued), and keeps an internal copy of the Fourier transform (MultidimArray of complex values). The most simple use of the Fourier transformer is the following (for instance, to compute the magnitude of the Fourier transform of a volume (the syntax for an image is exactly the same).
// Read volume Image V; V.read("myVolume.vol");
// Compute the Fourier transform FourierTransformer transformer; MultidimArray< std::complex > Vfft; transformer.FourierTransform(V(),Vfft);
// Compute magnitude MultidimArray Vmag; Vmag.resizeNoCopy(Vfft); FOR_ALL_ELEMENTS_IN_ARRAY3D(Vmag) A3D_ELEM(Vmag,k,i,j)=20*log10(abs(A3D_ELEM(Vfft,k,i,j)));
You can use the same transformer to compute several Fourier transforms. This is rather efficient if all input volumes are of the same size. This is because before computing a Fourier Transform, FFTW establish a plan for a computationally optimal algorithm. If the size changes, the plan has to be recomputed (all this is internally handled by the class FourierTransformer). If the size does not change, the plan is not recomputed and the Fourier transform is computed straightaway.
Image V; FourierTransformer transformer; MultidimArray< std::complex > Vfft; for (int i=0; i<10; i++) {
- read(formatString("volume%02d.vol",i)); // volume01.vol, volume02.vol, ... transformer.FourierTransform(V(),Vfft);
}
To perform an inverse Fourier transform you simply do:
Note that transformer.FourierTransform(V(),Vfft); makes a copy of the internal object for the Fourier transform of V into the local variable Vfft. If you do not want to pay for this extra copy, but you want that Vfft shares the same memory as the internal Fourier transformer, then you can do:
You can also use the transformer without explicitly setting the input and output every time as long as the memory assigned to the input object does not change.
Image V; V.initZeros(64,64,64);FourierTransformer transformer;MultidimArray< std::complex > Vfft; transformer.setReal(V()); for (int i=0; i<10; i++) { // Change V depending on your algorithm ... // Compute the Fourier transform transformer.FourierTransform(); transformer.getFourierAlias(Vfft); // Vfft is sharing the memory of the transformer
// Change the Fourier transform according to your algorithm ...
// Compute the inverse Fourier transform transformer.inverseFourierTransform(); // The output is in V() }
Normally, FFTW only computes one half of the spectrum for real valued inputs (because of the Hermitian symmetry). You may use getCompleteFourier to get the full Fourier transform. In any case, the Fourier transform is computed in the range from [0,1) or [0,2*pi) (depending on how you normalize digital frequencies; MATLAB does the same, the concept is similar to what MATLAB produces before doing fftshift). The first value of the Fourier transform corresponds to the frequency origin (0,0) (i.e., the sum of all pixels in the image). A typical loop is the following:
Image I; I.initZeros(64,64);FourierTransformer transformer;MultidimArray< std::complex > Ifft; transformer.FourierTransform(I(),Ifft);
Matrix1D idx(2); // Indexes in Fourier planeMatrix1D digitalFreq(2); // Corresponding digital frequency between 0 and 1 (cycles/pixel)Matrix1D continuousFreq(2); // Corresponding continuous frequency (Angstroms/pixel) const double Tm=1.5; // Sampling rate in Angstroms per pixel FOR_ALL_ELEMENTS_IN_ARRAY2D(Ifft) { XX(idx) = j; YY(idx) = i;
// Digital frequency FFT_idx2digfreq(I(), idx, digitalFreq);
// Continuous frequency digfreq2contfreq(digitalFreq, continuousFreq, Tm); }
All functions for index<->frequency conversions are in xmipp_fft.h.
Fourier filters are implemented through the class FourierFilter. You need to include <reconstruction/fourier_filter.h>.
This class can apply Fourier filters whose mask in Fourier space is binary or continuous, 2D or 3D. Normally the mask is not explicitly created in memory. For instance, for applying a low pass filter, the mask is calculated on the fly during the application of the filter to any image. However, for certain filter types considered to be difficult on the fly, the mask is explicitly created in memory before applying it. Then, during the application, the mask value is simply read from this internal variable. In the following we see a couple of examples of filter application.
Application of a Highpass filter to an image
Image I; I.read("image.xmp"); FourierFilter Filter;Filter.FilterBand=HIGHPASS; Filter.w1=0.05; Filter.raised_w=0.02; I().setXmippOrigin(); Filter.applyMaskSpace(I()); I.write("filtered_image.xmp");
Application of a Wedge filter to a volume