数字图像处理小作业2：图像扭曲与修复（高斯模糊、高斯噪声、维纳滤波）

数字图像处理小作业2

代码仓库

https://github.com/SingleZombie/WienerFilteringCpp

需求分析

老师给的需求如下：

• 读入一幅灰度图像
• 生成模糊图像：高斯模糊滤波+噪声
• 利用维纳滤波去模糊

总结来说，本次作业处理的图像对象是灰度图。灰度图要经过扭曲和修复两个过程。扭曲是高斯模糊和噪声（我选择了高斯噪声），修复用的是维纳滤波的方法。

相比上次的作业，本次需求清晰了不少。

技术学习

高斯模糊

模糊一个图像，可以通过卷积核来实现。具体来说，就是把一个像素的值加上它旁边几个像素的值。这样从直观上来看，像素之间的差异被减小了，图片就变模糊了。

问题是，一个像素值该如何加上旁边几个像素的值？每个像素相加的权重是什么？有人把二维高斯分布函数当成概率密度函数，权重就是概率密度函数的值。这样的利用高斯函数来计算卷积核权重的模糊方法就是高斯模糊。

高斯噪声

噪声就是把某个位置的像素值改变多少看成一个随机变量。如果随机变量符合高斯分布，那么这个噪声就是高斯噪声。随机变量和位置无关，只和随机变量自己的性质（均值、方差）有关。

维纳滤波

图像修复就是假设一个扭曲模型（输入为原图像的一个函数），之后设法求出这个模型中的参数，再利用这些参数复原图像。

维纳滤波就是假设原图像经过了一个卷积和一个噪声添加。该方法要的思想是让复原出的图像和原图像的”差异“尽可能小，这个”差异“被定义为均方误差。由于我没有学过信号处理这门课，而本课程又更关注算法的使用而非原理，这里不加证明地给出最终的式子：

$$\hat{F}(u,v) = [\frac{1}{H(u,v)} \cdot \frac{|H(u, v)|^2}{|H(u, v)|^2 + k}]G(u,v)$$

其中$\hat{F}(u,v)$是复原出来的图像的频域值，$H(u,v)$是退化函数在频域下的表示，$G(u,v)$是被扭曲过的图像。而$k$是一个可调的参数。

本次任务的难点应该就是调试$k$这个参数。

有一篇博客对维纳滤波的原理以及matlab下的实现做了讲解。

OpenCV用法

写着写着发现又碰到了很多不会用的函数。

cv::normalize

cv::normalize(src, dest, upper_bound = 1, lower_bound = 0, norm_type = cv::NORM_L2);

normalize用于归一化，避免一个向量的模变大。在图像中，卷积核做归一化可以保证图像的亮度不变。

我只会用这个函数的前5个参数，有几个参数的名字是我取的。

src：原图像

dest：输出图像

upper_bound：归一化后的上界

lower_bound：归一化后的下界

norm_type：归一化的范式

归一化的基本操作就是要让模为1，因此输出的下界是0，上界是1。归一化的操作是把向量的每个分量都除以向量的模。

范式就是模是如何计算的。第一范式把模计算成所有分量之和，第二范式把模计算成所有分量平方和开方…………无穷范式把模计算成最大分量。

cv::filter2D

cv::filter2D(srcMat, destMat, ddepth, kernel, anchor = (-1, -1))

这个函数用于做图像的卷积，省去了手写循环的无聊过程。

同样，我只列出了我会的参数，参数名字是随意写的。

srcMat：原图像

destMat：输出图像

ddepth：输出图像深度，-1为和原图像相同

kernel：卷积核

anchor：锚点，也就是被卷积的点在卷积核的哪个位置。(-1, -1)是卷积核正中心

cv::dft

cv::dft(srcMat, destMat, flag)

这个函数就是做离散傅里叶变换，内部实现方法肯定是用FFT做的。

srcMat：原图像

destMat：输出图像

flag：操作的一些参数。我知道的参数有：

​ cv::DFT_CCOMPLEX_OUTPUT: 输出结果以复数形式存储，通道数增加。我输入单通道的图像出来是双通道的。

​ cv::DFT_REAL_OUTPUT：仅保留结果复数的实数部分，通道数不变。

​ cv::DFT_SCALE：运算结束后除以一个数，用于IFFT中

​ cv::DFT_INVERSE：做的是逆变换。

