概述
（MATLAB/C/Python）快速中值滤波
- 一、中值滤波
- 二、快速中值滤波
- 介绍
- 原理
- 优化
- 三、代码
- MATLAB
- C
- Python
- 四、测试
- 其他
by HPC_ZY
最近一个项目中需要用到中值滤波（不能调库），但是核半径相当大，用传统的方法运行速度极慢。因此查阅文献找到快速中值滤波的方法。写成三种语言并分享给大家。
一、中值滤波
简单说就是，就是对某点邻域内所有像素进行排序，取序数在中间的值替代原始值。
这样做对脉冲噪声有良好的滤除作用，特别是在滤除噪声的同时，能够保护信号的边缘，使之不被模糊。
实现方法
step 1：通过从图像中的某个采样窗口取出奇数个像素进行排序
step 2：用排序后的中值取代要处理的像素即可
代码见后文。
注：本文主要讲图像处理，所有使用 “像素 ” 替代 “ 数据”。
二、快速中值滤波
介绍
中值滤波主要耗时就在于对窗口（邻域）内像素进行排序，因此针对如何获取中值进行优化。
最初——经典的中值滤波对窗口内所有像素进行排序，所以排序算法的选择就很重要。
接着——但后来有人意识到 “我们只需要中间值”，没必要将整个序列进行排序，因此出现了一些快速搜索数组中间值的方法。
然后——由于图像的像素值为0~255的整数，1979年有人就提出了用累积直方图寻找中值的方法。
最后——人们又在这基础上不断改进。
本文主要讲基于直方图的快速中值滤波方法
原理
- 实现方法
step 1：统计 N ∗ N N*N N∗N窗口内像素
step 2：计算累积直方图
step 3：找到累积直方图值刚好超过 N ∗ N 2 frac{N*N}{2} 2N∗N时对应的灰度值。
通过一个简单的例子更直观的说明。
这里假设一个只有8个灰度级的图像，某个 5 ∗ 5 5*5 5∗5的窗口内值如下，
[ 0 4 5 5 5 3 4 7 6 5 6 7 7 7 6 3 2 1 1 0 0 1 1 0 0 ] begin{bmatrix} 0 &4 &5&5&5 \ 3 & 4&7&6&5\6&7&7&7&6\3&2&1&1&0\0&1&1&0&0 end{bmatrix} ⎣⎢⎢⎢⎢⎡0363044721577115671055600⎦⎥⎥⎥⎥⎤
其统计直方图如下表。对于5*5的窗口，中值的序数为13，所以参照累积直方图能快速中值为4。
这里我们可以发现，其实并不需要算出整个累积直方图，只要某一次计算超过半数，就可以停止了。
| 灰度 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| 直方图 | 5 | 4 | 1 | 2 | 2 | 4 | 3 | 4 |
| 累积直方图 | 5 | 9 | 10 | 12 | 14 | 18 | 21 | 25 |
- 分析
直方图方法中，对于图像灰度级为 L L L、窗口尺寸为 N ∗ N N*N N∗N，我们只需要进行以下操作：
1）进行 N ∗ N N*N N∗N次统计（获得直方图）
2）最多进行 L − 1 L-1 L−1次累加+判断（计算累积直方图的同时判断是否超过半数）
通常图像的 L L L都为256，我们可以推断对于 3 ∗ 3 3*3 3∗3的窗口，传统排序方法速度更快；但当N越来越大，直方图法的优势就越来越大。
代码见后文
优化
进一步思考可以发现，我们在平移窗口时，新窗口里有
N
∗
(
N
−
1
)
N*(N-1)
N∗(N−1)的像素都是上一窗口里有的（如下图黄色），如果重新统计就浪费时间。
所以有人提出不必重新计算直方图，只用更新直方图——移除（如下图红色）+加入（如下图绿色）。
代码见后文
三、代码
本文首先用MATLAB讲解演示：传统排序中值法，快速直方图中值法，改进快速直方图法
然后直接分享C、Python的改进快速直方图法。
MATLAB
- 经典方法
为了方便转为C语言，没有使用MATLAB的函数
% 经典中值滤波（不处理边缘）（C格式）
function out = medianfilterC(in,r)
% 参数与初始化
[M,N] = size(in);
out = zeros(M,N,'uint8');
L = 2*r+1;
cidx = L*L/2+0.5; % 如果使用sort(),中值位置在这里
cache = zeros(L*L,1,'uint8');
% 遍历图片
for i = 1+r:M-r
for j = 1+r:N-r
% 获取邻域
for m = -r:r
for n = -r:r
cache((m+r)*L+n+r+1) = in(i+m,j+n);
end
end
% 选择排序(找中值，所以没必要全排序, 且没必要保存过程的值）
for p = 1:cidx
minval = cache(p);
idx = p;
for q = p+1:L*L
if cache(q)<minval
minval = cache(q);
idx = q;
end
end
cache(idx) = cache(p);
end
% 中值
out(i,j) = minval;
end
end
end
- 直方图方法
% 直方图法（不处理边缘）（C格式）
function out = medianfilterCHist(in,r)
% 参数与初始化
[M,N] = size(in);
out = zeros(M,N,'uint8');
L = 2*r+1;
cidx = L*L/2+0.5; % 如果使用sort(),中值位置在这里
% 遍历图片
for i = 1+r:M-r
for j = 1+r:N-r
% 获取邻域直方图
hlist = zeros(256,1,'uint32');
for m = -r:r
for n = -r:r
tmp = in(i+m,j+n)+1;
hlist(tmp) = hlist(tmp)+1;
end
end
% 累积直方图求中值
hsum = 0;
for k = 1:256
hsum = hsum+hlist(k);
if hsum>=cidx
out(i,j) = k-1;
break
end
end
end
end
end
- 改进（更新）直方图方法
% 改进直方图法（不处理边缘）（C格式）
function out = medianfilterCHistUpdata(in,r)
% 参数与初始化
[M,N] = size(in);
out = zeros(M,N,'uint8');
L = 2*r+1;
cidx = L*L/2+0.5; % 如果使用sort(),中值位置在这里
% 遍历图片
for i = 1+r:M-r
%% 计算第一列的
hlist = zeros(256,1,'uint32');
% 获取邻域直方图
for m = -r:r
for n = -r:r
tmp = in(i+m,1+r+n)+1;
hlist(tmp) = hlist(tmp)+1;
end
end
% 累积直方图求中值
hsum = 0;
for k = 1:256
hsum = hsum+hlist(k);
if hsum>=cidx
out(i,1+r) = k-1;
break
end
end
%% 计算后续的
for j = 2+r:N-r
% 更新直方图
for m = -r:r
% 剔除
tmp = in(i+m,j-r-1)+1;
hlist(tmp+1) = hlist(tmp+1)-1;
% 加入
tmp = in(i+m,j+r)+1;
hlist(tmp+1) = hlist(tmp+1)+1;
end
% 累积直方图求中值
hsum = 0;
for k = 1:256
hsum = hsum+hlist(k);
if hsum>=cidx
out(i,j) = k-1;
break
end
end
end
end
end
- 自带
当然了，MATLAB有他自带的中值滤波，贼快
out = medfitler2(in, [N,N]);
C
改进（更新）直方图方法
/*** 中值滤波 ***/
void medfilter2(
BYTE* pImgOut, // （out）滤波后图像（0~255）
const BYTE* pImg, // （in）原始灰度图像（0~255）
int imgWidth, // （in）图像宽（pixel）
int imgHeight, // （in）图像高（pixel）
int nR // （in）窗口核半径（pixel）
)
{
int nL = 2 * nR + 1;
int num = nL*nL;
int cidx = num / 2 + 1; // 中值所在的位置
// 初始化统计直方图
int *hist = (int*)malloc(256 * sizeof(int));
// 开始遍历
int tmp = 0;
for (int i = nR; i < imgHeight - nR; i++)
{
memset(hist, 0, 256 * sizeof(int));
// 第一列
for (int m = -nR; m <= nR; m++)
{
for (int n = -nR; n <= nR; n++)
{
tmp = pImg[(i + m)*imgWidth + (nR + n)];
hist[tmp]++;
}
}
int histsum = 0;
for (int k = 0; k < 256; k++)
{
histsum += hist[k];
if (histsum >= cidx)
{
pImgOut[i*imgWidth + nR] = k;
break;
}
}
// 计算后续
for (int j = 1 + nR; j < imgWidth - nR; j++)
{
for (int m = -nR; m <= nR; m++)
{
tmp = pImg[(i + m)*imgWidth + (j - 1 - nR)];
hist[tmp]--;
tmp = pImg[(i + m)*imgWidth + (j + nR)];
hist[tmp]++;
}
int histsum = 0;
for (int k = 0; k < 256; k++)
{
histsum += hist[k];
if (histsum >= cidx)
{
pImgOut[i*imgWidth + j] = k;
break;
}
}
}
}
}
Python
写时候遇到一个小问题，由于本人刚入门不久，不太懂底层的原理。
用python跑for循环的时候运行速度特别慢，如果有大佬知道原因希望留言告诉我，谢谢了。
代码还是放在这里。
def medfilter( img, r ):
# 参数设置
rows = img.shape[0]
cols = img.shape[1]
L = 2*r + 1
num = L*L
cidx = num/2+0.5
out = np.zeros([rows,cols])
# 循环求解
for i in range(r, rows-r):
# 计算第一列
hist = [0]*256
# 获取直方图
tmp = img[i-r:i+r+1,0:L].flatten()
for n in range(0,num):
hist[tmp[n]] += 1
hist[tmp[n]] += 1
# 累积直方图
histsum = 0
for k in range(0,256):
histsum += hist[k]
if histsum>=cidx:
out[i,r] = k
break
# 后续计算
for j in range(1+r,cols-r):
for m in range(-r,r+1):
tmp = img[i+m, j-1-r]
hist[tmp] -= 1
tmp = img[i+m, j+r]
hist[tmp] += 1
histsum = 0
for k in range(0,256):
histsum += hist[k]
if histsum>=cidx:
out[i,j] = k
break
return out;
四、测试
这里只展示matlab测试结果，
clear; close all; clc
% 制造含噪图像
im = imread('test.jpg');
gray = rgb2gray(im);
in = imnoise(gray,'salt & pepper',0.02); % 椒盐噪声
r = 4;
% MATLAB自带函数
tic
out0 = medfilt2(in,[2*r+1,2*r+1]);
toc
% 经典排序方法
tic
out2 = medianfilterC(in,r);
toc
% 直方图中值法
tic
out3 = medianfilterCHist(in,r);
toc
% 改进直方图中值法
tic
out4 = medianfilterCHistUpdata(in,r);
toc
结果如下图，由于我写的没有处理边缘的像素，所以结果有黑边。这里使用半径为4，如果更大效果就非常明显了。
其他
- 参考文献
《A Fast Two-Dimensional Median Filtering Algorithm》
《Median Filter in Constant Time》
《Fast Median Filtering by Use of Fast Localization of Median Value》 - 由于本人刚入门Python不久，不太懂底层的原理。用python跑for循环的时候运行速度特别慢，如果有大佬知道原因希望留言告诉我，谢谢了。
最后
以上就是糟糕心情为你收集整理的（MATLAB/C/Python）快速中值滤波一、中值滤波二、快速中值滤波三、代码四、测试其他的全部内容，希望文章能够帮你解决（MATLAB/C/Python）快速中值滤波一、中值滤波二、快速中值滤波三、代码四、测试其他所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复