OpenCV实现图像基础频率域滤波_opencv统计图像高频

写在前面： 刚开始接触数字图像处理频率域滤波时，很是陌生，感觉这个技术使用范围很窄，不如空域直接处理来的实在，最近看书发现有些情况又不得不在频率域中进行操作，个人感觉图像的复原与重建就是最大的应用点。特此实现一些基本的频率域滤波操作为后学习打下基础…

1. 频率域滤波步骤

前处理： 包括对图像边界填充，使之达到OpenCV傅里叶变换最佳尺寸，然后就是将$f(x,y)$乘以$-1^{(x+y)}$，使傅里叶变换位于填充后图像大小的频率矩形的中心。

后处理： 包括从反变换回来的图的左上象限提取一个大小为MXN(原图尺寸)的区域，得到与输入图像大小相同的滤波后的结果$g(x,y)$。

具体操作步骤详见《数字图像处理》第4章。

2.总结

1. $H(u,v)$是实函数，即虚部为0。
2. $G(u,v)$ = $H(u,v)\cdot F(u,v)$= $H(u,v)\cdot (R(u,v)+JI(u,v))$ = $H(u,v)\cdot R(u,v)+H(u,v)\cdot JI(u,v)$
3. 反变换获得图像$g(x,y)$取实部即可。
4. 在进行高通滤波时，滤波后的图像中会含有负值，这时直接将像素值归一化到[0, 1]，得到的结果同时标定了正值和负值，图像呈现深灰色，虽然也对，但是为了呈现另一种效果（底色为黑）我将负值进行了手动截断为0，然后在进行归一化。

3.程序实现

//***************************************tool***************************************//


//function:	使傅里叶频谱中心化
void fftshift(const cv::Mat& src, cv::Mat& dst)
{
	dst = src.clone();
	if (dst.type() != CV_32F) dst.convertTo(dst, CV_32F);
	
	for (int r = 0; r < dst.rows; ++r)
	{
		for (int c = 0; c < dst.cols; ++c)
		{
			if ((r + c) % 2 != 0)	//r + c 不为偶数，当前点变为负数
				dst.at<float>(r, c) = -dst.at<float>(r, c);
		}
	}
}


//function:	创建最佳傅里叶变换尺寸
cv::Size Make_BetterSize(const cv::Mat src)
{
	int w = cv::getOptimalDFTSize(src.cols);
	int h = cv::getOptimalDFTSize(src.rows);

	//width, height
	return Size(w, h);
}


//function:	将图像中负数截断为0
//opencv中应该有现成的工具函数，需要查一下文档
void MinusToZero(const cv::Mat& src, cv::Mat& dst)
{
	dst = src.clone();
	
	for (int r = 0; r < dst.rows; ++r)
	{
		for (int c = 0; c < dst.cols; ++c)
		{
			if (dst.at<float>(r, c) < 0)
				dst.at<float>(r, c) = 0.0;
		}
	}
}


//*****************************calculate filtering kernel*****************************//


//function:	创建理想低通滤波器
void calcILPF(const cv::Mat& src, int D0, Mat& dst)
{
	Size filterSize = Make_BetterSize(src);
	Mat H(filterSize, CV_32FC1, Scalar::all(0));

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//当距离小于截止频率时，设置为1
			if (d <= (double)D0)
				H_data[x] = 1.f;
		}
	}
	//改为双通道图像便于与F(u,v)相乘
	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}


//function:	巴特沃斯低通滤波器
void calcButterWorth(const cv::Mat& src, int D0, Mat& dst, float n)
{
	Size filterSize = Make_BetterSize(src);
	cv::Mat H = cv::Mat::zeros(filterSize, CV_32FC1);

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);		//指向第y行首个元素的指针
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//根据巴特沃斯低通滤波器计算H(u, v)
			H_data[x] = 1.0 / (1 + pow(d / D0, 2 * n));
		}
	}
	
	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}


