[Insight-users] FFTW library in ITK

[Jakub Bican]

Hi all,

i am trying to use FFTW library in ITK through
FFTWRealToComplexConjugateImageFilter class. I am using VS.NET compiler so i
downloaded several FFTW binaries from fftw.org site.

Then there were two kinds of errors, depending on the FFTW library:
1) if i used MinGW compiled lib or Intel compiled lib form fftw.org download
page, the compilation was OK, but when i executed the code, it crashed
inside a fftw execution function (during transforming the data).

2) if i tried to compile my own FFTW binaries using project files from
fftw.org. In this case, some symbols were not found by linker.


So i went into ITK's tests and compiled itkAlgorithmsTests4. First, it
wasn't able to link, but i found the error and the solution here:
Insight\Code\Algorithms\CMakeLists.txt#353:
        TARGET_LINK_LIBRARIES(itkAlgorithmsTests4 fftw3 fftw3f)
replace with:
        TARGET_LINK_LIBRARIES(itkAlgorithmsTests4)


Then, when compiled, all tests ran OK!!
So i examined case 1) and found following problem:
First i casted an image by CastImageFilter from some integral type to an
image of doubles. Output of a caster was passed to the FFTW
RealToComplexConjugate filter. The execution crashed inside fftw filter
execution if:
1) caster filter was not explicitly updated by Update() call before FFT
filter's Update(), AND
2) dimension of image is 3, AND
3) PixelType of image passed into fft filter is 'double'.

I have written an intensive testing code, it is attached. See the lines
marked by //* - those are the cases, that fails on my system. I tested that
with several different "harmless" filters (now it is set to AbsImageFilter)
instead of CastImageFilter. The results on my systems are:
1) dimension must be 3
2) any input pixel type
3) any output pixel type {float|double} (!! this differs from the previous
test on real data)
4) any filter
5) filter may not be explicitly updated before FFTW filter's Update() call


Please, do someone know, where is the problem?? It can be hacked by calling
Update() on a filter preceding FFTW filter each times before FFTW filter's
execution, but it is not very clear solution.

Thanks,
Jakub
-------------- next part --------------
PROJECT(ffterr)

FIND_PACKAGE(ITK)

IF(ITK_FOUND)
	INCLUDE(${USE_ITK_FILE})
ELSE(ITK_FOUND)
	MESSAGE(FATAL_ERROR
					"ITK not found. Please set ITK_DIR.")
ENDIF(ITK_FOUND)

SET(CMAKE_MODULE_PATH ${ITK_SOURCE_DIR}/CMake)

SET(USE_FFTW 1)
FIND_PACKAGE( FFTW )

ADD_EXECUTABLE(ffterr ffterr.cxx )
							 
TARGET_LINK_LIBRARIES(ffterr ITKCommon ITKBasicFilters ITKIO ITKAlgorithms)

-------------- next part --------------
#include <stdlib.h>

#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkAbsImageFilter.h" //some harmless filter
#include "itkFFTWRealToComplexConjugateImageFilter.h"

template < unsigned int Dimension,
           class IntegralPixelType,	
					 class RealPixelType >          