在我的程序里，FFT的写法是cv::dft(src, dest, cv::DFT_COMPLEX_OUTPUT);，IFFT的写法是cv::dft(src, dest, cv::DFT_SCALE | cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT);

C++标准库用法

本次程序中我竟然碰到了完全没见过的c++标准库函数。

std::normal_distribution

#include <random>

std::random_device divice;
std::default_random_engine rng(divice());
std::normal_distribution<> normDis(mu, sigma);

double randomNumber = normDis(rng);

正态分布函数在c++的随机数库中。使用c++的随机数库有一些麻烦。语句前两行貌似是一个分布函数的初始化步骤，我不求甚解地直接使用了。语句第三行用平均值和方差初始化一个正态分布。初始化了之后每次就能获取一个符合正态分布的值了。

用OpenCV完成的组合操作

本节讲的是完成某一功能的，使用到OpenCV的自定义函数。

我在完成这个作业的时候，参考了很多别人的博客。但是有些人用的是matlab实现的，里面出现了C++OpenCV里没有的函数。我只好去搜了一些这些函数的原理，参考了别人的C++实现并自己实现了一下。

circShift

circShift(mat, x, y)

该函数的作用是，把一个矩阵mat向x轴正方向移动x个单位，向y轴正方向移动y个单位。移动中溢出的数字会填补到空出来的地方，也就是整个矩阵是循环移动的。

该函数实现起来很简单，先把x,y变成模意义下的值，再用临时矩阵做比较繁琐的copyTo就行了。

cv::Mat circShift(const cv::Mat& mat, int x, int y)
{
	int mWidth = mat.cols;
	int mHeight = mat.rows;
	x = (x % mWidth + mWidth) % mWidth;
	y = (y % mHeight + mHeight) % mHeight;

	cv::Mat res(mat.size(), mat.type());
	cv::Mat m1(mat, cv::Rect(0, 0, mWidth - x, mHeight));
	cv::Mat m2(mat, cv::Rect(mWidth - x, 0, x, mHeight));
	m1.copyTo(res(cv::Rect(x, 0, mWidth - x, mHeight)));
	m2.copyTo(res(cv::Rect(0, 0, x, mHeight)));
	
	cv::Mat res2(mat.size(), mat.type());
	m1 = cv::Mat(res, cv::Rect(0, 0, mWidth, mHeight - y));
	m2 = cv::Mat(res, cv::Rect(0, mHeight - y, mWidth, y));
	m1.copyTo(res2(cv::Rect(0, y, mWidth, mHeight - y)));
	m2.copyTo(res2(cv::Rect(0, 0, mWidth, y)));

	return res2;
}

由于我对OpenCV的用法不算很熟，没有写性能更好的方法，我感觉我的写法copyTo用得有点多了。我写的时候比较匆忙，只顾着实现，没有考虑到性能的问题。

psf2otf

psf2otf(mat, outSize = Size(-1, -1))

该函数的作用是，把一个二维点扩散函数转换成指定大小频域形式。

点扩散函数也许在数学上有其他定义，但根据我对该函数效果的理解，这里的点扩散函数指的就是一个用矩阵表示的二维函数，矩阵某位置的值表示二维函数在此处的函数值，该矩阵中心点是原点$(0, 0)$。

该函数的直接作用非常不好理解，但知道该函数的应用后就容易理解了。一个图像用一个卷积核卷积在时域上的计算可能比较慢，但是把图像和卷积核转换到到频域后，只要相乘，再转换回时域就可以得到卷积的结果了。由于卷积核的原点在中心点而不在左上角，且其大小可能比图像小，因此在卷积操作前需要把卷积核调整位置，并且扩大到和原图像一样的大小。经过处理后的卷积核的频域表示就能和原图像的频域表示进行正确的计算了。

该函数作用的宏观解释在明白函数的应用后，就比较容易理解了。但要理解其具体实现方法，需要了解FFT本身原理。我是通过理解一维FFT的原理后推广到二维的。

这里不加解释地直接给出该函数的具体操作了。该函数先把卷积核扩展至指定大小，再把卷积核的中心移动到左上角（准确来说是原点，计算机里左上角是原点）。注意这一步移动是循环移动，卷积核的部分内容被循环移动到了 右上角、左下角，右下角。扩展和移动结束后，再做FFT就行了。

