[Insight-users] FFTW library in ITK
Jakub Bican
jakub.bican at matfyz.cz
Sat May 21 07:57:00 EDT 2005
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;
}
}
More information about the Insight-users
mailing list