'm facing with a code I'm implementing for an exam using the GPU. Specifically, the code I'm writing is in C++, and I'm using the CUFFT library to perform the Fast Fourier Transform (FFT). The purpose of the project is to resize images that have both even and odd dimensions. The current steps I am taking are as follows:
- Perform FFT with CUFFT.
- Shift low frequencies to the center of the spectrum.
- Crop the center of the Fourier spectrum using a square mask with dimensions equal to cropping dimensions of the image.
- Perform inverse FFT shift.
- Perform IFFT again with CUFFT.
The problem is: when processing square images with odd dimensions, even without using the GPU kernel for image resizing (executing only steps 1-2-4-5), the resulting image appears with colors that do not correspond to reality. I am attaching the code in which I perform the FFT shift on square images with odd dimensions and the obtained result.
Could someone help me identify the reason behind this issue?
#include <iostream>
#include <opencv2/opencv.hpp>
#include <cuda_runtime.h>
#include <cuda_runtime_api.h>
#include <cuda.h>
#include <stdlib.h>
#include <stdio.h>
#include <cufft.h>
#include <string.h>
#include <fstream>
using namespace std;
using namespace cv;
//bool forward is true in fftshift in forward sense, otherwise is false.
__global__ void fftshift(cufftDoubleComplex *im, int height, int width, bool forward)
{
int x, y, indshift;
cufftDoubleComplex c;
int idy = blockIdx.x*blockDim.x + threadIdx.x; //columns
int idx = blockIdx.y*blockDim.y + threadIdx.y; //rows
int ind = idx + idy*width;
// input sizes
int sx = width;
int sy = height;
// size of top-left quadrant
int cx = forward ? (sx + 1) / 2 : sx / 2;
int cy = forward ? (sy + 1) / 2 : sy / 2;
if(idx < cy && idy < width)
{
if(idx<cy && idy<cx)
{
x=idx+cy;
y=idy+ cx;
}
else if(idx<cy && idy>=cx)
{
x=idx+cy;
y=idy- cx;
}
indshift = x+y*width;
c.x = im[ind].x;
c.y = im[ind].y;
im[ind].x = im[indshift].x;
im[ind].y = im[indshift].y;
im[indshift].x = c.x;
im[indshift].y = c.y;
}
}
int main() {
int in_max,in_min;
cufftHandle planZ2Z, planIZ2Z;
ofstream fp;
// Read input image using OpenCV
cv::Mat inputImage = cv::imread("donnaD.jpg", 0);
if (inputImage.empty()) {
std::cerr << "Error: Couldn't read the input image." << std::endl;
return -1;
}
// Get input image dimensions
int inputWidth = inputImage.cols;
int inputHeight = inputImage.rows;
cufftDoubleComplex* h_input = (cufftDoubleComplex * )malloc(inputHeight*inputWidth*sizeof(cufftDoubleComplex));
cufftDoubleComplex* h_output = (cufftDoubleComplex * )malloc(inputHeight*inputWidth*sizeof(cufftDoubleComplex));
cufftPlan2d(&planZ2Z, inputHeight, inputWidth, CUFFT_Z2Z);
cufftPlan2d(&planIZ2Z, inputHeight, inputWidth, CUFFT_Z2Z);
// Find minimum and maximum values
in_max= static_cast<int>(inputImage.at<uchar>(0, 0));
in_min = static_cast<int>(inputImage.at<uchar>(0, 0));
for (int i = 0; i < inputImage.rows; ++i) {
for (int j = 0; j < inputImage.cols; ++j) {
int pixel_value = static_cast<int>(inputImage.at<uchar>(i, j));
if (pixel_value < in_min) {
in_min = pixel_value;
}
if (pixel_value > in_max) {
in_max = pixel_value;
}
}
}
printf("height : %d, width %d", inputHeight,inputWidth);
printf("\nMax value: %d\n",in_max);
printf("\nMin value: %d\n",in_min);
fp.open("input_donnaD.txt",ios::out);
for (int i = 0; i < inputImage.rows; ++i) {
for (int j = 0; j < inputImage.cols; ++j) {
fp<< (int)inputImage.at<uchar>(i,j)<<" ";
}
fp<<"\n";
}
fp.close();
for(int i=0;i<inputHeight;i++)
{
for(int j=0;j<inputWidth;j++)
{
h_input[i * inputWidth + j].x = (double) inputImage.at<uchar>(i,j);
h_input[i * inputWidth + j].y = 0.0f;
}
}
// Allocate memory for input and padded output on the device
cufftDoubleComplex* d_input, *d_output, *FI;
cudaMalloc((void**)&d_input, inputWidth * inputHeight * sizeof(cufftDoubleComplex));
cudaMalloc((void**)&d_output, inputWidth * inputHeight * sizeof(cufftDoubleComplex));
cudaMemset(d_output,0,inputWidth*inputHeight*sizeof(cufftDoubleComplex));
cudaMalloc((void**)&FI, inputWidth * inputHeight * sizeof(cufftDoubleComplex));
cudaMemset(FI,0,inputWidth*inputHeight*sizeof(cufftDoubleComplex));
// Copy input data from host to device using cudaMemcpy2D
cudaMemcpy2D(d_input, inputWidth * sizeof(cufftDoubleComplex), h_input, inputWidth * sizeof(cufftDoubleComplex), inputWidth * sizeof(cufftDoubleComplex), inputHeight, cudaMemcpyHostToDevice);
// Define block and grid dimensions
dim3 blockDim(1, 1);
dim3 gridDim;
gridDim.y = inputHeight/blockDim.y+((inputHeight%blockDim.y)==0? 0:1);
gridDim.x = inputWidth/blockDim.x+((inputWidth%blockDim.x)==0? 0:1);
cout<<"gridDim.x: "<<gridDim.x<<" gridDim.y: "<<gridDim.y<<"\n";
cout<<"blockDim.x: "<<blockDim.x<<" blockDim.y: "<<blockDim.y<<"\n";
// Launch the kernel
cufftExecZ2Z(planZ2Z, d_input, FI, CUFFT_FORWARD);
fftshift<<<gridDim,blockDim>>>(FI,inputHeight,inputWidth, true);
fftshift<<<gridDim,blockDim>>>(FI,inputHeight,inputWidth, false);
cufftExecZ2Z(planIZ2Z, FI, d_output, CUFFT_INVERSE);
// Copy the result back to the host using cudaMemcpy2D
cudaMemcpy(h_output,d_output,inputWidth*inputHeight*sizeof(cufftDoubleComplex), cudaMemcpyDeviceToHost);
printf("\nOutput dims of image: (%d,%d) \n",inputHeight,inputWidth);
double max = h_output[0].x;
double min = h_output[0].x;
for(int i=0;i<inputHeight;i++)
{
for(int j=0;j<inputWidth;j++)
{
if(h_output[j + i*inputWidth].x>max)
max = h_output[j + i*inputWidth].x;
if(h_output[j+ i*inputWidth].x<min)
min = h_output[j + i*inputWidth].x;
}
}
//print double max and min value
printf("\n\nMax value: %f\n",max);
printf("\nMin value: %f\n",min);
//write the value of output_res in file using fp file descriptor
fp.open("output.txt",ios::out);
for(int i=0;i<inputHeight;i++)
{
for(int j=0;j<inputWidth;j++)
fp<< round((h_output[i*inputWidth+j].x - min) * (in_max - in_min) / (max - min) + 0)<<" ";
fp<<"\n";
}
//close fp file descriptor
fp.close();
// Convert the padded output to OpenCV Mat for visualization or further processing
cv::Mat outputImage(inputHeight, inputWidth, CV_8UC1);
for (int i = 0; i < inputHeight; ++i) {
for (int j = 0; j < inputWidth; ++j) {
outputImage.at<uchar>(i, j) = round((h_output[i*inputWidth+j].x - min) * (in_max - in_min) / (max - min) + 0);
}
}
// Find minimum and maximum values
in_max= static_cast<int>(outputImage.at<uchar>(0, 0));
in_min = static_cast<int>(outputImage.at<uchar>(0, 0));
for (int i = 0; i < outputImage.rows; ++i) {
for (int j = 0; j < outputImage.cols; ++j) {
int pixel_value = static_cast<int>(outputImage.at<uchar>(i, j));
if (pixel_value < in_min) {
in_min = pixel_value;
}
if (pixel_value > in_max) {
in_max = pixel_value;
}
}
}
printf("height : %d, width %d", inputHeight,inputWidth);
printf("\nMax value: %d\n",in_max);
printf("\nMin value: %d\n",in_min);
// Display or save the padded image using OpenCV as needed
imwrite("output.jpg",outputImage);
// Cleanup
delete[] h_input;
delete[] h_output;
cudaFree(d_input);
cudaFree(d_output);
cudaFree(FI);
cufftDestroy(planZ2Z);
cufftDestroy(planIZ2Z);
return 0;
}
Input image - size 301x301:
Output image - size 301x301:
I would like to obtain the same image as input, but using only the FFT in the forward sense, double FFT shift to adjust the zero frequencies, and inverse FFT the result is wrong.
For the time being, I set aside the initial issue of resizing images as I am already encountering problems with just these steps.