利用维纳滤波对运动模糊图像进行复原

一、概要
本文首先给出图像的维纳滤波器的推导过程用以补充书中缺失的证明，而后在C++的编译环境下，利用opencv构造维纳滤波器实现运动模糊图像的复原，并结合代码给出算法流程，最后在本文末给出全部代码。
二、图像的维纳滤波原理



三、构造退化函数
这部分冈萨雷斯所著的数字信号处理第三版已经给出了建模过程，因此这里直接引用书本上给出的运动模糊函数的推导过程。如下所示：


在这里应当注意：实际操作中往往会将图像的频谱中心化，因此进行图像模拟退化时所使用的退化函数也应当中心化。设原始图像大小为M·N，经过扩充后的大小为(2M)·(2N)。此时退化函数的尺寸也应当为(2M)·(2N)，且退化函数应当改写成(中心化)：

四、算法流程
1、对原始图像进行傅里叶变换，计算出原始图像f(x,y)的幅度谱和相位谱。即：

该部分代码如下所示：

	Mat src = imread("T.jpg");
	Mat gray_src;
	cvtColor(src, gray_src, CV_BGR2GRAY);
	Mat padded = Mat::zeros(2 * gray_src.rows, 2 * gray_src.cols, CV_32FC1);
	//图像扩充
	for (int row = 0; row < gray_src.rows; row++)
	{
   
		for (int col = 0; col < gray_src.cols; col++)
		{
   
			padded.at<float>(row, col) = gray_src.at<uchar>(row, col);
		}
	}
	//中心化
	for (int row = 0; row < padded.rows; row++)
	{
   
		for (int col = 0; col < padded.cols; col++)
		{
   
			padded.at<float>(row, col) *= pow(-1, row + col);
		}
	}
	//傅里叶变换
	Mat planes[] = {
    Mat_<float>(padded),Mat::zeros(padded.size(),CV_32FC1) };
	Mat complexImg;
	merge(planes, 2, complexImg);
	dft(complexImg, complexImg);
	//求原始图像的幅度谱和相位谱
	Mat temp1[] = {
    Mat_<float>(padded),Mat::zeros(padded.size(),CV_32FC1) };
	split(complexImg, temp1);
	Mat amplitude1, angle1;
	amplitude1 = Magnitude(temp1[0], temp1[1]);
	angle1 = Phase(temp1[0], temp1[1]);

2、构造退化函数h(x,y)，并计算其幅度谱和相位谱，即：

该部分代码如下所示：

	//定义退化函数
	Mat degenerate[] = {
    Mat_<float>(padded),Mat::zeros(padded.size(),CV_32FC1) };
	for (int row = 0; row < degenerate[0].rows; row++)
	{
   
		for (int col = 0; col < degenerate[0].cols; col++)
		{
   
			float v = pi*(0.1*(row - 262) + 0.1*(col - 175));
			if (v == 0)
			{
   
				degenerate[0].at<float>(row, col) = 1;
				degenerate[1].at<float>(row, col) = 0;
			}
			else
			{
   
				degenerate[0].at<float>(row, col) = sin(v) / v*cos(v);
				degenerate[1].at<float>(row, col) = -sin(v) / v*sin(v);
			}
		}
	}
	//求退化函数的幅度谱和相位谱
	Mat amplitude2, angle2;
	amplitude2 = Magnitude(degenerate[0], degenerate[1]);
	angle2 = Phase(degenerate[0], degenerate[1]);

3、实现图像的退化
代码实现的原理为：

随后将上式分离出来的实部与虚部分别赋予两个通道，在这里是用photo[0]存储实部，photo[1]存储虚部。
代码如下所示：

	//图像退化
	Mat photo[] = {
    Mat_<float>(padded),Mat::zeros(padded.size(),CV_32FC1) };
	for (int row = 0; row < padded.rows; row++)
	{
   
		for (int col = 0; col < padded.cols; col++)
		{
   
			photo[0].at<float>(row, col) = amplitude1.at<float>(row, col)*amplitude2.at<float>(row, col)*cos(angle1.at<float>(row, col) + angle2.at<float>(row, col));
			photo[1].at<float>(row, col) = amplitude1.at<float>(row, col)*amplitude2.at<float>(row, col)*sin(angle1.at<float>(row, col) + angle2.at<float>(row, col));
		}
	}

4、得到运动模糊图像
以上都是在频域中进行操作，现将上一步得到的结果进行傅里叶逆变换，然后进行显示。运动模糊的效果如下：

该部分代码如下：

	//空域中显示退化图像
	Mat complex;
	merge(photo, 2, complex);
	idft(complex, complex);
	Mat temp2[] = {
    Mat_<float>(padded),Mat::zeros(padded.size(),CV_32FC1) };
	split(complex, temp2);
	for (int row = 0; row < temp2[0].rows; row++)
	{
   
		for (int col = 0; col < temp2[0].cols; col++)
		{
   
			temp2[0].at<float>(row, col) *= pow(-1, row + col);
		}
	}
	for (int row = 0; row < temp2[0].rows; row++)
	{