图像复原之Lucy-Richardson算法滤波

算法原理见：图像去卷积（二）——Richardson–Lucy算法

     算法代码来自：Deconvolution with OpenCV?

     具体原理到代码实现过程有待进一步学习，这里先验证一下算法效果。

     程序和图片打包：Lucy-Richardson算法（opencv版） stackoverflow.rar

#include<iostream>
#include<opencv.hpp>
using namespace std;
using namespace cv;
static void RLTikh_deconvolution(Mat img, Mat &result, int num_iterations);
float sigmaG = 4.0, EPSILON =0.01;
int main()
{
    Mat guassianImage= imread("E://image//guassian.jpg");
    Mat guassianImageF, resImage7, resImage21, resImage150;
    guassianImage.convertTo(guassianImageF, CV_64FC3, 1.0 / 255, 0);
 
    imshow("guassianImageF", guassianImageF);
    //GaussianBlur(originImage, guassianImage, Size(11, 11), 4, 4);
    //imwrite("E://image//guassian.jpg", guassianImage);这两句是制造由于图像复原的特定参数的高斯模糊图像，
    //现在看此算法的sigmaG和winSize参数只有和高斯模糊图像中的GaussianBlur函数参数完全对应，才能取得较好效果，此例中为4和11.
    RLTikh_deconvolution(guassianImageF, resImage7, 7);
    RLTikh_deconvolution(guassianImageF, resImage21, 21);
    RLTikh_deconvolution(guassianImageF, resImage150, 150);
    imshow("restorationImage7", resImage7);
    imshow("restorationImage21", resImage21);
    imshow("restorationImage150", resImage150);
    waitKey(0);
    return 0;
}
static void RLTikh_deconvolution(Mat img, Mat &result, int num_iterations)
{
    // Lucy-Richardson Deconvolution Function
    // input-1 img: NxM matrix image
    // input-2 num_iterations: number of iterations
    // input-3 sigma: sigma of point spread function (PSF)
    // output result: deconvolution result
 
    // Window size of PSF
    int winSize = 11;// 10 * sigmaG + 1;
 
    // Initializations
    Mat Y = img.clone();
    Mat J1 = img.clone();
    Mat J2 = img.clone();
    Mat wI = img.clone();
    Mat imR = img.clone();
    Mat reBlurred = img.clone();
 
    Mat T1, T2, tmpMat1, tmpMat2;
    T1 = Mat(img.rows, img.cols, CV_64FC3, 0.0);
    T2 = Mat(img.rows, img.cols, CV_64FC3, 0.0);
 
    // Lucy-Rich. Deconvolution CORE
 
    double lambda = 0;
    for (int j = 0; j < num_iterations; j++)
    {
        if (j > 1) {
            // calculation of lambda
            multiply(T1, T2, tmpMat1);
            multiply(T2, T2, tmpMat2);
            lambda = sum(tmpMat1)[0] / (sum(tmpMat2)[0] + EPSILON);
            // calculation of lambda
        }
 
        Y = J1 + lambda * (J1 - J2);
        Y.setTo(0, Y < 0);
 
        // 1)
        GaussianBlur(Y, reBlurred, Size(winSize, winSize), sigmaG, sigmaG);//applying Gaussian filter 
        reBlurred.setTo(EPSILON, reBlurred <= 0);
 
        // 2)
        divide(wI, reBlurred, imR);
        imR = imR + EPSILON;
 
        // 3)
        GaussianBlur(imR, imR, Size(winSize, winSize), sigmaG, sigmaG);//applying Gaussian filter 
 
        // 4)
        J2 = J1.clone();
        multiply(Y, imR, J1);
 
        T2 = T1.clone();
        T1 = J1 - Y;
    }
 
    // output
    result = J1.clone();
}
--------------------- 
作者：Coming_is_winter 
来源：CSDN 
原文：https://blog.csdn.net/coming_is_winter/article/details/92800150 
版权声明：本文为博主原创文章，转载请附上博文链接！