Abstract: 数字图像处理：第23天
Keywords: 频域滤波，频域滤波器，振铃，卷积
本文最初发表于csdn，于2018年2月17日迁移至此
开篇废话
这两天写了一下频域滤波的代码，并且发现以前博客里代码的一个BUG，产生BUG的原因是一维存储的二维矩阵行宽和列长被写反了，而且测试没有测出来，看来测试是必要的，而且还学习了下使用github，托管代码，因为是小菜，托管HelloWorld也不丢人，频域滤波打算写两到三篇，然后结束，进行下一大块问题的研究。
频域滤波的主要知识框架作了一个结构图，完全个人理解，如有理论上的问题，请及时指出：
频域滤波基础
对于空域滤波和时域滤波，其最重要的理论基础就是，频率域的高频部分对应于空域中变化陡峭的细节，而频域的低频部分对应于空域变化换面的平坦区域，为了得到平滑或陡峭的部分，我们就想到了频域滤波，在空域中对细节和非细节的提取和定义并没有频域那么直接与明确，通过频域这个实验室，可以得出很多针对不同问题，特殊的滤波器，然后将其转换到空域，生成小模板，使用卷积运算来完成快速的滤波，达到我们的目的。基本关系示意图：
没错，冰冰和凤姐就差一个滤波器！（此为示意图，没有这样的滤波器，请勿模仿0.0）
我们的基本原理就是高频对应变化陡峭的部分，低频对应变化缓慢的部分。低频和高频在频率域以数字度量，滤波器对频域的数字进行不同程度的增强或抑制，之后进行逆变换，得到想要的图像，频域滤波器必须在频域关于原点对称，不然对称部分频率的逆变换将是原图产生伪影效果。
影响频域滤波效果的只有输入和滤波器，对于未处理的输入，如果直接进行DFT和滤波的话，会产生纠缠效果，我们来看数学效果:
卷积定理表明：对于3x3的两个矩阵a和b的卷积，应该等于a的DFT和b的DFT对位相乘（对位相乘的意思就是对于a和b对位相乘的结果c，c(x,y)=a(x,y)*b(x,y)，而非一般的矩阵相乘），再进行IDFT变换，就可以得到和空域的卷积结果，但我们得到的结果是：1
2
3
4
5
6
7
8
9
10
11
12clear all;clc;
a=[1 2 1;1 1 1 ;0 2 1];
fa=fft2(a);
b=[3 2 1;3 1 1 ;2 0 1];
fb=fft2(b);
temp=eye(3);
for R=1:3
for C=1:3
temp(R,C)=fa(R,C)*fb(R,C);
end
end
res=ifft2(temp);
卷积 
等于结果： 
通过我们用手计算，发现结果不对，这是什么原因呢，我们暂时不说原因，我们把矩阵用0扩充到5x5：1
2
3
4
5
6
7
8
9
10
11
12clear all;clc;
a=[1 2 1 0 0;1 1 1 0 0;0 2 1 0 0;0 0 0 0 0;0 0 0 0 0];
fa=fft2(a);
b=[3 2 1 0 0;3 1 1 0 0;2 0 1 0 0;0 0 0 0 0;0 0 0 0 0];
fb=fft2(b);
temp=eye(5);
for R=1:5
for C=1:5
temp(R,C)=fa(R,C)*fb(R,C);
end
end
res=ifft2(temp);
卷积 
等于 
没错，看红色框框，和我们用笔算的3x3的结果一样；
原因是，如果对于一个有限宽度为A的信号使用宽度为B的窗进行卷积，那么结果将会是一个宽度为A+B-1的信号，但如果信号宽度为A但周期也是A即信号尾部和下一个周期的头是紧挨着的话，那么进行卷积的时候尾部的信号就会和头部的信号相互纠缠，使结果的头部和尾部受到污染，中间部分保持正确，但如果信号比较短，就像上面3x3的话就完全被污染了。
但上面的信号是并不是周期的，而是有限宽度的，问题就出在DFT上，因为DFT把信号扩充了，而且是周期性的扩充，DFT时信号其实是这样的：
所以DFT的两个信号就是周期的而不是我们之前输入的信号，所以就会产生纠缠
避免这个问题很简单，对两个矩阵都用用0进行填充，填充成5x5后，得到正确的卷积结果。
滤波器的特性，滤波器是对傅里叶谱的实部和虚部进行等比增强或抑制，也就是傅里叶谱经过滤波器后，相位是没有变化的，主要原因之前提到过，就是一个谱的相角决定了图片的主要结构特征，所以，如果改变相角，那图片将产生巨大问题，所以使用的滤波器都是0相移的，但如果使用非0相移可以得到其他目的的结果，比如让一幅图面目全非。
对于理想滤波器，我们必须考虑其扩充问题，原因是，离线理想滤波器是频域的，其扩充应该在空域进行，而理想滤波器的IDFT在空域不是有限的而是一个无限的sinc函数，进行阶段后再填充，返回频率域后填充后的滤波器将产生严重的振铃现象，这里我们的办法是只对图像进行填充，填充到原图的4倍大小，然后使用填充后大小的滤波器，来减轻缠绕问题，同时时使问题简单化，但只是缓解缠绕问题，振铃问题依然存在。
上图a为频率域一维理想滤波器，b为IDFT后的空间波形，c是进行填充，d是填充后DFT的结果
下面我来解释下振铃现象，来个官方的解释：
振铃效应（Ringingeffect）是影响复原图像质量的众多因素之一，其典型表现是在图像灰度剧烈变化的邻域出现类吉布斯（Gibbs）分布–(满足给定约束条件且熵最大的分布)的振荡。在图像盲复原中，振铃效应是一个不可忽视的问题，其严重降低了复原图像的质量，并且使得难于对复原图像进行后续处理。振铃效应是由于在图像复原中选取了不适当的图像模型造成的；在图像盲复原中如果点扩散函数选择不准确也是引起复原结果产生振铃效应的另一个原因，特别是选用的点扩散函数尺寸大于真实点扩散函数尺寸时，振铃现象更为明显；振铃效应产生的直接原因是图像退化过程中信息量的丢失，尤其是高频信息的丢失。
官方解释可能不好理解，在频域理想滤波器情况下，其实就是sinc函数产生的，只考虑上图的a和b，a是频域理想滤波器，b就是空域与之对应的滤波器，图像如果和b卷积，结果必然产生类似波浪一样的噪声，图像边缘将会产生抖动，边缘变宽，这就是振铃效果，我们来观察不同截止频率的理想滤波器以及其IDFT效果：
可以看出越宽（截止频率大）的滤波器频域振铃越步明显，相反，越窄的滤波器，空域振铃效果越大。
频域滤波过程
频率域滤波的过程：因为我们要得到的是图像，输入是图像，输出也是图像，这个过程叫图像处理，如果输出的时频谱或其他的，我们叫图像分析，滤波器的主要参数是类型和截止频率。
下面来看一下算法的完整计算过程：
原图：
空域填充：
用0对图片进行填充，得到4倍原图的新图像
频谱中心化
这是为了配合滤波器，因为滤波器一般被设计为中间为低频，四周为高频的模式
DFT：
过度到频域
设计滤波器
FIR滤波器，有限长单位冲激响应滤波器，又称为非递归型滤波器，图像处理中都使用这种滤波器产生一个高斯低通滤波器：
对位相乘：
IDFT：
将滤波结果重建为图像
空域取消填充：
将填充过的图像恢复
代码：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
static void showfilter(double *filter,int width,int height){
IplImage *show=cvCreateImage(cvSize(width, height),8,1);
for(int i=0;i<width;i++)
for(int j=0;j<height;j++){
cvSetReal2D(show, j, i, filter[i*width+j]*255.0);
}
cvNamedWindow("Filter", 1);
cvShowImage("Filter", show);
cvWaitKey(0);
cvSaveImage("/Users/Tony/DIPImage/step.jpg", show, 0);
cvReleaseImage(&show);
}
void MatrixMulti_R_C(double *src1,Complex *src2,Complex *dst,int size){//m(1,1)=a(1,1)*b(1,1);
for(int i=0;i<size;i++){
dst[i].real=src2[i].real*src1[i];
dst[i].imagin=src2[i].imagin*src1[i];
}
}
int ChangtoPower2(int size){
size--;//避免为2的幂的size会被扩大
int i=0;
while ((size/=2)>0) {
i++;
}
return 2<<i;
}
//将图像伸缩到2的幂次大小，并填充
void ResizeMatrix4FFT(IplImage *src,IplImage **dst){
int width=src->width;
int height=src->height;
int re_width=ChangtoPower2(width);
int re_height=ChangtoPower2(height);
IplImage *temp=cvCreateImage(cvSize(re_width, re_height), src->depth, src->nChannels);
cvResize(src, temp, 0);
*dst=cvCreateImage(cvSize(re_width*2, re_height*2), src->depth, src->nChannels);
cvZero(*dst);
for(int i=0;i<re_width;i++)
for(int j=0;j<re_height;j++)
cvSetReal2D(*dst, j, i, cvGetReal2D(temp, j, i));
cvReleaseImage(&temp);
}
//将扩充后的图像还原为左上角的
void CutImage421(IplImage *src,IplImage *dst){
//IplImage *temp=cvCreateImage(cvSize(src->width/2, src->height/2), src->depth, src->nChannels);
int width=dst->width;
int height=dst->height;
for(int i=0;i<width;i++)
for(int j=0;j<height;j++)
cvSetReal2D(dst, j, i, cvGetReal2D(src, j, i));
}
//频域滤波
void FrequencyFiltering(IplImage *src,IplImage *dst,int filter_type,double param1,int param2){
IplImage *temp=NULL;
//调整至2的幂，并用黑色填充，防止周期缠绕
ResizeMatrix4FFT(src, &temp);
int fft_width=temp->width;
int fft_height=temp->height;
//产生滤波器
double *filter=(double *)malloc(sizeof(double)*fft_height*fft_width);
if(filter==NULL){
printf("frequency filter malloc faile");
exit(0);
}
//生成滤波器
switch(filter_type){
case ILPF:
IdealLPFilter(filter, fft_width, fft_height, param1);
break;
case BLPF:
if(param2<0)
param2=2;
ButterworthLPfilter(filter, fft_width, fft_height, param1, param2);
break;
case GLPF:
GaussianLPFilter(filter, fft_width, fft_height, param1);
break;
case IHPF:
IdealHPFilter(filter, fft_width, fft_height, param1);
break;
case BHPF:
if(param2<0)
param2=2;
ButterworthHPfilter(filter, fft_width, fft_height, param1, param2);
break;
case GHPF:
GaussianHPFilter(filter, fft_width, fft_height, param1);
break;
}
//FFT
Complex *temp_complex=(Complex*)malloc(sizeof(Complex)*fft_height*fft_width);//fft结果
if(temp_complex==NULL){
exit(0);
}
ImageFFT(temp, temp_complex);
//相乘
MatrixMulti_R_C(filter,temp_complex,temp_complex,fft_width*fft_height);
//IFFT
ImageIFFT(temp_complex, temp, temp->width, temp->height);
//还原图像
IplImage *result2=cvCreateImage(cvSize(temp->width/2, temp->height/2), temp->depth, temp->nChannels);
CutImage421(temp, result2);
cvResize(result2, dst, 0);
free(filter);
free(temp_complex);
cvReleaseImage(&temp);
cvReleaseImage(&result2);
}
空域滤波与频域滤波
空域和频域的桥梁是傅里叶变换，而纽带是卷积定理，对于频域的特性，我们将其看做一个实验室，进行试验并产生小的滤波模板，将小的滤波模板对空域图像进行循环卷积，得到我们想要的结果，空域小模板多半是奇数大小的（3x3，5x5）的，其计算复杂度低，过程简单，用时较短，这就是之前提到的那个面试的问题的答案，如果我当时把这个过程给面试官讲清楚，估计我现在就是一名图像工作者了。
总结
这就是频域滤波的基础，后面的事情就是设计滤波器，满足不同的问题，有了基础，设计滤波器就要自己发挥了。