cv::Mat psf2otf(const cv::Mat& mat, cv::Size outSize)
{
	if (outSize.height == -1 || outSize.width == -1)
	{
		outSize = mat.size();
	}

	cv::Mat otf;
	cv::Size paddleSize = outSize - mat.size();
	cv::copyMakeBorder(mat, otf,
		paddleSize.height - paddleSize.height / 2, paddleSize.height / 2,
		paddleSize.width - paddleSize.width / 2, paddleSize.width / 2,
		cv::BORDER_CONSTANT,
		cv::Scalar::all(0));

	otf = circShift(otf, -otf.size().width / 2, -otf.size().height / 2);
	otf = grayImageFFT(otf);
	
	return otf;
}

我在写这个函数的时候参考了这篇博客。这篇博客的实现 在卷积结束以后，还判断了虚数部分是否比较小，可以忽略。但我没有多加这一步判断。

图像与卷积核在频域的计算

这部分没什么新的内容，就是把代码组合了一下：

// Get frequent field of origin image
cv::Mat fOrigin = grayImageFFT(origin);

// Get frequent field of kernel and rearrange it
cv::Mat h = getKernal(..........);
cv::Mat funcH = psf2otf(h, mat.size());

// Execute complex number multiply
for (int i = 0; i < fOrigin.rows; i++)
{
	for (int j = 0; j < fOrigin.cols; j++)
	{
		std::complex<double> g(fOrigin.at<cv::Vec2f>(i, j)[0], fOrigin.at<cv::Vec2f>(i, j)[1]);
		std::complex<double> h(funcH.at<cv::Vec2f>(i, j)[0], funcH.at<cv::Vec2f>(i, j)[1]);

		g *= h;

		fOrigin.at<cv::Vec2f>(i, j)[0] = g.real();
		fOrigin.at<cv::Vec2f>(i, j)[1] = g.imag();
	}
}

// Do IFFT
cv::Mat matRes = grayImageIFFT(fOrigin);
cv::Mat result = cv::Mat(matRes, cv::Rect(0, 0, origin.cols, origin.rows));

以上代码算半个伪代码，完成了普通的频域相乘操作。如果要做其它的操作（比如维纳滤波），改循环里的公式就行了。

结构设计

同上次任务。

其中图像处理分成两个模块：退化图像生成与退化图像修复。

退化图像生成又分两个模块：模糊和加噪声。

最终实现的时候又加了一个交互模块，可以动态调参。这个模块要调用退化图像生成模块和退化图像修复模块。模块化设计的好处在这里就体现出来了。

虽然我感觉我模块分得特别好，最终代码清晰简明，符合数据流的内在逻辑结构，但我还是不太想在这里炫耀了。

程序设计心路历程

我看完高斯模糊和高斯噪声的原理后，很快就把这两部分实现了。

在了解原理，知道了C++有正态分布生成器后，这部分应该非常容易实现。

先提一句，最开始高斯模糊卷积我是直接拿filter2D做的。

问题就出在维纳滤波上。

一开始，我把原图像做了FFT，卷积核直接做了FFT再扩展到和原图像同样大小，再逐像素用公式计算，最后把计算结果IFFT。最终，我得到了一幅乱七八糟的图像。

频域的操作难以可视化，难以调试。我甚至都不知道我哪一步错了。但是，我凭借着我冷静的大脑，准备把这个过程逐一调试。

首先，我打算验证FFT和IFFT是否正确。我对一个图像进行FFT再IFFT回来，发现图像没有复原。果然，FFT和IFFT这步我就做错了。为了解决这个问题，我查了好多东西，还瞎照网上写了个FFT的优化函数。最后我理解了一下cv::dft的参数的含义，才把FFT和IFFT做对。

接着，我开始怀疑我的卷积核是不是不太对。FFT的原点在左上角，而卷积核的原点在中间，这好像不太对劲。但我一时半会儿解决不了这个问题。于是我干脆把卷积核的原点改到了左上角，即每次卷积的时候只对右下角的部分卷积，忽略掉这个问题，先解决其他问题。

我现在想验证除了卷积核外，我整个维纳滤波公式的运算是否有问题。我知道，如果没有噪声的话，逆滤波可以直接把图像复原，此时维纳滤波k值取0。因此，我暂时去掉了噪声，把k值设成0，看算法效果如何。结果图像确实复原出来了，但四周有一些不和谐的条状噪声。这个结果搞得我一头雾水。