int
test_fft_caster(bool         useCaster    = true,
								bool         updateFilter = false,
								unsigned int Size					= 77 )
{

  typedef itk::Image< IntegralPixelType,  Dimension >   IntegralImageType;
  typedef itk::Image< RealPixelType,  Dimension >   RealImageType;
	typedef itk::CastImageFilter< IntegralImageType, RealImageType > IntegralToRealCastFilterType;
	typedef itk::AbsImageFilter< RealImageType, RealImageType > HarmlessFilterType;
	typedef itk::FFTWRealToComplexConjugateImageFilter<RealPixelType,Dimension> R2CImageFilter;


	IntegralToRealCastFilterType::Pointer caster = IntegralToRealCastFilterType::New();
	HarmlessFilterType::Pointer filter = HarmlessFilterType::New();


	if (useCaster)
	{
		//create image
		IntegralImageType::SizeType imageSize;
		IntegralImageType::IndexType imageIndex;
		IntegralImageType::RegionType region;
		IntegralImageType::Pointer img = IntegralImageType::New();
		for(unsigned int i = 0; i < Dimension; i++)
			{
			imageSize.SetElement(i,Size);
			imageIndex.SetElement(i,0);
			}
		region.SetSize(imageSize);
		region.SetIndex(imageIndex);
		img->SetLargestPossibleRegion(region);
		img->SetBufferedRegion(region);
		img->SetRequestedRegion(region);
		img->Allocate();

		//cast image from integral to real type
		caster->SetInput(img);
		
		if (updateFilter)
			caster->Update();
	}
	else
	{
		//create image
		RealImageType::SizeType imageSize;
		RealImageType::IndexType imageIndex;
		RealImageType::RegionType region;
		RealImageType::Pointer img = RealImageType::New();
		for(unsigned int i = 0; i < Dimension; i++)
			{
			imageSize.SetElement(i,Size);
			imageIndex.SetElement(i,0);
			}
		region.SetSize(imageSize);
		region.SetIndex(imageIndex);
		img->SetLargestPossibleRegion(region);
		img->SetBufferedRegion(region);
		img->SetRequestedRegion(region);
		img->Allocate();

		//pass image through some harmless filter
		filter->SetInput(img);
		
		if (updateFilter)
			filter->Update();
	}

	//fourier-transform image
	R2CImageFilter::Pointer fftFilter = R2CImageFilter::New();
	if (useCaster) {
		fftFilter->SetInput(caster->GetOutput());
	}
	else
	{
		fftFilter->SetInput(filter->GetOutput());
	}
	
	try 
		{
			fftFilter->Update();
		} 
	catch (...) 
		{
			return EXIT_FAILURE;
		} 

  return EXIT_SUCCESS;
}


int main( int argc, char * argv[] )
{
	unsigned nerr = 0;

	unsigned size = 77; //use some "clever" size - FFTW may have explicit routines for sizes like 32, etc

	/* Combinations:
				Parameter 1 - dimension         : {1|2|3}
				Parameter 2 - input pixel type  : {unsigned char|int|float|double}
				Parameter 3 - output pixel type : {double|float}
				Parameter 4 - use caster filter or some harmless filter?      : {true|false}
				Parameter 5 - update filter explicitly before fftw execution? : {true|false}
				Parameter 6 - size of data      : {1|2|3}
	*/

	if (test_fft_caster<3, unsigned char, double>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, int, double>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, float, double>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, double, double>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, unsigned char, float>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, int, float>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, float, float>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, double, float>(true, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<2, unsigned char, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, unsigned char, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, double>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, float>(true, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, unsigned char, double>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, int, double>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, float, double>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, double, double>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, unsigned char, float>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, int, float>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, float, float>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<3, double, float>(false, false, size) == EXIT_FAILURE) //*
		nerr++;

	if (test_fft_caster<2, unsigned char, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, unsigned char, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, double>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, float>(false, false, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, unsigned char, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, int, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, float, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, double, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, unsigned char, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, int, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, float, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, double, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, unsigned char, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, unsigned char, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, double>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, float>(true, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, unsigned char, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, int, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, float, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, double, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, unsigned char, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, int, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, float, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<3, double, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, unsigned char, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, unsigned char, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, int, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, float, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<2, double, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, double>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, unsigned char, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, int, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, float, float>(false, true, size) == EXIT_FAILURE)
		nerr++;

	if (test_fft_caster<1, double, float>(false, true, size) == EXIT_FAILURE)
		nerr++;


	if (nerr>0)
	{
		std::cout << nerr << " errors !!!" << std::endl;
		return EXIT_FAILURE;
	}
	else
	{
		std::cout << "Success !!!" << std::endl;
		return EXIT_SUCCESS;
	}
}