循序漸進之（十）圖像複原之Lucy-Richardson算法濾波

                   循序漸進之（十）圖像複原之Lucy-Richardson算法濾波（OpenCV版）

     算法原理見：圖像去卷積（二）——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();
}
           

      實作結果：