Lab 2 - 图像滤波和傅里叶变换

本次实验我们利用 OpenCV 进行简单的图像滤波。

本次实验内容为课程作业，计算成绩。你需要将 源代码、可执行程序（注明运行环境）和实验文档 打包并上传至 学在浙大 学*，压缩包名称为 Lab2-学号-姓名.zip/7z 。

本次作业提交截止时间： 2022年3月14日 23:59:59 ，逾期将要扣分。

实现盒状均值滤波

实现高斯滤波

实现中值滤波

实现简单的双边滤波

利用傅里叶变换完成图像的频域变换

空域滤波的基本思想是取像素的邻域，将这个邻域内的全部图像信息进行整合以得到对应像素的估计。

均值滤波是一种最简单的滤波方式。它以像素邻域内的平均值代替原先的像素：
$\overline{I}(x,y)=\frac{1}{|N|}\sum_{(u,v)\in{N}} I(x+u,y+v)\tag{1}$

从上面的式子可以看到，这里的卷积核K与图像域坐标(x,y)无关，因此我们说 均值滤波是线性的 。

在 OpenCV 中，我们有两种基本方法实现均值滤波：

直接依照(1)实现；

计算(2)的卷积核，然后按照(3)进行卷积。

其中，图像卷积的计算可以通过 cv::filter2D 完成

图像通常是有边界的，这里我们并没有考虑边界的情形。对于边界，常见的处理方法有下面几种：

填充固定像素 iiiiii|abcdefgh|iiiiiii

重复靠近边缘的某像素 aaaaaa|abcdefgh|hhhhhhh

重复边缘的倒影 fedcba|abcdefgh|hgfedcb

可以点击 cv::BorderType 了解各种常见的边界处理方法。
cv::filter2D 函数已经处理了边界，本次实验无需再手工处理边界。

除了手工实现外，如果我们的卷积核是矩形，还可以使用 OpenCV 的 cv:: boxFilter 使用均值滤波器。

本次实验需要你 利用卷积方法 完成图像的均值滤波。你需要手工实现一个名为 BoxFilter 的可执行程序，用法如下：

BoxFilter <input-image> <output-image> <w> <h>

它使用矩形的卷积核进行均值滤波，其中， w 和 h 是卷积核中心点到边界的距离，卷积核矩阵的大小应当为 (w*2+1)*(h*2+1) 。

通常我们认为图像像素之间的相关性随着距离增加应该不断减弱，但是(1)的均值并没有体现这一性质。在对图像进行均值滤波时，如果图像中有一些很显著的亮点，滤波后它的周围会形成光斑。这正是因为均值滤波无视了距离，对很远处的像素依旧采用同样的权重导致的。一些场合，我们为了美感会需要这种效果。另一些场合，这种结果是不利的，特别是在做图像处理时。

（从上到下：原始图像、高斯滤波和用圆形邻域进行的均值滤波，可以看到灯光在均值滤波下形成了圆形的光斑。）

高斯滤波是另一种均匀的线性滤波器，它具有如下形式：
$\overline{I}(x,y)=\sum_{u,v}G(u,v)I(x+u,y+v) \tag{4}$

本次实验需要你实现的第二个任务是高斯滤波，结合前面均值滤波的知识，你需要做的是生成高斯卷积核，然后用它对图片进行滤波。你实现的程序用法如下：

GaussianFilter <input-image> <output-image> <sigma>

其中， sigma 代表了高斯核的参数 ，是一个浮点数。高斯核理论上大小是无穷的，你可以只保留中心 $5\sigma$ 的大小，也就是卷积核矩阵大小为 $(2*\lfloor 5\sigma \rfloor + 1) \times (2*\lfloor 5 \sigma \rfloor + 1)$

OpenCV 提供了 cv::GaussianBlur 函数可以用于高斯滤波，本次实验需要你自己建立高斯卷积核，然后使用 cv::filter2D 进行卷积。

高斯滤波适用于图像带有高斯噪声情况下的去噪。在图像被非高斯噪声污染的情况下，高斯滤波不一定能得到理想的去噪效果。

中值滤波是一种序统计滤波器（Order-Statistic Filter），序统计滤波器是依据邻域的值在统计上的次序关系来进行过滤的。这里，中值滤波器用邻域内像素亮度的中值来取代原本的像素值，即：
$\overline{I}(x,y)=\text{Median}\{I(x+u,y+v)|(u,v)\in N\} \tag{5}$

中值是将样本集合划分为相同大小的两部分的样本点的值。

由定义可见中值滤波器是非线性的。

椒盐噪声（随机的01噪音）不是高斯的，利用前面的卷积滤波方法不能得到很好的结果。中值滤波对于椒盐噪声有比较好的抑制效果。直观上看，在邻域内的样本中，椒盐噪声平均地分布在最大和最小端，通过取中值可以较好地过滤椒盐噪声。

（中值滤波可以很好地过滤椒盐噪声，第二行的两幅图分别采用 3x3 邻域和 9x9 邻域，更大的邻域使得图像边界变得圆滑。）

本次实验的第三个任务需要完成 盒状邻域 的中值滤波。你需要实现名为 MedianFilter 的程序，用法如下：

MedianFilter <input-image> <output-image> <w> <h>

参数 w 和 h 的含义与均值滤波中对应参数的含义相同。

OpenCV 提供了 cv::medianBlur 函数可以用于中值滤波，本次实验需要你自己实现中值滤波，不应当使用该函数。

前面的三种滤波器都会破坏图像的边界，在卷积核很大的时候，均值滤波和高斯滤波都会让边界变得模糊，在邻域很大时，中值滤波会减小边界的曲率。由于物体边界是物体的一个重要特征，很多任务里我们不希望图像边界被破坏。

双边滤波提供了一种降噪同时保持边界的方法。它的思路很简单：如果邻域内像素的亮度差异很大，它在加权平均时的贡献也应当小。我们可以在高斯滤波加权平均的基础上引入一个新的项，反应亮度差带来的加权：

$\overline{I}(p)=\frac{1}{W_p}\sum_{q \in S}G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I(p)-I(q)|)I(q) \tag{6}$

$W_p=\sum_{q\in S}G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I(p)-I(q)|) \tag{7}$

这里 $G_{\sigma_s}$

为什么要有 $W_p$

本次实验需要你理解双边滤波的思想，并完成双边滤波。你需要实现名为 BilateralFilter 的程序，用法如下：

BilateralFilter <input-image> <output-image> <sigma-s> <sigma-r>

参数 <sigma-s> 和 <sigma-r> 分别对应 ${\sigma_s}$

OpenCV 同样提供了 cv::bilateralBlur 函数可以用于双边滤波，你不应当在本次作业中使用该函数。

傅里叶变换

大小为W*H的图像I的二维傅里叶变换定义为：
$\hat {I} (\omega_1,\omega_2)=\sum_{x=0}^{W-1}\sum_{y=0}^{H-1}I(x,y)e^{-j2\pi(\frac{x\omega_1}{W}+\frac{y\omega_2}{H})} \tag{8}$

上式定义在复数域中，因此傅里叶变换的输入和输出都是复数域上的。但我们用实数域表示图像，因此在操作时要额外注意类型：变换到频域时，采用复数表达，变换回图像时，只保留实数部分。

在 OpenCV 中，图像的傅里叶变换可以用 cv::dft 函数完成。下面的代码展示了如何进行图像的傅里叶变换：

Mat dft_result(image.size(), CV_32FC2);
dft(image, dft_result, DFT_COMPLEX_OUTPUT); 
当输入/输出的矩阵大小满足特殊条件时，cv::dft  函数性能最佳，参考 cv::getOptimalDFTSize。
变换后的结果中，左上角对应了零频率。在操作时我们要使用这种格式，但在可视化时我们希望将左上角置于中心。使用下面的函数  fftshift  将零频率移到中央：
// rearranges the outputs of dft by moving the zero-frequency component to the center of the array.
void fftshift(const Mat &src, Mat &dst) {
	dst.create(src.size(), src.type());
	int rows = src.rows, cols = src.cols;
	Rect roiTopBand, roiBottomBand, roiLeftBand, roiRightBand;
	if  (rows %  2  ==  0) {
		roiTopBand =  Rect(0,  0, cols, rows /  2);
		roiBottomBand =  Rect(0, rows /  2, cols, rows /  2);
	} else {
		roiTopBand =  Rect(0,  0, cols, rows /  2  +  1);
		roiBottomBand =  Rect(0, rows /  2  +  1, cols, rows /  2);
	if  (cols %  2  ==  0)  {
		roiLeftBand =  Rect(0,  0, cols /  2, rows);
		roiRightBand =  Rect(cols /  2,  0, cols /  2, rows);
	}  else  {
		roiLeftBand =  Rect(0,  0, cols /  2  +  1, rows);
		roiRightBand =  Rect(cols /  2  +  1,  0, cols /  2, rows);
	Mat srcTopBand =  src(roiTopBand);
	Mat dstTopBand =  dst(roiTopBand);
	Mat srcBottomBand =  src(roiBottomBand);
	Mat dstBottomBand =  dst(roiBottomBand);
	Mat srcLeftBand =  src(roiLeftBand);
	Mat dstLeftBand =  dst(roiLeftBand);
	Mat srcRightBand =  src(roiRightBand);
	Mat dstRightBand =  dst(roiRightBand);
	flip(srcTopBand, dstTopBand,  0);
	flip(srcBottomBand, dstBottomBand,  0);
	flip(dst, dst,  0);
	flip(srcLeftBand, dstLeftBand,  1);
	flip(srcRightBand, dstRightBand,  1);
	flip(dst, dst,  1);
为了简化实验，我们只使用单通道图像。我们首先创建一幅 512x512 的单通道浮点数图像  I，在中心放置一个 20x60 的矩形，如下图：

Mat I(512, 512, CV_32FC1);
I = 0;
I(Rect(256-10, 256-30, 20, 60)) = 1.0; 
我们将其进行傅里叶变换，得到 512x512 的复（CV_32FC2）矩阵J，然后利用  fftshift  函数将零频率放置到中心。
Mat J(I.size(), CV_32FC2);
dft(I, J, DFT_COMPLEX_OUTPUT);
fftshift(J, J); 
我们计算  J  的每个元素的模长：
Mat Mag;
vector<Mat> K;
split(J, K); // 将实数和虚数部分分解到 K[0] 和 K[1] 
pow(K[0], 2, K[0]); // 计算平方 
pow(K[1], 2, K[1]);
Mag = K[0]+K[1]; // 两个分量的平方和  
最后，我们将  可视化出来：
Mat logMag;
log(Mag+1, logMag);
normalize(logMag, logMag, 1.0, 0.0, NORM_MINMAX);
// ... 
imshow("Magnitude", logMag); 
你应当得到如下结果：

你可以尝试对上同的图像进行 DFT 变换，可视化它的结果。
如果我们将矩形进行旋转，结果会怎么变化？平移呢？如果改变矩形的长和宽又会怎么样？
在需要手动处理边界问题时，可以先cv::copyMakeBorder扩大边界，以使得卷积核可以处理图像最边缘的像素，在处理完之后，使用cv::Rect把原图像区域取出来。
[1] R. C. Gonzalez and R. E. Woods, Digital Image Processing, Third Edition.
[2]  Bilateral Filtering for Gray and Color Images
[3] Making your own linear filters!