Unconventional FFT output image observed when using cufft............???

Accelerated Computing CUDA CUDA Programming and Performance
1

Sir,

  1. I AM trying to find the FFT of a DICOM IMAGE using cuFFT and then display it.

THE STEPS DONE ARE:

-LOADING THE DICOM IMAGE USING VTK

-CONVERTING THE VTK DICOM IMAGE TO C ARRAY

-DO CUFFT ON THE C ARRAY

-GET BACK THE RESULT OF CUFFT AND CONVERT TO VTK IMAGE FOR VISUALIZING THE FFT

-FINALLY I WANT TO TRY THE ABOVE SAME USING VOLUME INPUT.

ADDITIONAL INFORMATIONS:

-SIZE OF IMAGE IS 512512(FOR 2D) ,512512*512 (FOR 3D)
-PIXEL SPACING IS .5mm
-loading and displaying operations are done in main.cpp.The fft magnitude calculation function is done in a .cu file.Then i am calling this function from main.cpp.

THE ISSUES:
FFT OUTPUT IMAGE IS NOT AS REQUIRED.FROM MY UNDERSTANDING A NON CENTERED FFT IMAGE IS DIAGONALLY

SYMMETRICAL AND VALUES ARE MORE AT THE CORNERS(AS THE IMAGE IS A LOWPASS SIGNAL).BUT I AM GETTING WRONG FFT DISPLAY.
ALSO WHEN I PRINTED THE FFT OUTPUT VALUES I AM GETTING 1.#QNAN VALUES

THE CODE USED IS GIVEN BELOW:

fftmagnitude.cu

global void magnitudeFFT(Complex* idata, float* odata)
{
const int i = blockIdx.x * blockDim.x + threadIdx.x;
const int j = blockIdx.y * blockDim.y + threadIdx.y;

odata[j*SIGNAL_SIZE + i] = cuCabsf(idata[j*SIGNAL_SIZE + i]) / (SIGNAL_SIZE*SIGNAL_SIZE);

}

void fftMagnitude(Complex* h_signal, float* h_out){

// Allocate device memory for the image
Complex *d_signal;
cudaMalloc((void **)&d_signal, SIGNAL_SIZE*SIGNAL_SIZE*sizeof(Complex));

// Copy host memory to device
cudaMemcpy(d_signal, h_signal, (SIGNAL_SIZE*SIGNAL_SIZE*sizeof(Complex)), cudaMemcpyHostToDevice);

// CUFFT 
cufftHandle plan;

cufftPlan2d(&plan, SIGNAL_SIZE, SIGNAL_SIZE, CUFFT_C2C); 		
cufftExecC2C(plan, (cufftComplex*)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD);

// Allocate device memory for signalout
float *d_signal_out;
cudaMalloc((void **)&d_signal_out,SIGNAL_SIZE*SIGNAL_SIZE*sizeof(float));

dim3 blockDim(16, 8, 1);
dim3 gridDim(SIGNAL_SIZE / 16, SIGNAL_SIZE / 8, 1);

//launching the kernal..
magnitudeFFT << <gridDim, blockDim >> >(d_signal, d_signal_out);

//Destroy CUFFT context
cufftDestroy(plan);

//cleanup memory
//free(h_signal);
cudaFree(d_signal);

cudaMemcpy(h_out, d_signal_out, SIGNAL_SIZE*SIGNAL_SIZE*sizeof(float), cudaMemcpyDeviceToHost);

}

main.cpp

int main()
{

//................reading image....................
vtkSmartPointer<vtkDICOMImageReader> fixedImageReader1 = vtkSmartPointer<vtkDICOMImageReader>::New();
fixedImageReader1->SetFileName("D:\\headOneSlices512\\Slice255.dcm");
fixedImageReader1->Update();



float * volume1 = (float*)malloc(sizeof(float)*SIGNAL_SIZE*SIGNAL_SIZE);

vtkSmartPointer<vtkImageExport> exporter = vtkSmartPointer<vtkImageExport>::New();
exporter->SetInputData(fixedImageReader1->GetOutput());
//exporter->ImageLowerLeftOn();//default
exporter->ImageLowerLeftOff();
exporter->Export(volume1);
exporter->Update();


//..............................body of the programme.................................//


// Allocate host memory for the signal
float* h_signal_real1 = (float*)malloc(sizeof(float)* SIGNAL_SIZE*SIGNAL_SIZE);
Complex* h_signal1 = (Complex*)malloc(sizeof(Complex)* SIGNAL_SIZE*SIGNAL_SIZE);

//making complex input
	for (int j = 0; j < SIGNAL_SIZE; ++j) {
		for (int i = 0; i < SIGNAL_SIZE; ++i) {

			h_signal1[j*SIGNAL_SIZE + i].x = volume1[j*SIGNAL_SIZE + i];
			h_signal1[j*SIGNAL_SIZE + i].y = 0;

			h_signal_real1[j*SIGNAL_SIZE + i] = volume1[j*SIGNAL_SIZE + i];
		}
	}


//Allocate host memory for the  signal
float* h_fftoutput_signal = (float*)malloc(sizeof(float)*SIGNAL_SIZE*SIGNAL_SIZE);


//calling fftmagnitude function from .cu file
fftMagnitude(h_signal1, h_fftoutput_signal);

//displaying fft values
for (int j = 0; j < SIGNAL_SIZE /32; ++j) {
	for (int i = 0; i < SIGNAL_SIZE /32; ++i) {	
		std::cout << h_fftoutput_signal[j*SIGNAL_SIZE + i] << "  ";
	}
	std::cout << std::endl;
}





//..................................image view part......................................// 



// Convert the c-style image to a vtkImageData


vtkSmartPointer<vtkImageImport> imageImport =
	vtkSmartPointer<vtkImageImport>::New();
imageImport->SetDataSpacing(.5, .5, 0);
imageImport->SetDataOrigin(0, 0, 0);
imageImport->SetWholeExtent(0, SIGNAL_SIZE-1, 0, SIGNAL_SIZE-1, 0, 0);
imageImport->SetDataExtentToWholeExtent();
imageImport->SetNumberOfScalarComponents(1);
imageImport->SetImportVoidPointer(h_fftoutput_signal);
imageImport->Update();



//interactor and viewer............................
vtkSmartPointer<vtkResliceImageViewer>viewer1 = vtkSmartPointer<vtkResliceImageViewer>::New();
vtkSmartPointer<vtkResliceImageViewer>viewer2 = vtkSmartPointer<vtkResliceImageViewer>::New();

vtkSmartPointer<vtkRenderWindowInteractor>interactor1 = vtkSmartPointer<vtkRenderWindowInteractor>::New();
vtkSmartPointer<vtkRenderWindowInteractor>interactor2 = vtkSmartPointer<vtkRenderWindowInteractor>::New();


interactor1->SetRenderWindow(viewer1->GetRenderWindow());
interactor2->SetRenderWindow(viewer2->GetRenderWindow());

viewer1->SetupInteractor(interactor1);
viewer2->SetupInteractor(interactor2);


viewer1->SetInputData(fixedImageReader1->GetOutput());
viewer2->SetInputData(imageImport->GetOutput());

viewer1->GetRenderer()->ResetCamera();
viewer1->Render();
viewer2->GetRenderer()->ResetCamera();
viewer2->Render();

interactor1->Start();
interactor2->Start();

return 0;

}

SOMEONE PLEASE HELP ME…