Performing convolution in frequency domain manually but getting wrong output image in CPP/Opencv

Question

I followed the following steps:- 1. Calculated dft of image 2. Calculated dft of kernel (but 1st padded it to size of image) 3. Multiplied real and imaginary parts of both dft individually 4. Calculated inverse dft I tried to display the images in each intermediate step but the final image comes out to be almost black except in corners. Image fourier transform output after multiplication and its inverse dft output

input image

enter code here

#include <iostream>
#include <stdlib.h>
#include <opencv2/opencv.hpp>
#include <stdio.h>
int r=100;
#define SIGMA_CLIP 6.0f
using namespace cv;
using namespace std;

void updateResult(Mat complex)
{
Mat work;
idft(complex, work);
Mat planes[] = {Mat::zeros(complex.size(), CV_32F), Mat::zeros(complex.size(), CV_32F)};
split(work, planes);                // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))

magnitude(planes[0], planes[1], work);    // === sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)
normalize(work, work, 0, 1, NORM_MINMAX);
imshow("result", work);
}
void shift(Mat magI) {

// crop if it has an odd number of rows or columns
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

int cx = magI.cols/2;
int cy = magI.rows/2;

Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

Mat tmp;                            // swap quadrants (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);                     // swap quadrant (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2);
}
Mat updateMag(Mat complex )
{

Mat magI;
Mat planes[] = {Mat::zeros(complex.size(), CV_32F), Mat::zeros(complex.size(), CV_32F)};
split(complex, planes);                // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))

magnitude(planes[0], planes[1], magI);    // sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)

// switch to logarithmic scale: log(1 + magnitude)
magI += Scalar::all(1);
log(magI, magI);

shift(magI);
normalize(magI, magI, 1, 0, NORM_INF); // Transform the matrix with float values into a
         return magI;                                 // viewable image form (float between values 0 and 1).
//imshow("spectrum", magI);
 }
 Mat createGausFilterMask(Size imsize, int radius) {

// call openCV gaussian kernel generator
double sigma = (r/SIGMA_CLIP+0.5f);
Mat kernelX = getGaussianKernel(2*radius+1, sigma, CV_32F);
Mat kernelY = getGaussianKernel(2*radius+1, sigma, CV_32F);
// create 2d gaus
Mat kernel = kernelX * kernelY.t();

int w = imsize.width-kernel.cols;
int h = imsize.height-kernel.rows;

int r = w/2;
int l = imsize.width-kernel.cols -r;

int b = h/2;
int t = imsize.height-kernel.rows -b;

Mat ret;
copyMakeBorder(kernel,ret,t,b,l,r,BORDER_CONSTANT,Scalar::all(0));

return ret;

}

//code reference https://docs.opencv.org/2.4/doc/tutorials/core/discrete_fourier_transform/discrete_fourier_transform.html
int main( int argc, char** argv )
{ 

String file;
file = "lena.png";

Mat image = imread(file, CV_LOAD_IMAGE_GRAYSCALE);

Mat padded;                             

int m = getOptimalDFTSize( image.rows );
int n = getOptimalDFTSize( image.cols );  
copyMakeBorder(image, padded, 0, m - image.rows, 0, n -image.cols, BORDER_CONSTANT, Scalar::all(0));//expand input image to optimal size , on the border add zero values

Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI);  
dft(complexI, complexI); //computing dft
split(complexI, planes); //image converted to complex and real dft here

Mat mask = createGausFilterMask(padded.size(),r );  // Forming the gaussian filter
Mat mplane[] = {Mat_<float>(mask), Mat::zeros(mask.size(), CV_32F)};
Mat kernelcomplex;
merge(mplane, 2, kernelcomplex); 

dft(kernelcomplex, kernelcomplex);

split(kernelcomplex, mplane);// splitting the dft of kernel to real and complex 
mplane[1]=mplane[0]; //overwriting imaginary values with real values of kernel dft
Mat kernel_spec;
merge(mplane, 2, kernel_spec);
mulSpectrums(complexI, kernel_spec, complexI, DFT_ROWS);
Mat magI=updateMag(complexI);

namedWindow( "image fourier", CV_WINDOW_AUTOSIZE );

imshow("spectrum magnitude", magI);


updateResult(complexI); //converting to viewable form, computing idft

waitKey(0); 

return 0;
}

Which step is going wrong? Or am i missing on to some concept?

---

Edited the code with help of Cris and it now works perfectly.

Answer 1

There are two immediately apparent issues:

1. The Gaussian is real-valued and symmetric. Its Fourier transform should be too. If the DFT of your kernel has a non-zero imaginary component, you're doing something wrong.

   Likely, what you are doing wrong is that your kernel has its origin in the middle of the image, rather than at the top-left sample. This is the same issue as in this other question. The solution is to use the equivalent of MATLAB's ifftshift, an implementation of which is shown in the OpenCV documentation ("step 6, Crop and rearrange").

2. To apply the convolution, you need to multiply the two DFTs together, not the real parts and imaginary parts of the DFTs. Multiplying two complex numbers a+ib and c+id results in ac-bd+iad+ibc, not ac+ibd.

   But since the DFT of your kernel should be real-valued only, you can simply multiply the real component of the kernel with both the real and imaginary components of the image: (a+ib)c = ac+ibc.

---

It seems very roundabout what you are doing with the complex-valued images. Why not let OpenCV handle all of that for you? You can probably* just do something like this:

Mat image = imread(file, CV_LOAD_IMAGE_GRAYSCALE);

// Expand input image to optimal size, on the border add zero values
Mat padded;                             
int m = getOptimalDFTSize(image.rows);
int n = getOptimalDFTSize(image.cols);  
copyMakeBorder(image, padded, 0, m - image.rows, 0, n -image.cols, BORDER_CONSTANT, Scalar::all(0));

// Computing DFT
Mat DFTimage;
dft(padded, DFTimage); 

// Forming the Gaussian filter
Mat kernel = createGausFilterMask(padded.size(), r);
shift(kernel);
Mat DFTkernel;
dft(kernel, DFTkernel);

// Convolution
mulSpectrums(DFTimage, DFTkernel, DFTimage, DFT_ROWS);

// Display Fourier-domain result
Mat magI = updateMag(DFTimage);
imshow("spectrum magnitude", magI);

// IDFT
Mat work;
idft(complex, work); // <- NOTE! Don't inverse transform log-transformed magnitude image!

Note that the Fourier-Domain result is actually a special representation of the complex-conjugate symmetric DFT, intended to save space and computations. To compute the full complex output, add the DFT_COMPLEX_OUTPUT to the call to dft, and DFT_REAL_OUTPUT to the call to idft (this latter then assumes symmetry, and produces a real-valued output, saving you the hassle of computing the magnitude).

^{* I say probably because I haven't compiled any of this... If there's something wrong, please let me know, or edit the answer and fix it.}