deblurring image by deconvolution using opencv

Question

I have two images o1 & o2, and I have blurred the two images using the same Gaussian blurring kernel. Then I have found kernel k1 = DFT(b1) / DFT (o1), where b1 is the image obtained by blurring o1.

I have used this kernal (k1) to perform deconvolution on b2, where b2 is obtained by blurring o2.

But deblurred output is not correct (the output image does not have any relation with original) what is the problem in my code ?

int main(int argc, char** argv) 
{
  Mat orig1 = imread(argv[1], 0);
  Mat orig2 = imread(argv[2], 0);

  Mat blur1, blur2;
  GaussianBlur(orig1, blur1, Size(11, 11), 0, 0 );
  GaussianBlur(orig2, blur2, Size(11, 11), 0, 0 );

  imshow("or1", orig1);
  imshow("bl1", blur1);
  imshow("or2", orig2);
  imshow("bl2", blur2);
  waitKey(0);



  deconvolution(orig1, blur1, orig2, blur2);

  return 0;
}
void deconvolution(Mat & o1, Mat & b1, Mat & o2, Mat & b2)
{
  Mat o1f, o2f, b1f, b2f;
  Mat o1dft, o2dft, b1dft, b2dft;

  o1.convertTo(o1f, CV_32F);
  b1.convertTo(b1f, CV_32F);
  o2.convertTo(o2f, CV_32F);
  b2.convertTo(b2f, CV_32F);

  computeDFT(o1f, o1dft);
  computeDFT(b1f, b1dft);
  computeDFT(o2f, o2dft);
  computeDFT(b2f, b2dft);

  Mat k1, k2, b1d, b2d;
  divide(b1dft, o1dft, k1);

  Mat r1, r2;
  divide(b1dft, k1, r1);
  divide(b2dft, k1, r2);

  Mat idftr1, idftr2;
  computeIDFT(r1, idftr1);
  computeIDFT(r2, idftr2);

  Mat r1_8u, r2_8u;
  idftr1.convertTo(r1_8u, CV_8U);
  idftr2.convertTo(r2_8u, CV_8U);

  imshow("r1", r1_8u);
  imshow("r2", r2_8u);
  waitKey(0);
  destroyAllWindows();
}

Images o1, o2, b1, b2, r1 and r2 are given in order below:

Answer 1

The problem is most likely that your blurring kernel has vanishing coefficients for certain frequencies. For each coefficient of the transform of your signal (f) and blurring kernel (h), you calculate f/h right now. This is effectively a division by zero for these coefficients, resulting in the strong noise you observe.

A quick solution for this would be pseudo-inverse filtering:

use f/h only for |h| > epsilon

set coefficient to 0 else

If this isn't smooth enough, you can get better results with wiener filtering.