使用OpenCV进行Weiner反卷积 [英] Weiner deconvolution using opencv

问题描述

我已经开发出一种方法来估计运动模糊的点扩散函数，但是我想使用PSF进行反卷积.我决定使用wiener方法.

I have developed a way to estimate the point spread function of a motion blur, however I'd like to use the PSF to perform deconvolution. I decided to use the wiener method.

  cv::Mat deconvolution(cv::Mat input, cv::Mat kernel){
     cv::Mat Fin, Fkern, padded_kern, Fdeblur,out;
     cv::normalize(kernel,kernel);
     cv::dft(input,Fin,cv::DFT_COMPLEX_OUTPUT);
     cv::copyMakeBorder(kernel,padded_kern,0,Fin.rows-kernel.rows,0,Fin.cols-kernel.cols,cv::BORDER_CONSTANT,cv::Scalar::all(0));
     cv::dft(padded_kern,Fkern,cv::DFT_COMPLEX_OUTPUT);
     cv::mulSpectrums(Fin,Fkern,Fdeblur,0,true);
     cv::dft(Fdeblur,out,cv::DFT_INVERSE|cv::DFT_REAL_OUTPUT);
     cv::normalize(out,out,0, 1, CV_MINMAX);
     return out;
   }

但是，即使将cv::mulSpectrums(Fin,Fkern,Fdeblur,0,true); last选项设置为true之后，我似乎仍在执行正常的卷积.最后一个true选项不应该意味着我要乘以共轭并因此除以内核吗?

However even after setting cv::mulSpectrums(Fin,Fkern,Fdeblur,0,true); last option to true, I still seem to be performing a normal convolution. Should'nt the last true option mean that I am multiplying by the conjugate and therefore dividing the kernel?

推荐答案

只需实现此过滤器即可:

Just implemented this filter:

#include <windows.h>
#include <iostream>
#include <vector>
#include <stdio.h>
#include "fstream"
#include "iostream"
#include <algorithm>
#include <iterator>
#include "opencv2/opencv.hpp"

using namespace std;
using namespace cv;

#include <iostream>

#include <vector>

void Recomb(Mat &src, Mat &dst)
{
    int cx = src.cols >> 1;
    int cy = src.rows >> 1;
    Mat tmp;
    tmp.create(src.size(), src.type());
    src(Rect(0, 0, cx, cy)).copyTo(tmp(Rect(cx, cy, cx, cy)));
    src(Rect(cx, cy, cx, cy)).copyTo(tmp(Rect(0, 0, cx, cy)));
    src(Rect(cx, 0, cx, cy)).copyTo(tmp(Rect(0, cy, cx, cy)));
    src(Rect(0, cy, cx, cy)).copyTo(tmp(Rect(cx, 0, cx, cy)));
    dst = tmp;
}

void convolveDFT(Mat& A, Mat& B, Mat& C)
{
    // reallocate the output array if needed
    C.create(abs(A.rows - B.rows) + 1, abs(A.cols - B.cols) + 1, A.type());
    Size dftSize;
    // compute the size of DFT transform
    dftSize.width = getOptimalDFTSize(A.cols + B.cols - 1);
    dftSize.height = getOptimalDFTSize(A.rows + B.rows - 1);
    // allocate temporary buffers and initialize them with 0's
    Mat tempA(dftSize, A.type(), Scalar::all(0));
    Mat tempB(dftSize, B.type(), Scalar::all(0));
    // copy A and B to the top-left corners of tempA and tempB, respectively
    Mat roiA(tempA, Rect(0, 0, A.cols, A.rows));
    A.copyTo(roiA);
    Mat roiB(tempB, Rect(0, 0, B.cols, B.rows));
    B.copyTo(roiB);
    // now transform the padded A & B in-place;
    // use "nonzeroRows" hint for faster processing
    dft(tempA, tempA, 0, A.rows);

    dft(tempB, tempB, 0, A.rows);

    // multiply the spectrums;
    // the function handles packed spectrum representations well
    mulSpectrums(tempA, tempB, tempA, 0);
    // transform the product back from the frequency domain.
    // Even though all the result rows will be non-zero,
    // you need only the first C.rows of them, and thus you
    // pass nonzeroRows == C.rows

    dft(tempA, tempA, DFT_INVERSE + DFT_SCALE);
    // now copy the result back to C.
    C = tempA(Rect((dftSize.width - A.cols) / 2, (dftSize.height - A.rows) / 2, A.cols, A.rows)).clone();
    // all the temporary buffers will be deallocated automatically
}

