[ITK-users] 1D Complex To Complex FFT of a 2D Image

[DVigneault]

All--

I would like to compute a 1D complex to complex FFT on each line of a 2D
image.  I'm working from  Luis's response to a similar question
<http://www.itk.org/pipermail/insight-users/2007-May/022250.html>  ,  this
basic example from the FFTW manual
<http://www.fftw.org/fftw3_doc/Complex-One_002dDimensional-DFTs.html#Complex-One_002dDimensional-DFTs> 
, and the FFTWComplexToComplexImageFilter.hxx source code.  However, when I
try to assign the plan, I get a segmentation fault.

So, my questions are:

1.  Is there an existing filter or method for doing this in ITK that I've
missed?
2.  If there is not a filter, is there a better way of approaching the
problem than I've laid out in the code below?
3.  If not, can anyone see why the code is seg faulting?

Thank you all very much in advance!

###########
## The Code ##
###########

// Headers
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <fftw3.h>
#include
"/Users/Davis/Developer/ITK/dv-headers/dvVisualizeComplexSpectrumImageFilter.h"

// Image definitions
const unsigned int             Dimension = 2;
typedef unsigned short         OutputPixelType;
typedef std::complex< double > ComplexPixelType;
typedef itk::Image< ComplexPixelType, Dimension > ComplexImageType;
typedef itk::Image< OutputPixelType, Dimension >  OutputImageType;

// Input and output
typedef itk::ImageFileWriter< OutputImageType > WriterType;
typedef itk::ImageFileReader< ComplexImageType > InputReaderType;

typedef itk::VisualizeComplexSpectrumImageFilter< ComplexImageType,
OutputImageType > VisualizeSpectrumType;

int main( int argc, char * argv[] )
{
	
	if( argc != 2 )
    {
        std::cerr << "Usage: " << argv[0] << " path/to/input.png" <<
std::endl;
        std::cerr << "Example: ../data/input/tagged_image.png" << std::endl;
        return EXIT_FAILURE;
    }
	
	const char * inputFileName = argv[1];
	unsigned int FFTDimension = 0;
	unsigned int NonFFTDimension = 1;
	
	// Read in the input image
    InputReaderType::Pointer inputReader = InputReaderType::New();
    inputReader->SetFileName( inputFileName );
    inputReader->Update();
    
    ComplexImageType::Pointer input = ComplexImageType::New();
    input = inputReader->GetOutput();
    input->Update();

    // Create an empty image to paste the output in
    ComplexImageType::Pointer output = ComplexImageType::New();
	output->CopyInformation( input );
	output->Allocate();
	output->FillBuffer( std::complex<double> (0, 0) );
	output->Update();
	
	// Create a region that is ONE ROW of the image
	// You will modify the index to loop through the image
	ComplexImageType::SizeType size;
	size.Fill( 1 );
	size[FFTDimension] =
input->GetLargestPossibleRegion().GetSize()[FFTDimension];
	
	ComplexImageType::IndexType index;
	index.Fill( 0 );
	
	ComplexImageType::RegionType region;
	region.SetSize( size );
	// Note that the index is not set until the loop
	
	unsigned short i = 0;
	fftw_complex *in, *out;
	fftw_plan p;
	
	while( i <= size[NonFFTDimension] - 1) {
	
		index[NonFFTDimension] = i;
		region.SetIndex( index );
		
		input->SetRequestedRegion( region );
		input->SetBufferedRegion( input->GetRequestedRegion() );
		
		output->SetRequestedRegion( region );
		output->SetBufferedRegion( input->GetRequestedRegion() );
		
		in = reinterpret_cast<fftw_complex *>(input->GetBufferPointer());
		out = reinterpret_cast<fftw_complex *>(output->GetBufferPointer());
		
		std::cout << size[FFTDimension] << std::endl; // 160
		std::cout << in << std::endl; // 0x103aef000
		std::cout << out << std::endl; // 0x7f8af48f2d20
		std::cout << FFTW_FORWARD << std::endl; // -1
		std::cout << FFTW_ESTIMATE << std::endl; // 64
		
		p = fftw_plan_dft_1d(size[FFTDimension], in, out, FFTW_FORWARD,
FFTW_ESTIMATE); // Segmentation Fault: 11
		// When I run this with gdb, I get:
		// Program received signal EXC_BAD_ACCESS, Could not access memory.
		// Reason: 13 at address: 0x0000000000000000
		// 0x00000001009f4f8e in n1_8 ()
        
        fftw_execute(p); // Doesn't run because of the segfault
		
		i = i + 1;
	
	}
    
    //////////////
    // Write it //
    //////////////
    
    VisualizeSpectrumType::Pointer vis = VisualizeSpectrumType::New();
    vis->SetInput( output );
    vis->Update();
    
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput( vis->GetOutput() );
    writer->SetFileName( "../data/output/out.png" );
    writer->Update();
 
  return EXIT_SUCCESS;
  
}



--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/1D-Complex-To-Complex-FFT-of-a-2D-Image-tp7585494.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.