FFT mismatch ITK/matlab

Jonas · March 2, 2020, 11:28am

Hi,

I would like to know how to match the FFT implementation of matlab fftn function using ITK.

I tried both itkForwardFFTImageFilter alone and the PadIlmageFilter + ForwardFFTImageFilter. taken from https://itk.org/ITKExamples/src/Filtering/FFT/ComputeForwardFFT/Documentation.html

So far the results are different in both cases for a 3D image. I’m in the process of implementing a filter, and I’m worried about these coefficient differences. I couldn’t find explicit references of potential difference in the documentation nor on Mathwork…

thanks

phcerdan · March 2, 2020, 2:32pm

How big are the differences? Could you provide an example?

Also, try to compile ITK with ITK_USE_FFTWD and/or ITK_USE_FFTWF, and check if there are differences too.
I think Matlab uses the FFTW library.

Jonas · March 2, 2020, 3:54pm

I checked my last CMakeCache.txt. FFTW is already activated on my system.

//Use double precision fftw if found
ITK_USE_FFTWD:BOOL=ON

//Use single precision fftw if found
ITK_USE_FFTWF:BOOL=ON

For now I’m testing my code on a tiny volume of (20,20,20).

I’m sorry for the data format but this is what I get from the matlab prompt

As for the itk output it is in the text file.The code used is the one from the FFT example. The format is not easy to read once again… Emacs in csv-mode gives a bit more readability. Results are printed using simple nested loops.

itk_fft_examples.txt (230.3 KB)

phcerdan · March 2, 2020, 5:06pm

How are you comparing the results then? You are providing a screenshot and a csv file. I am afraid this is unreadable.
You can visualize images in the frequency domain with any visualization software. Please upload the input image and the frequency images (after FFT) generated with both, Matlab and ITK, and point to C++ ITK script used to generate it.

Jonas · March 3, 2020, 3:21pm

Sorry for the late reply,

Here are more infos :
For matlab, I’m using the tools for reading NIFTI/Analyse images, which is annoying because it switches axes in the matrix transform compared to the itk implementation. But I did manage to have equivalent axis now.

Here is the first slice visualization on slicer of the real part image after FFT with automatic range for grey level. imaginary parts are also different.

testFFT.cpp (3.7 KB)
tinyCylinder.nii (31.6 KB)
itk_fft_imagPart.nii (31.6 KB)
itk_fft_realPart.nii (31.6 KB)
matlab_fftReal_flippedLRTD.nii (31.6 KB)
fftImag_flippedLRTD.nii (31.6 KB)

phcerdan · March 4, 2020, 12:07pm

Thanks for providing more details to reproduce it.

Please post the minimal set of matlab commands you are using to read the image, apply the fftn, and write the output image, with the switch of axis you commented.
The output image generated by a fft has positive and negative frequencies, and also Nyquist frequencies depending on the oddness in size of the axis, so after switching of axis, the output has to be interpreted with care.

Jonas · March 4, 2020, 1:19pm

data = load_nii(‘ImageName.nii’);

fftImg = fftn(data.img);

realPart = make_nii(real(fftmg),[1 1 1]);
imagPart = make_nii(imag(fftImg),[1 1 1]);

save_nii(realPart,‘fftRealPart.nii’);
save_nii(imagPart,‘fftImagPart.nii’);

flip_lrtd('fftRealPart.nii",“fftRealPart_flipped_lrtd”);
flip_lrtd('fftImagPart.nii",“fftImagPart_flipped_lrtd”);

Here are the infos on the toolbox.

Last information that might be useful, ITK transformation matrix (displayed by slicer in the volume module) is [-1 0 0; 0 -1 0; 0 0 1]. The matlab transformation matrix generated before the flip is [1 0 0; 0 1 0; 0 0 1].

phcerdan · March 4, 2020, 4:06pm

I don’t know what kind of transformation that plugin is doing with the function load_nii, but it doesn’t seem right, some comments in that plugin page seems to share this opinion.

You get the same results in ITK and Matlab if you use the load_untouch_nii function instead, which seems to not perform any extra magic.

fftReal_untouch_tinyCylinder.nii (31.6 KB) itk_fftReal_tinyCylinder_nowrappad.nii (31.6 KB)

Also, I don’t think you are using the WrapPadImageFilter in your itk script correctly. Use the itkFFTPadImageFilter instead.
I have used no pad filter (to match your matlab usage of fftn) and both, the VNL fft and the FFTW fft in ITK, with same results.

I am sure it is an issue in how that plugin is reading/wrtiting the image in matlab. Not in the ITK side.

Jonas · March 4, 2020, 5:26pm

I was going to ask about those prime numbers, the link is very welcomed.

I could replicate the results of your answer. Thank you very much !

Unfortunately this toolbox is the only one available to read nifti with a matlab version inferior to 2019 which we don’t have here…

I will be able to continue my work. Thanks again for the time you spent on this.

Jonas · March 6, 2020, 4:09pm

Just a quick update.

I was working on making a C++ version of Optimally oriented flux. The algorithm provided by Benmansour http://insight-journal.org/browse/publication/885 had some weird artefacts that I couldn’t explain. It was probably my attempt at updating the filter from itk V.3 ish to 5.1 that introduced weird things (middle view of the image).

Following the original matlab implementation from Law, seems much better (left view). Thanks again for the help.

A last question, is there an ITK implementation of Bessels functions with float parameters ? VNL bessels implementation only have int parameters. I would like to avoid using Boost::math for this filter.

PS: The Insight Journal comments are full of spam bot messages, it is becoming worst and worst.

blowekamp · March 6, 2020, 4:15pm

C++17 now has several Bessel function:
https://en.cppreference.com/w/cpp/numeric/special_functions

Do any of them meet your needs?

Corrected: Those are C++17 functions not C++11. ITK currently only requires C++11, but perhaps for your application C++17 as a requirement would be OK.

phcerdan · March 6, 2020, 4:53pm

Unrelated to the original topic, please open a new one.