本节主要介绍了FPGA图像处理中常用的图像降噪算法，包括均值滤波算法、中值滤波算法、高斯滤波算法、双边滤波算法，以及它的MATLAB与FPGA实现。
图像噪声简介
在图像传感器成像过程中， 光电转换及数模放大时， 不可避免地会产生噪声； 在图像传输过程中，也将二次引入噪声。 假设原始图像的噪声为$I(x,y)$，噪声noise为随机干扰噪声， 则真实图像的噪声可以用如下公式表示：
$$
f(x,y) = I(x,y)+noise
$$所以， 降噪过程中， 如何有效地去除叠加在原始图像中的噪声， 又尽可能地减少对原始图像数据的影响， 尤其是对原始图像的细节及纹理的保留， 非常重要。
噪声有很多种，常见的噪声有椒盐噪声、高斯噪声、伽马噪声、指数噪声等。其中椒盐噪声是图像信号传输、解码等过程产生的或黑或白的噪声点，通常随机分布。而高斯噪声则是指与当前像素服从正态分布的噪声，通常是由于亮度不足或高温引起的图像传感器噪声。
从图中可见， 松鼠测试图添加椒盐噪声后叠加了随机的黑白点， 符合椒盐噪声异常突出的属性， 类似于图像传感器的坏点； 而添加高斯噪声后则是满屏的噪点， 模拟图像传感器因为照度/散热引起的全幅画面的噪声。
2D降噪思维导图
均值滤波算法
所有滤波算法都是通过当前像素周边的像素， 以一定的权重加权计算滤波后的结果。
因此主要涉及两个变量： 窗口内像素的权重值， 以及窗口的大小。 滤波的窗口有3×3、 5×5、 7×7、 11×11等， 窗口尺度越大， 相应的计算量也越大， 效果也越明显
均值滤波计算流程：
均值滤波可以简单地表示为在邻域窗3×3内， 所有像素权重相同， 简单加权后求平均值。 那么可想而知， 噪声并没有被去除， 只是被平均了而已。 此外， 均值滤波除了抑制噪声， 也有平滑纹理的效果。 如果窗口很大， 也会产生模糊的效果。
1.均值滤波的MATLAB实现
matlab代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21% 灰度图像均值滤波算法实现
% IMG为输入的灰度图像
% n为滤波的窗口大小，为奇数
function Q=avg_filter(IMG,n)
% IMG = rgb2gray(imread('../../0_images/Scart.jpg')); % 读取jpg图像
% n=3;
[h,w] = size(IMG);
win = zeros(n,n);
Q = zeros(h,w);
for i=1 : h
for j=1:w
if(i<(n-1)/2+1 || i>h-(n-1)/2 || j<(n-1)/2+1 || j>w-(n-1)/2)
Q(i,j) = IMG(i,j); %边缘像素取原值
else
win = IMG(i-(n-1)/2:i+(n-1)/2, j-(n-1)/2:j+(n-1)/2);
Q(i,j)=sum(sum(win)) / (n*n); %n*n窗口的矩阵，求和再均值
end
end
end
Q=uint8(Q);
均值滤波对椒盐噪声的处理并不是那么理想的。 均值滤波算法是较简单的滤波算法， 计算量较小， 其成效显而易见， 对噪声能起到一个平滑的作用。
2.均值滤波的FPGA实现**·
均值滤波FPGA实现的关键在于3*3滑窗的缓存以及图像边缘的处理
为3×3窗口生成的模块， 为了生成以目标像素为中心的3×3窗口， 需要缓存3行像素， 但在设计时只需要缓存两行像素， 第3行像素实时输入即可
计算均值时在FPGA中实现除法器会占用较多的逻辑资源。由于计算公式中分母时一个常数，因此可以将除法运算转化为乘法运算和移位操作，从而减少逻辑资源的消耗，计算公式如下：
$$
P=uint8(\frac{sum}{9})=round[\frac{sum\times round(\frac{2^{15}}{9})}{2^{15}}]=round(sum\times 3641 >> 15)
$$根据上述计算公式，可将均值滤波运算分解为以下几个步骤
计算$3\times 3$窗口中每行3个像素的累加和
$$
data_sum1=matrix_p11+matrix_p12+matrix_p13\
data_sum2=matrix_p21+matrix_p22+matrix_p23\
data_sum3=matrix_p31+matrix_p32+matrix_p33\
$$计算$3\times 3$窗口中所有9个像素的累加和
$$
data_sum=data_sum1+data_sum2+data_sum3
$$计算$data_mult=data_sum\times 3641$，其中$data_mult[22:15]$为整数部分，$data_mult[14:0]$为小数部分
对$data_mult$进行四舍五入计算，即$avg_data=data_mult[22:15]+data_mult[14]$，即得到均值滤波的结果
判断$3\times 3$窗口的中心像素是否位于图像边界，如果位于图像边界，则直接将中心像素作为均值滤波的结果输出，否则，将$avg_data$作为均值滤波的结果输出。
Verilog源代码（没有全粘，只粘贴了我认为的关键代码）：
Matrix_Generate_3X3_8Bit.v
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
233
234
235
236
237
238
239
240
241
242module Matrix_Generate_3X3_8Bit
#(
parameter [10:0] IMG_HDISP = 11'd640, // 640*480
parameter [10:0] IMG_VDISP = 11'd480,
parameter [10:0] DELAY_NUM = 11'd10 // Interval period from the penultimate row to the last row
)
(
// global clock & reset
input wire clk ,
input wire rst_n ,
// Image data prepared to be processed
input wire per_img_vsync , // Prepared Image data vsync valid signal
input wire per_img_href , // Prepared Image data href vaild signal
input wire [7:0] per_img_gray , // Prepared Image brightness input
// Image data has been processed
output wire matrix_img_vsync , // processed Image data vsync valid signal
output wire matrix_img_href , // processed Image data href vaild signal
output wire matrix_top_edge_flag , // processed Image top edge
output wire matrix_bottom_edge_flag , // processed Image bottom edge
output wire matrix_left_edge_flag , // processed Image left edge
output wire matrix_right_edge_flag , // processed Image right edge
output reg [7:0] matrix_p11 , // 3X3 Matrix output
output reg [7:0] matrix_p12 ,
output reg [7:0] matrix_p13 ,
output reg [7:0] matrix_p21 ,
output reg [7:0] matrix_p22 ,
output reg [7:0] matrix_p23 ,
output reg [7:0] matrix_p31 ,
output reg [7:0] matrix_p32 ,
output reg [7:0] matrix_p33
);
//----------------------------------------------------------------------
// href & vsync counter
reg [10:0] hcnt;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
hcnt <= 11'b0;
else
begin
if(per_img_href == 1'b1)
hcnt <= hcnt + 1'b1;
else
hcnt <= 11'b0;
end
end
reg per_img_href_dly;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
per_img_href_dly <= 1'b0;
else
per_img_href_dly <= per_img_href;
end
wire img_href_neg = ~per_img_href & per_img_href_dly; // falling edge of per_img_href
reg [10:0] vcnt;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
vcnt <= 11'b0;
else
begin
if(per_img_vsync == 1'b0)
vcnt <= 11'b0;
else if(img_href_neg == 1'b1)
vcnt <= vcnt + 1'b1;
else
vcnt <= vcnt;
end
end
//----------------------------------------------------------------------
// two fifo for raw data buffer
reg [10:0] extend_last_row_cnt;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
extend_last_row_cnt <= 11'b0;
else
begin
if((per_img_href == 1'b1)&&(vcnt == IMG_VDISP - 1'b1)&&(hcnt == IMG_HDISP - 1'b1))
extend_last_row_cnt <= 11'd1;
else if((extend_last_row_cnt > 11'b0)&&(extend_last_row_cnt < DELAY_NUM + IMG_HDISP))
extend_last_row_cnt <= extend_last_row_cnt + 1'b1;
else
extend_last_row_cnt <= 11'b0;
end
end
wire extend_last_row_en = (extend_last_row_cnt > DELAY_NUM) ? 1'b1 : 1'b0;
wire fifo1_wenb;
wire [7:0] fifo1_wdata;
wire fifo1_renb;
wire [7:0] fifo1_rdata;
wire fifo2_wenb;
wire [7:0] fifo2_wdata;
wire fifo2_renb;
wire [7:0] fifo2_rdata;
assign fifo1_wenb = per_img_href;
assign fifo1_wdata = per_img_gray;
assign fifo1_renb = per_img_href & (vcnt > 11'b0) | extend_last_row_en;
assign fifo2_wenb = per_img_href & (vcnt > 11'b0);
assign fifo2_wdata = fifo1_rdata;
assign fifo2_renb = per_img_href & (vcnt > 11'b1) | extend_last_row_en;
sync_fifo
#(
.C_FIFO_WIDTH (8 ),
.C_FIFO_DEPTH (1024 )
)
u1_sync_fifo
(
.rst (~rst_n ),
.clk (clk ),
.wr_en (fifo1_wenb ),
.din (fifo1_wdata),
.full ( ),
.rd_en (fifo1_renb ),
.dout (fifo1_rdata),
.empty ( ),
.data_count ( )
);
sync_fifo
#(
.C_FIFO_WIDTH (8 ),
.C_FIFO_DEPTH (1024 )
)
u2_sync_fifo
(
.rst (~rst_n ),
.clk (clk ),
.wr_en (fifo2_wenb ),
.din (fifo2_wdata),
.full ( ),
.rd_en (fifo2_renb ),
.dout (fifo2_rdata),
.empty ( ),
.data_count ( )
);
//----------------------------------------------------------------------
// Read data from fifo
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
begin
{matrix_p11, matrix_p12, matrix_p13} <= 24'h0;
{matrix_p21, matrix_p22, matrix_p23} <= 24'h0;
{matrix_p31, matrix_p32, matrix_p33} <= 24'h0;
end
else
begin
{matrix_p11, matrix_p12, matrix_p13} <= {matrix_p12, matrix_p13, fifo2_rdata}; // 1st row input
{matrix_p21, matrix_p22, matrix_p23} <= {matrix_p22, matrix_p23, fifo1_rdata}; // 2nd row input
{matrix_p31, matrix_p32, matrix_p33} <= {matrix_p32, matrix_p33, per_img_gray}; // 3rd row input
end
end
reg extend_last_row_en_dly;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
extend_last_row_en_dly <= 1'b0;
else
extend_last_row_en_dly <= extend_last_row_en;
end
reg [1:0] vsync;
reg [1:0] href;
reg [1:0] top_edge_flag;
reg [1:0] bottom_edge_flag;
reg [1:0] left_edge_flag;
reg [1:0] right_edge_flag;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
vsync <= 2'b0;
else
begin
if((per_img_href == 1'b1)&&(vcnt == 11'd1)&&(hcnt == 11'b0))
vsync[0] <= 1'b1;
else if((extend_last_row_en == 1'b0)&&(extend_last_row_en_dly == 1'b1))
vsync[0] <= 1'b0;
else
vsync[0] <= vsync[0];
vsync[1] <= vsync[0];
end
end
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
begin
href <= 2'b0;
top_edge_flag <= 2'b0;
bottom_edge_flag <= 2'b0;
left_edge_flag <= 2'b0;
right_edge_flag <= 2'b0;
end
else
begin
href[0] <= per_img_href & (vcnt > 11'b0) | extend_last_row_en;
href[1] <= href[0];
top_edge_flag[0] <= per_img_href & (vcnt == 11'd1);
top_edge_flag[1] <= top_edge_flag[0];
bottom_edge_flag[0] <= extend_last_row_en;
bottom_edge_flag[1] <= bottom_edge_flag[0];
left_edge_flag[0] <= per_img_href & (vcnt > 11'b0) & (hcnt == 11'b0) | (extend_last_row_cnt == DELAY_NUM + 1'b1);
left_edge_flag[1] <= left_edge_flag[0];
right_edge_flag[0] <= per_img_href & (vcnt > 11'b0) & (hcnt == IMG_HDISP - 1'b1) | (extend_last_row_cnt == DELAY_NUM + IMG_HDISP);
right_edge_flag[1] <= right_edge_flag[0];
end
end
assign matrix_img_vsync = vsync[1];
assign matrix_img_href = href[1];
assign matrix_top_edge_flag = top_edge_flag[1];
assign matrix_bottom_edge_flag = bottom_edge_flag[1];
assign matrix_left_edge_flag = left_edge_flag[1];
assign matrix_right_edge_flag = right_edge_flag[1];
endmodulemean_filter_proc.v
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
159module mean_filter_proc
#(
parameter [10:0] IMG_HDISP = 11'd640, // 640*480
parameter [10:0] IMG_VDISP = 11'd480
)
(
input wire clk ,
input wire rst_n ,
// Image data prepared to be processed
input wire per_img_vsync , // Prepared Image data vsync valid signal
input wire per_img_href , // Prepared Image data href vaild signal
input wire [7:0] per_img_gray , // Prepared Image brightness input
// Image data has been processed
output reg post_img_vsync , // processed Image data vsync valid signal
output reg post_img_href , // processed Image data href vaild signal
output reg [7:0] post_img_gray // processed Image brightness output
);
//----------------------------------------------------------------------
// Generate 8Bit 3X3 Matrix
wire matrix_img_vsync;
wire matrix_img_href;
wire matrix_top_edge_flag;
wire matrix_bottom_edge_flag;
wire matrix_left_edge_flag;
wire matrix_right_edge_flag;
wire [7:0] matrix_p11;
wire [7:0] matrix_p12;
wire [7:0] matrix_p13;
wire [7:0] matrix_p21;
wire [7:0] matrix_p22;
wire [7:0] matrix_p23;
wire [7:0] matrix_p31;
wire [7:0] matrix_p32;
wire [7:0] matrix_p33;
Matrix_Generate_3X3_8Bit
#(
.IMG_HDISP (IMG_HDISP ),
.IMG_VDISP (IMG_VDISP )
)
u_Matrix_Generate_3X3_8Bit
(
// global clock & reset
.clk (clk ),
.rst_n (rst_n ),
// Image data prepared to be processed
.per_img_vsync (per_img_vsync ), // Prepared Image data vsync valid signal
.per_img_href (per_img_href ), // Prepared Image data href vaild signal
.per_img_gray (per_img_gray ), // Prepared Image brightness input
// Image data has been processed
.matrix_img_vsync (matrix_img_vsync ), // processed Image data vsync valid signal
.matrix_img_href (matrix_img_href ), // processed Image data href vaild signal
.matrix_top_edge_flag (matrix_top_edge_flag ), // processed Image top edge
.matrix_bottom_edge_flag(matrix_bottom_edge_flag), // processed Image bottom edge
.matrix_left_edge_flag (matrix_left_edge_flag ), // processed Image left edge
.matrix_right_edge_flag (matrix_right_edge_flag ), // processed Image right edge
.matrix_p11 (matrix_p11 ), // 3X3 Matrix output
.matrix_p12 (matrix_p12 ),
.matrix_p13 (matrix_p13 ),
.matrix_p21 (matrix_p21 ),
.matrix_p22 (matrix_p22 ),
.matrix_p23 (matrix_p23 ),
.matrix_p31 (matrix_p31 ),
.matrix_p32 (matrix_p32 ),
.matrix_p33 (matrix_p33 )
);
//----------------------------------------------------------------------
// calc sum of [p11,p12,p13;p21,p22,p23;p31,p32,p33]
reg [ 9:0] data_sum1;
reg [ 9:0] data_sum2;
reg [ 9:0] data_sum3;
reg [11:0] data_sum;
always @(posedge clk)
begin
data_sum1 <= matrix_p11 + matrix_p12 + matrix_p13;
data_sum2 <= matrix_p21 + matrix_p22 + matrix_p23;
data_sum3 <= matrix_p31 + matrix_p32 + matrix_p33;
data_sum <= data_sum1 + data_sum2 + data_sum3;
end
//----------------------------------------------------------------------
// avg_data = round(data_sum/9.0) -> avg_data = round(data_sum*3641 >> 15)
reg [22:0] data_mult;
always @(posedge clk)
begin
data_mult <= data_sum * 12'd3641;
end
reg [7:0] avg_data;
always @(posedge clk)
begin
avg_data <= data_mult[22:15] + data_mult[14];
end
//----------------------------------------------------------------------
// lag 4 clocks signal sync
reg [3:0] matrix_img_vsync_r1;
reg [3:0] matrix_img_href_r1;
reg [3:0] matrix_edge_flag_r1;
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
begin
matrix_img_vsync_r1 <= 4'b0;
matrix_img_href_r1 <= 4'b0;
matrix_edge_flag_r1 <= 4'b0;
end
else
begin
matrix_img_vsync_r1 <= {matrix_img_vsync_r1[2:0],matrix_img_vsync};
matrix_img_href_r1 <= {matrix_img_href_r1[2:0],matrix_img_href};
matrix_edge_flag_r1 <= {matrix_edge_flag_r1[2:0],matrix_top_edge_flag | matrix_bottom_edge_flag | matrix_left_edge_flag | matrix_right_edge_flag};
end
end
reg [7:0] matrix_p22_r1 [0:3];
always @(posedge clk)
begin
matrix_p22_r1[0] <= matrix_p22;
matrix_p22_r1[1] <= matrix_p22_r1[0];
matrix_p22_r1[2] <= matrix_p22_r1[1];
matrix_p22_r1[3] <= matrix_p22_r1[2];
end
//----------------------------------------------------------------------
// result output
always @(posedge clk)
begin
if(matrix_edge_flag_r1[3] == 1'b1)
post_img_gray <= matrix_p22_r1[3];
else
post_img_gray <= avg_data;
end
always @(posedge clk or negedge rst_n)
begin
if(!rst_n)
begin
post_img_vsync <= 1'b0;
post_img_href <= 1'b0;
end
else
begin
post_img_vsync <= matrix_img_vsync_r1[3];
post_img_href <= matrix_img_href_r1[3];
end
end
endmodule仿真结果：
中值滤波算法
顾名思义， 中值滤波算法就是取滤波窗口内的中间值进行计算的算法。
对椒盐噪声， 均值滤波并不能很好地去除它， 噪声只是被平均化了； 但中值滤波很好地去除了异常的椒盐噪声。 不过仔细观察原图与中值滤波后的图， 处理后的图在边缘纹理上， 也有一定程度的丢失。
如何快速实现窗口内像素的中间值获取？并行3步中值法
- 分别对每行3个像素进行两两比较， 得到最大值Max、 中间值Mid、 最小值Min
- 求3个最大值的最小值Max_Min、 3个中间值的中间值Mid_Mid， 以及3个最小值的最大值Min_Max
- 求Max_Min、 Mid_Mid、 Min_Max的中间值， 即为最终结果
1.中值滤波的MATLAB实现
matlab实现代码：
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% 灰度图像中值滤波算法实现
% IMG为输入的灰度图像
% n为滤波的窗口大小，为奇数
function Q=med_filter(IMG,n) %目前n只能等于3
% IMG = rgb2gray(imread('../../0_images/Scart.jpg')); % 读取jpg图像
% n=3;
[h,w] = size(IMG);
win = zeros(n,n);
Q = zeros(h,w);
for i=1 : h
for j=1:w
if(i<(n-1)/2+1 || i>h-(n-1)/2 || j<(n-1)/2+1 || j>w-(n-1)/2)
Q(i,j) = IMG(i,j); %边缘像素取原值
else
win = IMG(i-(n-1)/2:i+(n-1)/2, j-(n-1)/2:j+(n-1)/2); %n*n窗口的矩阵
% Q(i,j)=median(median(win)); %求中值
max1 = max(win(1,1:3)); mid1 = median(win(1,1:3)); min1 = min(win(1,1:3));
max2 = max(win(2,1:3)); mid2 = median(win(2,1:3)); min2 = min(win(2,1:3));
max3 = max(win(3,1:3)); mid3 = median(win(3,1:3)); min3 = min(win(3,1:3));
max_min = min([max1, max2, max3]);
mid_mid = median([mid1, mid2, mid3]);
min_max = max([min1, min2, min3]);
Q(i,j) = median([max_min, mid_mid, min_max]);
end
end
end
Q=uint8(Q);
% subplot(121);imshow(IMG1);title('【1】原图');
% subplot(122);imshow(IMG4);title('【2】3*3中值滤波');
2.中值滤波的FPGA实现*
中值滤波中生成$3\times 3$窗口方案跟均值滤波一致，这里主要介绍FPGA如何实现并行3步中值法
中值滤波计算分解为以下3步（即3clk delay）：
分别对3×3窗口中每行3个像素进行比较， 得到3个最大值row1_max_data、 row2_max_data、row3_max_data， 3个中间值row1_med_data、row2_med_data、 row3_med_data， 3个最小值row1_min_data、 row2_min_data、 row3_min_data。
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
29always @(posedge clk)
begin
if((matrix_p11 <= matrix_p12)&&(matrix_p11 <= matrix_p13))
row1_min_data <= matrix_p11;
else if((matrix_p12 <= matrix_p11)&&(matrix_p12 <= matrix_p13))
row1_min_data <= matrix_p12;
else
row1_min_data <= matrix_p13;
end
always @(posedge clk)
begin
if((matrix_p11 <= matrix_p12)&&(matrix_p11 >= matrix_p13)||(matrix_p11 >= matrix_p12)&&(matrix_p11 <= matrix_p13))
row1_med_data <= matrix_p11;
else if((matrix_p12 <= matrix_p11)&&(matrix_p12 >= matrix_p13)||(matrix_p12 >= matrix_p11)&&(matrix_p12 <= matrix_p13))
row1_med_data <= matrix_p12;
else
row1_med_data <= matrix_p13;
end
always @(posedge clk)
begin
if((matrix_p11 >= matrix_p12)&&(matrix_p11 >= matrix_p13))
row1_max_data <= matrix_p11;
else if((matrix_p12 >= matrix_p11)&&(matrix_p12 >= matrix_p13))
row1_max_data <= matrix_p12;
else
row1_max_data <= matrix_p13;
end求3个最大值的最小值min_of_max_data、 3个中间值的中间值med_of_med_data， 以及3个最小值的最大值max_of_min_data
求min_of_max_data、 med_of_med_data、 max_of_min_data的中间值pixel_data， 得到中值滤波的结果
仿真结果：
高斯滤波算法
一维、二维高斯分布函数：
$$
\begin{cases}
\begin{aligned}
G(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^2}{2\sigma^2}}(一维高斯分布函数)\
G(x,y)=\frac{1}{\sqrt{2\pi}\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}(二维高斯分布函数)
\end{aligned}
\end{cases}
$$$\sigma$为高斯分布的强度，$\sigma$越大，数据越分散；反之，数据越向中心集中分布。
$\sigma = 3$的高斯滤波， 窗口扩大后， 模糊程度非常大。 可见， 滤波窗口对滤波强度的影响比$\sigma$的大小对滤波强度的影响更大
1.高斯滤波的MATLAB实现
根据二维高斯分布函数生成定点化模板G3：
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
59clear all; close all; clc;
% -------------------------------------------------------------------------
% ------------------------------生成G3模板-------------------------------
% -------------------------------------------------------------------------
% 计算5*5高斯模板
sigma = 3;
G1=zeros(5,5); %5*5高斯模板
for i=-2 : 2
for j=-2 : 2
% G1(i+3,j+3) = exp(-(i.^2 + j.^2)/(2*sigma^2)) / (2*pi*sigma^2);
G1(i+3,j+3) = exp(-(i^2 + j^2)/(2*sigma^2)) ;
end
end
% 归一化5*5高斯模板
temp = sum(sum(G1));
G2 = G1/temp;
% 5*5高斯模板 *1024定点化
G3 = floor(G2*1024);
% -------------------------------------------------------------------------
% ------------------------------调用-------------------------------
% -------------------------------------------------------------------------
clear all; close all; clc;
% -------------------------------------------------------------------------
% Read PC image to Matlab
IMG1 = imread('../../0_images/Scart.jpg');
IMG1 = rgb2gray(IMG1);
h = size(IMG1,1); % 读取图像高度
w = size(IMG1,2); % 读取图像宽度
subplot(131);imshow(IMG1);title('【1】原图');
% -------------------------------------------------------------------------
g = fspecial('gaussian',[5,5],3);
IMG2 = imfilter(IMG1, g, 'replicate');
subplot(132);imshow(IMG2);title('【2】Matlab自带高斯滤波');
% -------------------------------------------------------------------------
G =[32 38 40 38 32; ...
38 45 47 45 38; ...
40 47 50 47 40; ...
38 45 47 45 38; ...
32 38 40 38 32];
IMG3 = zeros(h,w);
n=5;
for i=1 : h
for j=1:w
if(i<(n-1)/2+1 || i>h-(n-1)/2 || j<(n-1)/2+1 || j>w-(n-1)/2)
IMG3(i,j) = IMG1(i,j); %边缘像素取原值
else
IMG3(i,j) = conv2(double(IMG1(i-(n-1)/2:i+(n-1)/2, j-(n-1)/2:j+(n-1)/2)), double(G), 'valid')/1024;
end
end
end
IMG3 = uint8(IMG3);
subplot(133);imshow(IMG3);title('【3】手动编写高斯滤波');
2.高斯滤波的FPGA实现**
为了生成以目标像素为中心的5×5窗口， 需要缓存5行像素， 但在设计时只需要缓存4行像素， 第5行像素实时输入即可。
fifo1_wenb与per_img_href的时序保持一致
fifo1_renb在per_img_href的第二行开始时有效， 并在一帧结束后多读取一行（只有在前一个fifo写完一行数据后，下一个fifo才会写入前一个fifo读出的数据）
fifo2_wenb在per_img_href的第二行开始时有效， 并在一帧结束后多写一行（因为$5\times5$窗口有两层边缘情况，此时需要再竖直往下移一格，这样才能重新计算最后一层边缘的中心值）
fifo2_renb在per_img_href的第三行开始时有效， 并在一帧结束后多读取两行
fifo3_wenb在per_img_href的第三行开始时有效， 并在一帧结束后多写一行
fifo3_renb在per_img_href的第四行开始时有效， 并在一帧结束后多读取两行
fifo4_wenb在per_img_href的第四行开始时有效， 并在一帧结束后多写一行
fifo4_renb在per_img_href的第五行开始时有效， 并在一帧结束后多读取两行
这里以$5\times 5$的滤波窗口为例，高斯滤波计算过程分为以下几个步骤：
将$5\times 5$窗口的像素与$5\times 5$窗口的高斯模板相乘得到$5\times 5$窗口的矩阵$mult_result$。因为高斯模板已经放大了1024倍，所以$mult_result$也放大了1024倍
计算$mult_result$矩阵每行5个数值的累加和
$$
sum_result1=mult_result11+mult_result12+mult_result13+mult_result14+mult_result15\
sum_result2=mult_result21+mult_result22+mult_result23+mult_result24+mult_result25\
sum_result3=mult_result31+mult_result32+mult_result33+mult_result34+mult_result35\
sum_result4=mult_result41+mult_result42+mult_result43+mult_result44+mult_result45\
sum_result5=mult_result51+mult_result52+mult_result53+mult_result54+mult_result55\
$$计算$mult_result$矩阵25个数值的累加和得到$sum_result$，其中$sum_result[17:10]$是整数部分，$sum_result[9:0]$是小数部分
$$
sum_result=sum_result1+sum_result2+sum_result3+sum_result4+sum_result5
$$对$sum_result$进行四舍五入计算，即$pixel_data = sum_result[17:10]+sum_result[9]$，即得到高斯滤波的结果
判断5×5窗口的中心像素是否位于图像边界， 即中心像素是否处于图像上边界2行、 下边界2行、 左边界2列、 右边界2列范围内。 如果中心像素位于图像边界， 则直接将中心像素作为高斯滤波的结果输出； 否则， 将pixel_data作为高斯滤波的结果输出。
仿真结果：
双边滤波算法
双边滤波是一种非线性滤波器，它既可以达到降噪平滑，同时又保持边缘的效果。
双边滤波的权重，不仅考虑了像素的空间距离（如高斯滤波），还考虑了像素范围的辐射差异（如像素与中心像素的相似程度，也是高斯分布的）。结合空间距离与相似度，计算得到最终的权重（空间距离与相似度的高斯分布）
计算方法：首先计算当前像素的高斯权重，然后计算当前像素的相似度权重，最后两者结合计算得到最终的权重
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% 灰度图像双边滤波算法实现
% I为输入的灰度图像
% n为滤波的窗口大小，为奇数
function B=bilateral_filter_gray(I,n,sigma_d, sigma_r)
% ---------------------------------------------------
% 仅供function自测使用
% clear all; close all; clc;
% I = rgb2gray(imread('../../0_images/Scart.jpg')); % 读取jpg图像
% n = 3; sigma_d = 3; sigma_r = 0.8;
dim = size(I); %读取图像高度、宽度
w=floor(n/2); %窗口 [-w, w]
% ---------------------------------------------------
% 计算n*n高斯模板
G1=zeros(n,n); %n*n高斯模板
for i=-w : w
for j=-w : w
G1(i+w+1, j+w+1) = exp(-(i^2 + j^2)/(2*sigma_d^2)) ;
end
end
% 归一化n*n高斯滤波模板
temp = sum(G1(:));
G2 = G1/temp;
% n*n高斯模板 *1023定点化
G3 = floor(G2*1023);
I = double(I);
% ---------------------------------------------------
% 计算n*n双边滤波模板+ 滤波结果
h = waitbar(0,'Speed of bilateral filter process...'); %创建进度条
B = zeros(dim);
for i=1 : dim(1)
for j=1 : dim(2)
if(i<w+1 || i>dim(1)-w || j<w+1 || j>dim(2)-w)
B(i,j) = I(i,j); %边缘像素取原值
else
A = I(i-w:i+w, j-w:j+w);
% H = exp( -(I(i,j)-A).^2/(2*sigma_r^2) ) ;
H = exp( -((A-I(i,j))/255).^2/(2*sigma_r^2)) ;%计算以当前像素为中心的n×n窗口内的相似度权重
F = G3.*H; %将高斯权重与相似度权重相乘， 得到同时考虑当前像素空间距离与相似度的双边滤波权重
F2=F/sum(F(:));
B(i,j) = sum(F2(:) .*A(:));
end
end
waitbar(i/dim(1));
end
close(h); % Close waitbar.
I = uint8(I);
B = uint8(B);
% ---------------------------------------------------
% 仅供function自测使用
% subplot(121);imshow(I);title('【1】原始图像');
% subplot(122);imshow(B);title('【2】双边滤波结果');不同sigma与相同窗口的高斯滤波与双边滤波效果对比：双边滤波对边缘的保护程度更好
不同sigma_d下双边滤波的效果对比：sigma_d对滤波的影响并不大
不同sigma_r下双边滤波的效果对比：sigma_r越大，图像的双边滤波强度越大
不同窗口大小的双边滤波效果对比：n 变大时， 磨皮的效果都要出来了， 可见窗口大小对双边滤波强度的影响最大
综上分析， 对双边滤波强度的影响， 首先是窗口的大小， 其次是sigma_r， 最后才是sigma_d。 所以固定窗口下， sigma_r有较大的权重。 也正是因此，像素相似度是双边滤波权重的一个重要因素。
1.双边滤波算法的MATLAB实现
MATLAB中双边滤波的实现是为FPGA实现服务的，不难看出，相邻两个像素的差值的绝对值，一定属于$[0,255]$，所以可以提前计算好$e^{\frac{i^2}{2\sigma_r^2}}$，这样就可以用查找表的方法来实现像素相似度权重的计算了（妙！）
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% 灰度图像双边滤波算法实现
% I为输入的灰度图像
% n为滤波的窗口大小，为奇数
function B=bilateral_filter_gray_INT(I,n,sigma_d, sigma_r)
% clear all; close all; clc;
% I = rgb2gray(imread('../../0_images/Scart.jpg')); % 读取jpg图像
% n = 3; sigma_d = 3; sigma_r = 0.1;
dim = size(I); %读取图像高度、宽度
w=floor(n/2); %窗口 [-w, w]
% ---------------------------------------------------
% 计算n*n高斯模板
G1=zeros(n,n); %n*n高斯模板
for i=-w : w
for j=-w : w
G1(i+w+1, j+w+1) = exp(-(i^2 + j^2)/(2*sigma_d^2)) ;
end
end
% 归一化n*n高斯滤波模板
temp = sum(G1(:));
G2 = G1/temp;
% n*n高斯模板 *1024定点化
G3 = floor(G2*1024);
% ---------------------------------------------------
% 计算相似度指数*1023结果
% H = zeros(1, 256);
for i=0 : 255
H(i+1) = exp( -(i/255)^2/(2*sigma_r^2));
end
H = floor(H *1023);
I = double(I);
% ---------------------------------------------------
% 计算n*n双边滤波模板+ 滤波结果
h = waitbar(0,'Speed of bilateral filter process...'); %创建进度条
B = zeros(dim);
for i=1 : dim(1)
for j=1 : dim(2)
if(i<w+1 || i>dim(1)-w || j<w+1 || j>dim(2)-w)
B(i,j) = I(i,j); %边缘像素取原值
else
A = I(i-w:i+w, j-w:j+w);
% H = exp( -(I(i,j)-A).^2/(2*sigma_r^2) ) ;
F1 = reshape(H(abs(A-I(i,j))+1), n, n); %计算相似度权重(10bit)
F2 = F1 * G3; %计算双边权重（20bit）
F3=floor(F2*1024/sum(F2(:))); %归一化双边滤波权重（扩大1024）
B(i,j) = sum(F3(:) .*A(:))/1024; %卷积后得到最终累加的结果（缩小1024）
end
end
waitbar(i/dim(1));
end
close(h); % Close waitbar.
I = uint8(I);
B = uint8(B);
% subplot(121);imshow(I);title('【1】原始图像');
% subplot(122);imshow(B);title('【2】双边滤波结果');经常说双边滤波算法是一种磨皮算法，既去掉了斑纹又保持了边缘
1
2
3
4
5
6
7
8
9
10
11
12
13clear all; close all; clc;
% -------------------------------------------------------------------------
% Read PC image to Matlab
IMG1= imread('../../0_images/girl.jpg'); % 读取jpg图像
% -------------------------------------------------------------------------
subplot(121);imshow(IMG1);title('【1】原图');
IMG2(:,:,1) = bilateral_filter_gray(IMG1(:,:,1), 3, 3, 0.1);
IMG2(:,:,2) = bilateral_filter_gray(IMG1(:,:,2), 3, 3, 0.1);
IMG2(:,:,3) = bilateral_filter_gray(IMG1(:,:,3), 3, 3, 0.1);
subplot(122);imshow(IMG2);title('【2】双边滤波3*3, sigma = [3, 0.1]');
2.双边滤波算法的FPGA实现**
相似度查找表在matlab中直接写入Verilog文件用以索引：
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
45clear all; close all; clc;
sigma_d = 3; sigma_r=0.3;
n=3;
w=floor(n/2);
% ---------------------------------------------------
% 计算n*n高斯模板
G1=zeros(n,n);
for i=-w : w
for j=-w : w
G1(i+w+1, j+w+1) = exp(-(i^2 + j^2)/(2*sigma_d^2)) ;
end
end
% 归一化n*n高斯模板
temp = sum(sum(G1));
G2 = G1/temp;
% n*n高斯模板 *1024定点化
G3 = floor(G2*1024);
% ----------------------------------------------------------------------
fp_gray = fopen('.\Similary_LUT.v','w');
fprintf(fp_gray,'//Sigma_r = %f\n', sigma_r);
fprintf(fp_gray,'module Similary_LUT\n');
fprintf(fp_gray,'(\n');
fprintf(fp_gray,' input\t\t[7:0]\tPre_Data,\n');
fprintf(fp_gray,' output\treg\t[9:0]\tPost_Data\n');
fprintf(fp_gray,');\n\n');
fprintf(fp_gray,'always@(*)\n');
fprintf(fp_gray,'begin\n');
fprintf(fp_gray,'\tcase(Pre_Data)\n');
% -----------------------------
% 计算相似度指数*1023结果
Similary_ARRAY = zeros(1,256);
for i = 0 : 255
temp = exp(-(i/255)^2/(2*sigma_r^2));
Similary_ARRAY(i+1) = floor(temp*1023);
fprintf(fp_gray,'\t8''h%s : Post_Data = 10''h%s; \n',dec2hex(i,2), dec2hex(Similary_ARRAY(i+1),3));
end
fprintf(fp_gray,'\tendcase\n');
fprintf(fp_gray,'end\n');
fprintf(fp_gray,'\nendmodule\n');
fclose(fp_gray);双边滤波计算分解为以下步骤
计算当前像素为中心的$3\times 3$窗口内的相似度权重：计算两个数差的绝对值$c=abs(a-b)$，可转化为：$if ,a >b ,c= a-b ; if ,a\le b ,c=b-a$
将$3\times 3$窗口的高斯权重与$3\times 3$窗口的相似度权重相乘，得到同时考虑距离与相似度的双边滤波权重
对$3\times 3$窗口的双边滤波权重进行归一化处理
将$3\times 3$窗口的像素与$3\times 3$窗口的双边滤波权重进行卷积并累加，得到$sum_result$，其中$sum_result[17:10]$为整数部分、$sum_result[9:0]$作为小数部分
对$sum_result$四舍五入，即$sum_result[17:10]+sum_result[9]$作为双边滤波的结果
最后判断$3\times 3$创就的中心像素是否位于图像边界，如果位于图像边界，直接将中心像素作为结果输出
仿真结果：
Reference
- 《基于MATLAB与FPGA的图像处理教程》图书及其配套资料
