[Insight-users] I have a problem With FFT in ITK

[Luis Ibanez]

Hi Simon,

We know that a problem exists at the level of using FFT filters
in ITK that are based on the VNL ones. However, we didn't get to
find the source of the problem.

It may well be that the ITK FFT filters are missing some operation
that is required for the good use of the vnl FFT; or it may be an
actual bug in the VNL implementation.

The observation is that when you take an image, computes its FFT
and then apply the inverse FFT, you don't recover the original
image. The returned image is structured in different blocks and
it has non-zero values in places where it should be zero.

The simplest way to test this is with the examples:

    Insight/Examples/Filtering/
                  FFTDirectInverse.cxx
                  FFTDirectInverse2.cxx


The first one uses VNL,
while the second one uses FFTW.


If you have any suggestions on what may be the source
of the problem we will be happy to give them a try.


Please let us know,


      Thanks



         Luis




--------------------
Simon Warfield wrote:
> Luis Ibanez wrote:
> 
>>
>> Hi Vahid,
>>
>>
>> Yes, this is a known problem with the VNL implementation of FFT.
> 
> 
> Luis,
> 
>  Is it known what the problem is ?
> 
>  Does VNL simply produce wrong results ?
>  Or is it some common difference in convention e.g. not apply a 1/N 
> scaling to the inverse FFT ?
> 
> -- 
> Simon
> 
>>
>>
>> You may want to install FFTW in your system and rerun CMake in
>> ITK in order to specify that you want to use FFTW. This is
>> currently an advanced option in the CMake configuration of ITK.
>>
>>
>>   Regards,
>>
>>
>>      Luis
>>
>>
>>
>> ---------------------------
>> vahid taimouri wrote:
>>
>>> Hi all
>>> I tried to use FFT to build an ideal lowpass filter
>>> But I think FFT doesn't work correctly,
>>> Also I applied the FFT to an image and then the inverse FFT to its 
>>> result
>>> But surprisingly it didn't retrieve my image!!!
>>>  
>>> I added my codes here, please guild to solve this problem
>>> Thanks
>>>  
>>>  void main()
>>> {   char InputImageName[30] = "Inputimage.bmp";
>>>  char OutputImageName[30] = "OutputImage.bmp";
>>>  
>>>  typedef itk::Image< unsigned char, 2 > InputImageType;
>>>  typedef itk::Image< unsigned char, 2 > OutputImageType;
>>>  
>>>  typedef itk::Image< float, 2 > InputFilterType;
>>>  typedef itk::Image< float, 2 > OutputFilterType;
>>>  
>>>  typedef itk::ImageFileReader< InputImageType > ReaderType;
>>>  typedef itk::RescaleIntensityImageFilter<InputImageType, 
>>> InputFilterType > RescaleFilterType;
>>>  typedef itk::VnlFFTRealToComplexConjugateImageFilter< float, 2 > 
>>> FFTFilterType;
>>>  typedef itk::VnlFFTComplexConjugateToRealImageFilter< float, 2 > 
>>> IFFTFilterType;
>>>  typedef itk::RescaleIntensityImageFilter<InputFilterType, 
>>> OutputImageType > RescaleFilterType;
>>>  typedef itk::ImageFileWriter< OutputImageType > WriterType;
>>>  ReaderType::Pointer reader = ReaderType::New();
>>>  InputRescaleFilterType::Pointer inputintensityrescaler = 
>>> InputRescaleFilterType::New();
>>>  FFTFilterType::Pointer fftFilter = FFTFilterType::New(); 
>>>  IFFTFilterType::Pointer fftInverseFilter = IFFTFilterType::New();
>>>  OutputRescaleFilterType::Pointer outputintensityrescaler = 
>>> OutputRescaleFilterType::New();
>>>  WriterType::Pointer writer = WriterType::New();
>>>  
>>>  Inputintensityrescaler->SetOutputMinimum(  0  );
>>>  Inputintensityrescaler->SetOutputMaximum( 255 );
>>>  
>>>  Outputintensityrescaler->SetOutputMinimum(  0  );
>>>  Outputintensityrescaler->SetOutputMaximum( 255 );
>>>  
>>>  reader->SetFileName( InputImageName );
>>>  writer->SetFileName( OutputImageName );
>>>  reader->Update();
>>>  
>>>  inputintensityrescaler->SetInput( reader->GetOutput() );
>>>  fftFilter->SetInput( inputintensityrescaler->GetOutput() );
>>>  fftFilter->Update();
>>>  
>>>  fftInverseFilter->SetInput( fftFilter->GetOutput() );
>>>  
>>>  typedef itk::FlipImageFilter< InputFilterType > FlipperType;
>>>  FlipperType::Pointer flipper = FlipperType::New();
>>>  bool fliparray[2] = {true,true};
>>>  FlipperType::FlipAxesArrayType flipAxes( fliparray );
>>>  flipper->SetFlipAxes(flipAxes);
>>>  flipper->SetInput(fftInverseFilter->GetOutput());
>>>
>>>  
>>>  Outputintensityrescaler->SetInput( flipper->GetOutput() );
>>>
>>>  writer->SetInput(Outputintensityrescaler->GetOutput());
>>>  writer->Update();
>>> }
>>>
>>> ------------------------------------------------------------------------
>>> Do you Yahoo!?
>>> Everyone is raving about the all-new Yahoo! Mail Beta. 
>>> <http://us.rd.yahoo.com/evt=42297/*http://advision.webevents.yahoo.com/handraisers> 
>>>
>>>
>>>
>>> ------------------------------------------------------------------------
>>>
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>
>>
>>
>> _______________________________________________
>> Insight-users mailing list
>> Insight-users at itk.org
>> http://www.itk.org/mailman/listinfo/insight-users
> 
> 
>