去年写过手撸BMP文件头，直接操作图像文件的版本，这次因为课程作业，刚好就用OpenCV重新实现一遍。

实验内容

  • 实现盒状均值滤波
  • 实现高斯滤波
  • 实现中值滤波
  • 实现简单的双边滤波
  • 利用傅里叶变换完成图像的频域变换

理论及实验细节及效果

均值滤波

原理

均值滤波是一种线性滤波，它的卷积核在3x3的kernel中应该是：
$$
\begin{bmatrix}{1/9 ,1/9, 1/9 \1/9 ,1/9, 1/9 \1/9 ,1/9, 1/9 }\end{bmatrix}
$$

总之就是每个像素点滤波后的值是周围像素值的平均，即卷积核的值为$1/(w \times h)$.

实现

先是命令行参数的解析，后面每一个程序的基本结构都如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
int main(int argc, char** argv) {
	if (argc != 5) {
		cout << "Input Illegal!" << endl;
	}
	char* input = argv[1];
	char* output = argv[2];
	int w = atoi(argv[3]);
	int h = atoi(argv[4]);
	Mat mat = imread(input);
	Mat dst;
	myBoxFilter(mat, dst, 3, 3);
	imwrite(output, dst);
	return 0;
}

BoxFilter实现起来很简单，只需要把卷积核生成好，然后调用filter2D进行计算就可以了。

1
2
3
4
5
6
7
8
9
10
#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>
using namespace std;
using namespace cv;

void myBoxFilter(Mat& src, Mat& dst, int w, int h) {
	Mat mask(Size(h, w), CV_32FC1, Scalar(1.0 / w / h));
	filter2D(src, dst, src.depth(), mask);
}

结果

原图：

用3x3的均值滤波后的结果：

高斯滤波

原理

通常我们认为图像像素之间的相关性随着距离增加应该不断减弱。 均值滤波不能体现这一性质。在对图像进行均值滤波时，如果图像中有一些很显著的亮点，滤波后它的周围会形成光斑。这正是因为均值滤波无视了距离，对很远处的像素依旧采用同样的权重导致的。

因此，采用高斯卷积核，可以用距离来控制周围像素对中心的影响。

高斯卷积核理论上使无穷大的，但是距离远的点对于中心的影响比较小，同时为了方便计算，只取高斯核的中心部分进行卷积运算。

二维的高斯分布公式为：
$$
G(x,y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}
$$
可以由此公式构建高斯卷积核。

实现

先实现一个计算高斯分布的函数，此处x2相当于$x^2+y^2$

1
2
3
4
double G(double x2, double& sigma) {
	return (1.0 / 2 / PI / sigma) * exp(-1.0*(x2) / (2.0*sigma*sigma));
}

然后计算出高斯核中每个元素的值，并进行归一化。

此处默认取高斯核大小为5。

1
2
3
4
5
6
7
8
9
10
11
12
13
void myGaussianFilter(Mat& src, Mat& dst, double sigma) {
	int size = 5;
	int mid = size / 2;
	Mat mask(Size(size, size), CV_32FC1);
	for (int i = 0; i < size; ++i) {
		for (int j = 0; j < size; ++j) {
			mask.at<float>(i, j) = (float)G(pow(i - mid, 2) + pow(j - mid, 2), sigma);
		}
	}
	double s = sum(mask)[0];
	mask = mask / s;
	filter2D(src, dst, src.depth(), mask);
}

结果

原图：

取sigma=10时的结果：

中值滤波

原理

中值滤波是一种基于统计的非线性滤波。它取周围像素的种植作为当前的像素的值，可以有效地抵抗椒盐噪声的影响。

实现

因为属于非线性滤波，所以无法使用统一的线性卷积核来运算。

对于每一个点，都要取周围的点来排序取中值。

这里，我对于边缘的点采取的处理是仅取卷积核和图像重叠的部分参与运算。

对于每个点的操作是如下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int medianProcess(Mat& src, int y, int x, int w, int h, int ch) {
	vector<int> pixels;
	int mid_w = w / 2;
	int mid_h = h / 2;
	Size size = src.size();
	for (int i = y - mid_h; i <= y + mid_h; ++i) {
		for (int j = x - mid_w; j <= x + mid_w; ++j) {
			if (i >= 0 && i < size.height && j >= 0 && j < size.width) {
				pixels.push_back(src.at<Vec3b>(i, j)[ch]);
			}
			else {
				continue;
			}
		}
	}
	sort(pixels.begin(), pixels.end());
	return pixels.at(pixels.size() / 2);
}

然后，单独运算每个点的像素值就可以了。

