How to deblur image using fourier transform in open-cv or emgu-cv?

Question

i saw this video about debluring images using fourier transform in matlab video

and i want to convert the code in emgu cv

my code in emgucv :

  string path = Environment.GetFolderPath(Environment.SpecialFolder.Desktop);
        Image<Bgr, byte> img = new Image<Bgr, byte>(@"lal.png");


        //blur the image
        Image<Gray, byte> gray = img.Convert<Gray, byte>().SmoothBlur(31,31);

        //convert image to float and get the fourier transform
        Mat g_fl = gray.Convert<Gray, float>().Mat;
        Matrix<float> dft_image = new Matrix<float>(g_fl.Size);
        CvInvoke.Dft(g_fl, dft_image, Emgu.CV.CvEnum.DxtType.Forward, 0);

        //here i make an image of kernel with size of the original 


        Image<Gray, float> ker = new Image<Gray, float>(img.Size);
        ker.SetZero();

        for (int x = 0; x < 31; x++)
        {
            for (int y = 0; y < 31; y++)
            {

                //31 * 31= 961
                    ker[y, x] = new Gray(1/961);


            }
        }

        //get the fourier of the kernel
        Matrix<float> dft_blur = new Matrix<float>(g_fl.Size);
        CvInvoke.Dft(ker, dft_blur, Emgu.CV.CvEnum.DxtType.Forward, 0);


        // fouier image / fourier blur
        Matrix<float> res = new Matrix<float>(g_fl.Size);
        for (int x=0;x<g_fl.Cols;x++)
        {
            for (int y = 0; y < g_fl.Rows; y++)
            {

                    res[y, x] = dft_image[y, x] / dft_blur[y, x];
            }
        }

        //get the inverse of fourier
        Image<Gray, float> ready = new Image<Gray, float>(g_fl.Size);

        CvInvoke.Dft(res, ready, Emgu.CV.CvEnum.DxtType.Inverse, 0);

        CvInvoke.Imshow("deblur", ready.Convert<Gray,byte>());


        CvInvoke.Imshow("original", gray);



        CvInvoke.WaitKey(0);

but the result is black and not working , where is the mistake in my code if you have a code in opencv python you can post it :)??

Thanks :)

Answer 1

My old implementation of wiener filter:

#include "stdafx.h"
#pragma once
#pragma comment(lib, "opencv_legacy220.lib")
#pragma comment(lib, "opencv_core220.lib")
#pragma comment(lib, "opencv_highgui220.lib")
#pragma comment(lib, "opencv_imgproc220.lib")
#include "c:\Users\Andrey\Documents\opencv\include\opencv\cv.h"
#include "c:\Users\Andrey\Documents\opencv\include\opencv\cxcore.h"
#include "c:\Users\Andrey\Documents\opencv\include\opencv\highgui.h"
#include <string>
#include <iostream>
#include <complex>

using namespace std;
using namespace cv;