//----------------------------------------------------------
// Compute Re and Im planes of FFT from 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));
    Mat planes[] = { Mat_<double>(padded), Mat::zeros(padded.size(), CV_64FC1) };
    Mat complexImg;
    merge(planes, 2, complexImg);
    dft(complexImg, complexImg);
    split(complexImg, planes);
    // crop result
    planes[0] = planes[0](Rect(0, 0, Src.cols, Src.rows));
    planes[1] = planes[1](Rect(0, 0, Src.cols, Src.rows));
    FImg[0] = planes[0].clone();
    FImg[1] = planes[1].clone();
}
//----------------------------------------------------------
// Compute image from Re and Im parts of FFT
//----------------------------------------------------------
void InverseFFT(Mat *FImg, Mat &Dst)
{
    Mat complexImg;
    merge(FImg, 2, complexImg);
    dft(complexImg, complexImg, DFT_INVERSE + DFT_SCALE);
    split(complexImg, FImg);
    Dst = FImg[0];
}
//----------------------------------------------------------
// wiener Filter
//----------------------------------------------------------
void wienerFilter(Mat &src, Mat &dst, Mat &_h, double k)
{
    //---------------------------------------------------
    // Small epsilon to avoid division by zero
    //---------------------------------------------------
    const double eps = 1E-8;
    //---------------------------------------------------
    int ImgW = src.size().width;
    int ImgH = src.size().height;
    //--------------------------------------------------
    Mat Yf[2];
    ForwardFFT(src, Yf);
    //--------------------------------------------------
    Mat h = Mat::zeros(ImgH, ImgW, CV_64FC1);

    int padx = h.cols - _h.cols;
    int pady = h.rows - _h.rows;

    copyMakeBorder(_h, h, pady / 2, pady - pady / 2, padx / 2, padx - padx / 2, BORDER_CONSTANT, Scalar::all(0));

    Mat Hf[2];
    ForwardFFT(h, Hf);


    //--------------------------------------------------
    Mat Fu[2];
    Fu[0] = Mat::zeros(ImgH, ImgW, CV_64FC1);
    Fu[1] = Mat::zeros(ImgH, ImgW, CV_64FC1);

    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 < h.rows; i++)
    {
        for (int j = 0; j < h.cols; 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 work to avoid division by zero
            Fu[0].at<double>(i, j) = (c.real()*A);
            Fu[1].at<double>(i, j) = (c.imag()*A);
        }
    }
    InverseFFT(Fu, dst);
    Recomb(dst, dst);
}

// ---------------------------------
// 
// ---------------------------------

int main(int argc, char** argv)
{
    namedWindow("Image");
    namedWindow("Kernel");
    namedWindow("Result");
    Mat Img = imread("F:\\ImagesForTest\\lena.jpg", 0); // Source image
    Img.convertTo(Img, CV_32FC1, 1.0 / 255.0);
    Mat kernel = imread("F:\\ImagesForTest\\Point.jpg", 0); // PSF
    //resize(kernel, kernel, Size(), 0.5, 0.5);
    kernel.convertTo(kernel, CV_32FC1, 1.0 / 255.0);
    float kernel_sum = cv::sum(kernel)[0];
    kernel /= kernel_sum;

    int width = Img.cols;
    int height = Img.rows;
    Mat resim;
    convolveDFT(Img, kernel, resim);
    Mat resim2;
    kernel.convertTo(kernel, CV_64FC1);
    // Apply filter
    wienerFilter(resim, resim2, kernel, 0.01);
    imshow("Результат фильтрации", resim2);
    imshow("Kernel", kernel * 255);
    imshow("Image", Img);
    imshow("Result", resim);
    cvWaitKey(0);
}

结果如下所示(您看到的不是100％恢复):

Results look like this (as you see it not 100% restoration):

这篇关于使用OpenCV进行Weiner反卷积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！