我进一步思考，逆滤波能够复原图像，是因为噪声卷积操作和逆变换互为逆操作。但我卷积是拿filter2D做的，逆变换是在频域做的，这二者会不会工作原理不一样？于是，我暂时把卷积操作也变成频域操作。

终于，经过我这次的瞎折腾后，露娜sama美丽的面庞终于清晰地出现在了屏幕上。我又对退化图像加上了噪声。这次，k值取0的时候画面一片模糊，但k值只要稍稍增大，依旧能生成较退化图像更清晰的复原图像。至此，调试第一阶段结束，我高兴了好一会儿。

顺带一提，这个阶段为了调试方便，我把调k值放进了GUI里面，这样可以动态看到调参的结果。由于Debug模式下性能捉鸡，我又把OpenCV的Release版本编译了一遍。

问题是，现在我忽略了两个问题，一个是卷积核中心位置的问题，一个是用Filter2D的问题。我准备先解决第一个问题。

我把卷积核的原点又放回了中心。这次修改后，复原图像倒是能够生成了，但是图像的位置完全乱了：图像的左上角在右下角，右下角在左上角……。整个图像的四周都颠倒了。我去写了个OpenCV的fft2shift，试图把图像移正。但移正之后发现图像复原是失败的，最终图像和退化图像长得差不多。

在反复阅读了网上的代码后，我发现我代码唯一不太对的地方就是少了个psf2otf函数。我去搜了一下这个函数，竟然惊喜地搜到了这个函数的C++实现。我赶紧学了一下，看懂了函数的工作原理。我把这个函数加入了代码中，但怎么调都调不对。无奈之下，我只好双手离开键盘，选择冷静思考一波。

吃饭的时候，我仔细想了一想FFT的原理。我通过我以前对一维FFT的认识，以及昨天搜索到的零碎知识，渐渐搞清了psf2otf的原理以及目的。原来，图像FFT溢出的部分会“往回填”，图像靠右的部分表示的是负频率。为了保证卷积操作的正确进行，卷积核也要在扩充至和图像一样大小后，循环移动到左上角。

理清楚了每一条逻辑后，我修改了代码，在卷积核原点在中心时，图像终于能够正常地复原了。

最后是Filter2D的问题。如果卷积操作是拿这个函数做，得到的复原图像四周就会有条纹。我开始分析这种问题出现的原因。从“只有四周有问题”这个事实上，我想到很可能是Filter2D对边界处理的方式和用FFT处理不同。FFT是把整幅图像当成一个循环的图像来处理，因此会出现最右边像素受到最左边的像素的影响这种问题。而Filter2D对边界的处理有多种方式。我查了一下函数说明，果然Filter2D默认是把边缘反射处理，和在频域运算的方法不相同。但我尝试把边缘循环处理的时候，函数却报错了。源代码Assert掉了循环处理，不允许使用。无奈之下，我只好不使用Filter2D。反正我已经完全理解了这些操作的细节了。

最终程序的详细内容我就不讲了。感觉这份程序结构化做得非常棒，功能实现得非常完全。由于是作业，我还不得不加上了很多注释。相信大部分人都能看懂我那清晰易懂的代码。具体代码可以去代码仓库查看。

处理结果

原图

噪声图

修复结果

附原图的彩图版本：

感想

这个作业实在是太烦人了。原来我以为一个下午就能做完的东西，硬是做了两天。

我遇到的问题实在是有些令人恼火。我既不是程序代码不会写，也不是算法原理没搞懂，而是API怎么调，某个操作的实现细节这些边界知识没有掌握。本次作业没增加多少我对复原算法的认识，倒是进一步强化了我搜索信息和理解信息的能力。

之前，包括课本在内，很多介绍模糊和复原的文章都在用一个女性照片做为原图像。我一开始还没明白为什么。后来，在写程序碰到瓶颈之余，我好好分析了一下这个问题，总算得到了一个非常漂亮的结论：科研工作者大多为单身男性，每日面对的只有无情的代码和冷冰冰的文字。在这种恶劣的条件下，偶尔出现的女性照片，能给你黯淡的内心点上一丝救赎的光芒。所以，人们更多地选择那张女性照片。

可惜我的内心对于这种照片不为所动。我想，凭什么只用三次元的照片？于是，我随便在电脑里找了一张露娜sama的图片。果然，看着这幅CG被一步一步复原，我被各种BUG蹂躏得受伤的心灵得到了安慰。我为自己的机智和随机应变能力而赞叹不已。