（MATLAB/C/Python）快速中值滤波

一、中值滤波
二、快速中值滤波
- 介绍
- 原理
- 优化
三、代码
- MATLAB
- C
- Python
四、测试
其他

by HPC_ZY

最近一个项目中需要用到中值滤波（不能调库），但是核半径相当大，用传统的方法运行速度极慢。因此查阅文献找到快速中值滤波的方法。写成三种语言并分享给大家。

一、中值滤波

简单说就是，就是对某点邻域内所有像素进行排序，取序数在中间的值替代原始值。
这样做对脉冲噪声有良好的滤除作用，特别是在滤除噪声的同时，能够保护信号的边缘，使之不被模糊。

实现方法
step 1：通过从图像中的某个采样窗口取出奇数个像素进行排序
step 2：用排序后的中值取代要处理的像素即可

代码见后文。

注：本文主要讲图像处理，所有使用 “像素 ” 替代 “ 数据”。

二、快速中值滤波

介绍

中值滤波主要耗时就在于对窗口（邻域）内像素进行排序，因此针对如何获取中值进行优化。

最初——经典的中值滤波对窗口内所有像素进行排序，所以排序算法的选择就很重要。
接着——但后来有人意识到 “我们只需要中间值”，没必要将整个序列进行排序，因此出现了一些快速搜索数组中间值的方法。
然后——由于图像的像素值为0~255的整数，1979年有人就提出了用累积直方图寻找中值的方法。
最后——人们又在这基础上不断改进。

本文主要讲基于直方图的快速中值滤波方法

原理

实现方法
step 1：统计N∗NN*NN∗N窗口内像素
step 2：计算累积直方图
step 3：找到累积直方图值刚好超过N∗N2\frac{N*N}{2}2N∗N时对应的灰度值。
通过一个简单的例子更直观的说明。
这里假设一个只有8个灰度级的图像，某个5∗55*55∗5的窗口内值如下，
[0455534765677763211001100]\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

分析
直方图方法中，对于图像灰度级为LLL、窗口尺寸为N∗NN*NN∗N，我们只需要进行以下操作：
1）进行N∗NN*NN∗N次统计（获得直方图）
2）最多进行L−1L-1L−1次累加+判断（计算累积直方图的同时判断是否超过半数）
通常图像的LLL都为256，我们可以推断对于3∗33*33∗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-rfor j = 1+r:N-r% 获取邻域for m = -r:rfor n = -r:rcache((m+r)*L+n+r+1) = in(i+m,j+n);endend% 选择排序(找中值，所以没必要全排序, 且没必要保存过程的值）for p = 1:cidxminval = cache(p);idx = p;for q = p+1:L*Lif cache(q)<minvalminval = cache(q);idx = q;endendcache(idx) = cache(p);end% 中值out(i,j) = minval;end
endend

直方图方法

% 直方图法（不处理边缘）（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-rfor j = 1+r:N-r% 获取邻域直方图hlist = zeros(256,1,'uint32');for m = -r:rfor n = -r:rtmp = in(i+m,j+n)+1;hlist(tmp) = hlist(tmp)+1;endend% 累积直方图求中值hsum = 0;for k = 1:256hsum = hsum+hlist(k);if hsum>=cidxout(i,j) = k-1;breakendendend
endend

改进（更新）直方图方法

% 改进直方图法（不处理边缘）（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:rfor n = -r:rtmp = in(i+m,1+r+n)+1;hlist(tmp) = hlist(tmp)+1;endend% 累积直方图求中值hsum = 0;for k = 1:256hsum = hsum+hlist(k);if hsum>=cidxout(i,1+r) = k-1;breakendend%% 计算后续的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:256hsum = hsum+hlist(k);if hsum>=cidxout(i,j) = k-1;breakendendend
endend

自带
当然了，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 + 1num = L*Lcidx = num/2+0.5out = 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]] += 1hist[tmp[n]] += 1# 累积直方图histsum = 0for k in range(0,256):histsum += hist[k]if histsum>=cidx:out[i,r] = kbreak# 后续计算for j in range(1+r,cols-r):for m in range(-r,r+1):tmp = img[i+m, j-1-r]hist[tmp] -= 1tmp = img[i+m, j+r]hist[tmp] += 1histsum = 0for k in range(0,256):histsum += hist[k]if histsum>=cidx:out[i,j] = kbreak                  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）快速中值滤波相关推荐

快速中值滤波——Python实现
原理中值滤波是空域中常用的一种滤波方式,是一种非线性的滤波.它的原理就是将窗口像素排序,取中值,然后移动窗口,不断重复取中值的过程. 但是,可以发现,每次移动窗口,都需要对像素点进行排序,从而选取中 ...
matlab实现中值滤波去除基线漂移,快速中值滤波在滤除心电信号基线漂移中的应用...
[摘要]文中给出了一种非线性的滤除心电信号基线漂移的滤波方法,把基于排序统计理论的快速中值滤波方法应用于处理心电信号,通过多次对心电信号中选择的窗口数据进行排序,然后取中值的方法来达到滤波的效果.试验 ...
直方图实现快速中值滤波opencv
中值滤波能够有效去除图像中的异常点,具有去除图像噪声的作用.传统中值滤波的算法一般都是在图像中建立窗口,然后对窗口内的所有像素值进行排序,选择排序后的中间值作为窗口中心像素滤波后的值.由于这个做法在每 ...
快速中值滤波在心电图ECG中的应用
1.算法介绍和实现首先来搞明白,什么是快速中值滤波? 快速中值滤波非常简单,就是用过去连续N个数据,再对这N个数据进行排序,取排序后的中间那个数据,做为当前的输出,N即为窗口的长度. 算法实现: 1 ...
OpenCV图像处理专栏九 | 基于直方图的快速中值滤波算法
转载自:https://zhuanlan.zhihu.com/p/98092747 侵删前言这是OpenCV图像处理专栏的第9篇文章,主要介绍一个基于直方图的快速中值滤波算法,希望对大家有帮助. ...
【老生谈算法】matlab实现车牌识别中值滤波算法——车牌识别中值滤波算法
基于Matlab的车牌识别中值滤波算法的研究与实现 1.原文下载: 本算法原文如下,有需要的朋友可以点击进行下载序号原文(点击下载) 本项目原文 [老生谈算法]基于Matlab的车牌识别中值滤波算 ...
matlab设计自适应中值滤波,matlab课程设计(自适应中值滤波).doc
matlab课程设计(自适应中值滤波).doc 10信息工程系课程设计报告课程MATLAB课程设计专业通信工程班级2级本科二班学生姓名1景学号114学生姓名2学号1414学生姓名3王学号6学生姓名4学 ...
python实现中值滤波_Python实现中值滤波去噪方式
中值滤波器去噪: 中值滤波的主要原理是将数字图像中的某点用该点的邻域中各个像素值的中值所来代替,这样就能让目标像素周围能够更好的接近真实值,比如一张白纸上有一个黑点时,黑点的像素值比较大,经过中值滤波 ...
高效快速中值滤波算法c语言,快速中值滤波及c语言实现.docx
. .. 快速中值滤波及c语言实现学生姓名: 刘勇学号: 6100410218 专业班级: 数媒101 [摘要]本文讨论了用c语言在微机上实现中值滤波及快速算法,在程序设计的过程中充分考虑到程 ...

（MATLAB/C/Python）快速中值滤波

（MATLAB/C/Python）快速中值滤波

一、中值滤波

二、快速中值滤波

介绍

原理

优化

三、代码

MATLAB

C

Python

四、测试

其他

（MATLAB/C/Python）快速中值滤波相关推荐

最新文章

热门文章