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

最新推荐文章于 2023-01-17 08:39:12 发布

今天不飞了

最新推荐文章于 2023-01-17 08:39:12 发布

阅读量3.4w

点赞数 10

分类专栏：算法实现 MATLAB C/C++ 文章标签：中值滤波快速 MATLAB C/C++ Python

本文链接：https://blog.csdn.net/xsz591541060/article/details/108280557

版权

MATLAB 同时被 3 个专栏收录

80 篇文章 100 订阅

订阅专栏

算法实现

11 篇文章 6 订阅

订阅专栏

C/C++

10 篇文章 1 订阅

订阅专栏

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

一、中值滤波
二、快速中值滤波
三、代码
四、测试
其他

by HPC_ZY

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

一、中值滤波

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

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

代码见后文。

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

二、快速中值滤波

介绍

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

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

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

原理

实现方法
step 1：统计 $N * N$ 窗口内像素
step 2：计算累积直方图
step 3：找到累积直方图值刚好超过 $\frac{N*N}{2}$ 时对应的灰度值。
通过一个简单的例子更直观的说明。
这里假设一个只有8个灰度级的图像，某个 $5 * 5$ 的窗口内值如下，
$\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}$
其统计直方图如下表。对于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$ 、窗口尺寸为 $N * N$ ，我们只需要进行以下操作：
1）进行 $N * N$ 次统计（获得直方图）
2）最多进行 $L - 1$ 次累加+判断（计算累积直方图的同时判断是否超过半数）
通常图像的 $L$ 都为256，我们可以推断对于 $3 * 3$ 的窗口，传统排序方法速度更快；但当N越来越大，直方图法的优势就越来越大。

代码见后文

优化

进一步思考可以发现，我们在平移窗口时，新窗口里有 $N * (N - 1)$ 的像素都是上一窗口里有的（如下图黄色），如果重新统计就浪费时间。
所以有人提出不必重新计算直方图，只用更新直方图——移除（如下图红色）+加入（如下图绿色）。图像来自论《Median Filter in Constant Time》}
代码见后文

三、代码

本文首先用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循环的时候运行速度特别慢，如果有大佬知道原因希望留言告诉我，谢谢了。

今天不飞了

关注

10
点赞
踩
56

收藏

觉得还不错? 一键收藏
0
评论
（MATLAB/C/Python）快速中值滤波

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

扫一扫