Sir,
- 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…