1
2
3
4
5
6
7
8
9
10
11
void myMedianFilter(Mat& src, Mat& dst, int w, int h) {
	Size size = src.size();
	dst.create(size, src.type());
	for (int i = 0; i < size.height; ++i) {
		for (int j = 0; j < size.width; ++j) {
			dst.at<Vec3b>(i, j)[0] = medianProcess(src, i, j, w, h, 0);
			dst.at<Vec3b>(i, j)[1] = medianProcess(src, i, j, w, h, 1);
			dst.at<Vec3b>(i, j)[2] = medianProcess(src, i, j, w, h, 2);
		}
	}
}

效果

原图

取w=3,h=3的中值滤波结果：

双边滤波

原理

高斯滤波可以平滑图像，但是会模糊边界。为了可以保边，可以使用双边滤波。

在高斯滤波中，像素之间的欧式距离影响了其对于中心像素的影响。

在双边滤波中，不仅仅考虑欧式距离，还考虑像素值之间的距离的影响。对于像素值相差较远的点，认为此处是边界，所以其对于中心像素的影响较小。

将两个影响因素相乘，就是该像素对于中心像素的影响。最后进行归一化即可。

实现

对于每个点的操作如下：

对于每个点的周围mask范围，计算欧式距离，带入G函数，计算像素距离，带入G函数，作为该点的mask，进行卷积运算。

运算完后做一下归一化，保证影响的factor加起来等于1，使图片不会过亮或者过暗。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
double bilateralProcess(Mat& src, int y, int x, int kernel_size, int ch, double sigma_s, double sigma_r) {
	Mat mask = Mat(Size(kernel_size, kernel_size), CV_32FC1);
	Size size = src.size();
	double result = 0;
	double factor = 0;
	double total_factor = 0;
	Vec3b color = src.at<Vec3b>(y, x);
	for (int i = y - kernel_size / 2; i <= y + kernel_size / 2; ++i) {
		for (int j = x - kernel_size / 2; j <= x + kernel_size / 2; ++j) {
			if (i >= 0 && i < src.rows && j >= 0 && j < src.cols) {
				factor = G(pow(i - y, 2) + pow(j - x, 2), sigma_s) * G(pow(src.at<Vec3b>(i, j)[ch] - color[ch], 2), sigma_r);
				result += factor * src.at<Vec3b>(i, j)[ch];
				total_factor += factor;
			}
			else {
				continue;
			}
		}
	}
	return result / total_factor;
}

每个点的运算函数写好后，只要对于每个点进行运算即可：

1
2
3
4
5
6
7
8
9
10
11
12
void myBilateralFilter(Mat& src, Mat& dst, int kernel_size, double sigma_s, double sigma_r) {
	dst.create(src.size(), src.type());
	for (int i = 0; i < src.rows; i++)
	{
		for (int j = 0; j < src.cols; j++)
		{
			dst.at<Vec3b>(i, j)[0] = (char)bilateralProcess(src, i, j, kernel_size, 0, sigma_s, sigma_r);
			dst.at<Vec3b>(i, j)[1] = (char)bilateralProcess(src, i, j, kernel_size, 1, sigma_s, sigma_r);
			dst.at<Vec3b>(i, j)[2] = (char)bilateralProcess(src, i, j, kernel_size, 2, sigma_s, sigma_r);
		}
	}
}

结果

原图：

取sigma_s = 10, sigma_r = 10的结果：

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

原理

图像的存储是在空间域上的，可以通过傅里叶将其转换到频域并进行可视化。

实现

先把图像扩展到DFT效率最高的尺寸上，

然后用一个两通道的图像矩阵作为dft的结果。

然后调整dft结果的位置.

最后取对数，normalize之后可视化。

比较需要注意的是读入和输出的时候图像的type的转换。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
void myDFT(Mat& src, Mat& dst) {
	int m = getOptimalDFTSize(src.rows);
	int n = getOptimalDFTSize(src.cols);
	Mat padded;
	copyMakeBorder(src, padded, 0, m - src.rows, 0, n - src.rows, BORDER_CONSTANT, Scalar::all(0));

	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
	Mat complexMat;
	merge(planes, 2, complexMat);

	dft(complexMat, complexMat, DFT_COMPLEX_OUTPUT);

	fftshift(complexMat, complexMat);

	split(complexMat, planes);
	magnitude(planes[0], planes[1], planes[0]);
	Mat mag = planes[0];
	mag += Scalar::all(1);

	log(mag, mag);

	normalize(mag, mag, 0, 1, CV_MINMAX);
	mag.convertTo(dst, CV_8U, 255);
}

调整dft结果位置的代码给出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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);
}

效果

原图：

结果：