//function:	创建高斯低通滤波器
void calcGLPF(const cv::Mat& src, int D0, Mat& dst)
{
	Size filterSize = Make_BetterSize(src);
	cv::Mat H = cv::Mat::zeros(filterSize, CV_32FC1);

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);		//指向第y行首个元素的指针
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//根据gauss低通滤波器计算H(u, v)
			H_data[x] = expf((-d * d) / (2 * D0 * D0));
		}
	}

	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}


//function:	创建高斯高通滤波器
//这里未做简化，只需要用1减去高斯低通核即可创建高通核
void calcGHPF(const cv::Mat& src, int D0, Mat& dst)
{
	Size filterSize = Make_BetterSize(src);
	cv::Mat H = cv::Mat::ones(filterSize, CV_32FC1);

	int cy = H.rows / 2;
	int cx = H.cols / 2;

	for (int y = 0; y < H.rows; ++y)
	{
		float* H_data = H.ptr<float>(y);		//指向第y行首个元素的指针
		for (int x = 0; x < H.cols; ++x)
		{
			double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
			//根据gauss低通滤波器计算H(u, v)
			H_data[x] = 1.f - expf((-d * d) / (2 * D0 * D0));
		}
	}

	Mat planes[] = { H.clone(), H.clone() };
	Mat complexH;
	merge(planes, 2, complexH);

	complexH.copyTo(dst);
}

//********************************执行G(u,v) = F(u,v) * H(u,v)*******************************//


//function:		将输入图像在频率域进行处理，然后输出
//InputImage:	输入单通道图像
//H:			频率域滤波核
void frequency_filter(const cv::Mat& InputImage, const cv::Mat& H, cv::Mat& OutputImage)
{
	cv::Mat padded;
	int m = getOptimalDFTSize(InputImage.rows);
	int n = getOptimalDFTSize(InputImage.cols);
	//on the border add zero values
	cv::copyMakeBorder(InputImage, padded, 0, m - InputImage.rows, 0, n - InputImage.cols, BORDER_CONSTANT, Scalar::all(0));
	cv::Mat planes[] = { Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F) };
	//使傅里叶变换中心化
	fftshift(planes[0], planes[0]);
	//创建一个复数图像F(u, v)
	cv::Mat complexF;
	cv::merge(planes, 2, complexF);
	//执行傅里叶变换
	cv::dft(complexF, complexF);

	cv::Mat complexFH;
	cv::Mat ifft;
	//采用对应像素相乘G(u, v) = H(u, v) * F(u, v)
	cv::multiply(complexF, H, complexFH);
	cv::idft(complexFH, ifft, DFT_REAL_OUTPUT);
	
	//最后不要忘了再完成一次移动
	fftshift(ifft, ifft);
	//为了显示效果将负数截断为0
	MinusToZero(ifft, ifft);
	normalize(ifft, ifft, 0, 1, NORM_MINMAX);
	ifft.convertTo(ifft, CV_8UC1, 255);
	//裁剪回原尺寸输出
	ifft(cv::Rect(0, 0, InputImage.cols, InputImage.rows)).copyTo(OutputImage);
}


//*********************************main()*******************************//

int main()
{
	string path = "blurring_test.tif";

	Mat SrcImage = imread(path, IMREAD_GRAYSCALE);
	if (!SrcImage.data) {
		std::cout << "Could not open or find the image" << std::endl;
		return -1;
	}

	//高通滤波结果和高通核
	Mat H_outimg, H_kernel_H;
	//低通滤波结果和低通核
	Mat L_outimg, L_kernel_H;

	calcGLPF(SrcImage, 60, L_kernel_H);
	frequency_filter(SrcImage, L_kernel_H, L_outimg);

	calcGHPF(SrcImage, 80, H_kernel_H);
	frequency_filter(SrcImage, H_kernel_H, H_outimg);
	
	imwrite("low_blurring_test.png", L_outimg);
	imwrite("high_blurring_test.png", H_outimg);
	
	cv::waitKey(0);
	return 0;
}
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232

4. 测试结果

原图：

高斯低通滤波：

高斯高通滤波：