//----------------------------------------------------------
// Compute real and implicit parts of FFT for given image
//----------------------------------------------------------
void ForwardFFT(Mat &Src, Mat *FImg)
{
    int M = getOptimalDFTSize( Src.rows );
    int N = getOptimalDFTSize( Src.cols );
    Mat padded;    
    copyMakeBorder(Src, padded, 0, M - Src.rows, 0, N - Src.cols, BORDER_CONSTANT, Scalar::all(0));
    // Create complex representation of our image
    // planes[0] Real part, planes[1] Implicit part (zeros)
    Mat planes[] = {Mat_<double>(padded), Mat::zeros(padded.size(), CV_64F)};
    Mat complexImg;
    merge(planes, 2, complexImg); 
    dft(complexImg, complexImg);    
    // As result, we also have Re and Im planes
    split(complexImg, planes);
    // Crop specter, if it have odd number of rows or cols
    planes[0] = planes[0](Rect(0, 0, planes[0].cols & -2, planes[0].rows & -2));
    planes[1] = planes[1](Rect(0, 0, planes[1].cols & -2, planes[1].rows & -2));
    FImg[0]=planes[0].clone();
    FImg[1]=planes[1].clone();
}
//----------------------------------------------------------
// Restore our image using specter
//----------------------------------------------------------
void InverseFFT(Mat *FImg,Mat &Dst)
{
    Mat complexImg;
    merge(FImg, 2, complexImg);
    // Apply inverse FFT
    idft(complexImg, complexImg);
    split(complexImg, FImg);
    Dst=FImg[0];
}
//----------------------------------------------------------
// Wiener filter 
//----------------------------------------------------------
void wienerFilter(Mat &src,Mat &dst,Mat &_h,double k)
{
//---------------------------------------------------
// small number for numeric stability
//---------------------------------------------------
    const double eps=1E-8;
//---------------------------------------------------
    int ImgW=src.size().width;
    int ImgH=src.size().height;
//--------------------------------------------------
    Mat Yf[2];
    ForwardFFT(src,Yf); 
//--------------------------------------------------
    Mat h;
    h.create(ImgH,ImgW,CV_64F);
    h=0;    
    _h.copyTo(h(Rect(0, 0, _h.size().width, _h.size().height)));
    Mat Hf[2];
    ForwardFFT(h,Hf);
//--------------------------------------------------
    Mat Fu[2];
    Fu[0].create(ImgH,ImgW,CV_64F);
    Fu[1].create(ImgH,ImgW,CV_64F);

    complex<double> a;
    complex<double> b;
    complex<double> c;

    double Hf_Re;
    double Hf_Im;
    double Phf;
    double hfz;
    double hz;
    double A;

    for (int i=0;i<Hf[0].size().height;i++)
    {
        for (int j=0;j<Hf[0].size().width;j++)
        {
            Hf_Re=Hf[0].at<double>(i,j);
            Hf_Im=Hf[1].at<double>(i,j);
            Phf = Hf_Re*Hf_Re+Hf_Im*Hf_Im;
            hfz=(Phf<eps)*eps;
            hz =(h.at<double>(i,j)>0);
            A=Phf/(Phf+hz+k);
            a=complex<double>(Yf[0].at<double>(i,j),Yf[1].at<double>(i,j));
            b=complex<double>(Hf_Re+hfz,Hf_Im+hfz);
            c=a/b; // Deconvolution
                   // Other we do to remove division by 0
            Fu[0].at<double>(i,j)=(c.real()*A);
            Fu[1].at<double>(i,j)=(c.imag()*A); 
        }
    }
//--------------------------------------------------    
    Fu[0]/=(ImgW*ImgH);
    Fu[1]/=(ImgW*ImgH);
//--------------------------------------------------
    InverseFFT(Fu,dst);
    // remove out of rane values
    for (int i=0;i<Hf[0].size().height;i++)
    {
        for (int j=0;j<Hf[0].size().width;j++)
        {
            if(dst.at<double>(i,j)>215){dst.at<double>(i,j)=215;}
            if(dst.at<double>(i,j)<(-40)){dst.at<double>(i,j)=(-40);}
        }
    }
}
//----------------------------------------------------------
// Main
//----------------------------------------------------------
int _tmain(int argc, _TCHAR* argv[])
{
// Input image
    Mat img;
// Load it from drive
    img=imread("data/motion_fuzzy_lena.bmp",0);
//---------------------------------------------
    imshow("Src image", img);
// Image size
    int ImgW=img.size().width;
    int ImgH=img.size().height;
// Deconvolution kernel (coefficient sum must be 1)
// Image was blurred using same kernel
    Mat h;
    h.create(1,10,CV_64F);
    h=1/double(h.size().width*h.size().height);
// Apply filter
    wienerFilter(img,img,h,0.05);
    normalize(img,img, 0, 1, CV_MINMAX);
    imshow("Result image", img);
    cvWaitKey(0);
    return 0;
}

The result: