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

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

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

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

任务清单

空域滤波

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

均值滤波

均值滤波是一种最简单的滤波方式。它以像素邻域内的平均值代替原先的像素:
I(x,y)=1N(u,v)NI(x+u,y+v)(1) \overline{I}(x,y)=\frac{1}{|N|}\sum_{(u,v)\in{N}} I(x+u,y+v)\tag{1}
均值滤波是线性的,即
I+J=I+J \overline{I+J}=\overline{I}+\overline{J}
。因此我们可以引入一个卷积核,用空域卷积来表示它,我们定义:
K(u,v)={1Nif(u,v)N,0else.(2) K(u,v)=\left\{ \begin{array}{rcl} \frac{1}{|N|} && {\text{if}(u,v)\in{N},}\\ 0 && {else.} \end{array} \right. \tag{2}
那么均值滤波又可以写成
I(x,y)=u,vK(u,v)I(x+u,y+v)(3) \overline{I}(x,y)=\sum_{u,v}K(u,v)I(x+u,y+v) \tag{3}

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

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

其中,图像卷积的计算可以通过 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>

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

高斯滤波

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

Filter Compare
(从上到下:原始图像、高斯滤波和用圆形邻域进行的均值滤波,可以看到灯光在均值滤波下形成了圆形的光斑。)

高斯滤波是另一种均匀的线性滤波器,它具有如下形式:
I(x,y)=u,vG(u,v)I(x+u,y+v)(4) \overline{I}(x,y)=\sum_{u,v}G(u,v)I(x+u,y+v) \tag{4}
其中 G(u,v)=12πσ2expu2+v22σ2G(u,v)=\frac{1}{2\pi \sigma^2}\text{exp}-\frac{u^2+v^2}{2\sigma^2}是二维高斯核函数,它的形状如下图:

Gaussian Kernel

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

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

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

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

中值滤波

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

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

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

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

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

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

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

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

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

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

双边滤波

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

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

I(p)=1WpqSGσs(pq)Gσr(I(p)I(q))I(q)(6) \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}

Wp=qSGσs(pq)Gσr(I(p)I(q))(7) W_p=\sum_{q\in S}G_{\sigma_s}(\|p-q\|)G_{\sigma_r}(|I(p)-I(q)|) \tag{7}

这里GσsG_{\sigma_s}是空间距离上的高斯加权,GσrG_{\sigma_r}是灰度距离上的高斯加权。

为什么要有WpW_p ?两个高斯核还需要12πσ2\frac{1}{2\pi \sigma^2}的归一项么?

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

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

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

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

傅里叶变换

大小为W*H的图像I的二维傅里叶变换定义为:
I^(ω1,ω2)=x=0W1y=0H1I(x,y)ej2π(xω1W+yω2H)(8) \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 的矩形,如下图:
Rect

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 of rect

你可以尝试对上同的图像进行 DFT 变换,可视化它的结果。

如果我们将矩形进行旋转,结果会怎么变化?平移呢?如果改变矩形的长和宽又会怎么样?

Tips

参考资料