DISCLAIMER UP FRONT: I know I can smooth/blur images using OpenCV, but my assignment is more under the hood so I can understand the process. I'm new to image processing and DFT's in general, so this is proving to be a challenge.
DESCRIPTION:
I'm doing some image processing using the FFTW library and OpenCV to read an image into an OpenCV::Mat object and grayscale color space. I convert the data to type double so I can create a double pointer to it which FFTW's function for my FFT requires. This is in a Mat object named imageDouble and the pointer to it is pImageDouble.
std::string filename = "<path_to_your_image>";
Mat image = imread(filename, IMREAD_GRAYSCALE);
/***********************************************************************
*** creating Mat object of type *doube* and copying image contents to it
*** then creating a pointer to it.
***********************************************************************/
Mat imageDouble;
image.convertTo(imageDouble, CV_64FC1);
double* pImageDouble = imageDouble.ptr<double>(0);
I also have a Gaussian kernel that I've read into a Mat object. I perform a circular shift where the center (and highest) value of the Gaussian kernel is shifted to the top left (the (0,0) position) of the kernel, and then I zero-pad the kernel to the size of the original image that I read in. The results (type double) are stored in a Mat object named paddedGKernel and the pointer is named pPaddedGKernel.
Mat paddedGKernel = padGKernelWithZeros(nRows, nCols, rolledGKernel);
double* pPaddedGKernel = paddedGKernel.ptr<double>(0);
I initializefftw_complex objects to output the FFT results into, an fftw_plan, and then I allocate the memory for the fftw_complex objects and execute the plan. I use FFTW's fftw_plan_dft_r2c_2d() function to perform an FFT on both 2d objects (image and shifted/padded Gaussian kernel), then perform a pointwise multiplication in Fourier space to apply the Gaussian filter to the original image.
fftw_complex* out; //for result of FFT for original image
fftw_complex* outGaussian; //for result of FFT for Gaussian kernel
fftw_plan p;
/*************************************************
*** Allocating memory for the fftw_complex objects
*************************************************/
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * nCols);
outGaussian = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * nCols);
/****************************************
*** FFT on imageDouble, outputting to out
****************************************/
p = fftw_plan_dft_r2c_2d(imageDouble.rows, imageDouble.cols, pImageDouble, out, FFTW_ESTIMATE);
fftw_execute(p);
/**************************************************
*** FFT on paddedGKernel, outputting to outGuassian
**************************************************/
p = fftw_plan_dft_r2c_2d(paddedGKernel.rows, paddedGKernel.cols, pPaddedGKernel, outGaussian, FFTW_ESTIMATE);
fftw_execute(p);
/****************************************************
*** Pointwise multiplication to apply Gaussian kernel
****************************************************/
for (size_t i = 0; i < nCCols * nRows; i++)
{
out[i][0] = out[i][0] * (1 / outGaussian[i][0]);
}
I then perform the inverse FFT fftw_plan_dft_c2r_2d() on the result of the pointwise multiplication and output that to an OpenCV::Mat object imageCloneDouble and convert the data back to uchar, placing the data in OpenCV::Mat imageBack.
/***********************************************************************
*** Making Mat object (type *double*) to put data into and pointer to it
***********************************************************************/
Mat imageCloneDouble = Mat::zeros(nRows, nCols, CV_64FC1);
double* pImageCloneDouble = imageCloneDouble.ptr<double>(0);
/*****************************************************************************
*** Making and executing plan for inverse FFT, outputting to imageCloneDouble,
*** then normalizing since this function puts out unnormalized values.
*****************************************************************************/
fftw_plan pp = fftw_plan_dft_c2r_2d(imageCloneDouble.rows, imageCloneDouble.cols, out, pImageCloneDouble, FFTW_BACKWARD | FFTW_ESTIMATE);
fftw_execute(pp);
imageCloneDouble = imageCloneDouble / (nCols * nRows *2);
/***************************************************************************
*** Making Mat object (type *uchar*) and copying data inverse FFT data to it
***************************************************************************/
Mat imageBack;
imageCloneDouble.convertTo(imageBack, CV_8UC1);
When I show the transformed image imageBack I expect to get the original with the Gaussian kernel applied, but it looks like the original with potentially some modification overlaid with a rotated version of the original image that appears to have some modification. I can't understand where and why this flip/rotation is happening and why it is overlaid with what appears to be the original image, but I suspect it happens when I'm in Fourier space and performing a pointwise or element-wise multiplication. Either that or I'm omitting a process I need to perform on the result of the pointwise multiplication.
What do I need to do to this data in Fourier space to get back to the original image?
The key to fixing this issue started with what Cris pointed out as the second problem. So the first thing to do was to implement the code provided in that part of the reply.
auto* out_p = reinterpret_cast<std::complex<double>*>(out);
auto* outGaussian_p = reinterpret_cast<std::complex<double>*>(outGaussian);
for (std::size_t i = 0; i < N; i++)
{
out_p[i] *= outGaussian_p[i];
}
I had the Gaussian kernel normalized the usual way (dividing by the sum of all the values). After making the changes above from Cris, the image had the blur applied and the only issue was the range of pixel values had changed.
To fix that, I applied the formula for contrast stretching to match the original range and that solved the problem.
/* Finding the original image's minimum and maximum value */
double oMin, oMax;
minMaxLoc(image, &oMin, &oMax);
/* Finding min and max values for our 'imageBack' and rounding them. */
double min, max;
minMaxLoc(imageCloneDouble, &min, &max);
min = round(min);
max = round(max);
/* Code version of formula for contrast stretching */
if (oMin != min | oMax != max) {
double numerator = oMax - oMin;
double denominator = max - min;
double scaledMin = -(min)*numerator;
double scaledOMin = oMin * denominator;
double innerTerm = scaledMin + scaledOMin;
imageCloneDouble = ((numerator * imageCloneDouble) + innerTerm ) / denominator;
}