快速中值滤波算法

中值滤波算法 ：

在图像处理中，在进行如边缘检测这样的进一步处理之前，通常需要首先进行一定程度的降噪。 中值滤波 是一种非线性数字滤波器技术，经常用于去除图像或者其它信号中的噪声。这个设计思想就是检查输入信号中的采样并判断它是否代表了信号，使用奇数个采样组成的观察窗实现这项功能。观察窗口中的数值进行排序，位于观察窗中间的中值作为输出。然后，丢弃最早的值，取得新的采样，重复上面的计算过程。中值滤波是图像处理中的一个常用步骤，它对于 斑点噪声和椒盐噪声来说尤其有用。保存边缘的特性使它在不希望出现边缘模糊的场合也很有用。

为了演示中值滤波器的工作过程，我们给下面的数组加上观察窗 3 ，重复边界的数值：

x = [2 80 6 3]

y[1] = Median[2 2 80] = 2
y[2] = Median[2 80 6] = Median[2 6 80] = 6
y[3] = Median[80 6 3] = Median[3 6 80] = 6
y[4] = Median[6 3 3] = Median[3 3 6] = 3

于是
y = [2 6 6 3]

其中 y 是 x 的中值滤波输出。

普通中值滤波算法伪代码：

Input: image X of size m*n, kernel radius r.

output: image Y as X.

for i = r to m - r do

for j = r to n - r do

initialize list A[]

for a = i-r to i+r

for b = j-r to j+r

add X(a, b) to  A[]

sort A[] then Y(i ,j) = A[A.size/2]

但是，上述算法在像素处理处的复杂度为O(r ² ).

OpenCV实现代码：

#include "cv.h"
#include "highgui.h"
#include <iostream>
using namespace std;
using namespace cv;
int main(int argc, char* argv[])
        Mat src = imread("beauty.jpg");
        Mat dst;
        //参数是按顺序写的
        //高斯滤波
        //src:输入图像
        //dst:输出图像
        //Size(5,5)模板大小，为奇数
        //x方向方差
        //Y方向方差
        GaussianBlur(src,dst,Size(5,5),0,0);
        imwrite("gauss.jpg",dst);
        //中值滤波
        //src:输入图像
        //dst::输出图像
        //模板宽度，为奇数
        medianBlur(src,dst,3);
        imwrite("med.jpg",dst);
        //均值滤波
        //src:输入图像
        //dst:输出图像
        //模板大小
        //Point(-1,-1):被平滑点位置，为负值取核中心
        blur(src,dst,Size(3,3),Point(-1,-1));
        imwrite("mean.jpg",dst);
        //双边滤波
        //src:输入图像
        //dst:输入图像
        //滤波模板半径
        //颜色空间标准差
        //坐标空间标准差
        bilateralFilter(src,dst,5,10.0,2.0);//这里滤波没什么效果，不明白
        imwrite("bil.jpg",dst);
        waitKey();
        return 0;
View Code
　　快速中值滤波算法：
　　　　　O(r)复杂度的Huang算法：<A Fast Two-Dimensional Median Filtering Algorithm>
　　　　这个代码的核心在于维护一个kernel直方图，可以实现快速的读取和删除扫描区域的像素值。
int StopAmount = (2 * Radius + 1) * (2*Radius +1) * Percentile;
int Sum = 0;
for (int I = 0; I <= 255; I++)
    Sum += HistBlue(I);
    if (Sum >= StopAmount)      // 满足停止条件
        Value = I;              // 得到中值
        break; 
 　　　在MFC中实现的huang算法：
/******************************************************************
* 功能: 彩色图像的快速中值滤波平滑处理
* 参数: image0为原图形，image1平滑结果，
*        w、h为图象的宽和高
*       size为进行平滑的邻域边长
******************************************************************/
void FastSmoothMedianCl(BYTE* image0, BYTE* image1, unsigned int w, unsigned int h, unsigned int size)
    //将图像转化为矩阵形式
    BYTE** imageBuf0 = CreatImage(image0, w, h);
    BYTE** imageBuf1 = CreatImage(image1, w, h);
    //设定模板
    int x,y,c;
    int a;
    //int scale;
    int ** templt;
    //初始化 x = size/2 为中心的点的直方图，
    //根据邻域大小设定模板
    templt = new int*[5];
    for(c=0; c<3; c++)
        templt[c] = new int [256];
    for(y=size/2; y<h-size/2; y++)
        for(c=0; c<3; c++)
            for (int i= 0; i< 256; i ++){
                templt[c][i] = 0;
            for(int i = y-size/2; i < y+size/2; i++){
                for (int j = 0; j < size; j++){
                    int k = imageBuf0[i][j*4+c];
                    templt[c][k] ++;
        for(x=size/2+1; x<w-size/2; x++)
            for(c=0; c<3; c++)
                //取采样窗口中像素灰度的中值
                a=FastMedianValueCl(imageBuf0,w,h,templt,size,x,y,c);
                //a/= scale;
                //过限处理
                a = a>255?255:a;    
                a = a<0?0:a;
                imageBuf1[y][x*4+c]=a;
    //清理内存
    for(int i = 0; i < 3; i ++){
        delete templt[i];
    free(imageBuf0);
    free(imageBuf1);
/**************************************************
* 功能: 使用直方图对彩色图邻域获取中值
* 参数: imageBuf为目标图像 w、h为图像大小
*       templt为模板 tw为邻域大小 
*        x,y为当前采样窗口中心像素的坐标
*       cn为颜色分量编号 0为蓝色 1为绿色 2为红色
**************************************************/
int FastMedianValueCl(BYTE** imageBuf0, int w, int h, int**templt, int tw, int x, int y, int cn)
    int i,j,k;                   
    int px,py, mid;
    int count = 0;                 //用来保存采样窗口的像素数量
    int count_temp = 0;
    k=0;
    //从采样窗口中取得像素灰度
    for(i=0; i<tw; i++)
            py=y-tw/2+i;
            px=x+tw/2+1;
            //如果该像素位于采样窗口
            //保存像素灰度
            k = imageBuf0[py][px*4+cn];
            templt[cn][k]++;
            px=x-tw/2;
            //保存像素灰度
            k = imageBuf0[py][px*4+cn];
            templt[cn][k]--;
    for(int a = 0; a < 256; a++){
        mid = a;
        count_temp += templt[cn][a];
        if(count_temp > (tw*tw)/2) break;
    return mid;
View Code
　　 　与普通算法的比较：
　　　　O(1)复杂度的CTMF算法：<Median Filter in Constant Time.pdf>
　　首先，对于每一列图像，我们都为其维护一个直方图（对于8位图像，该直方图有256个元素），在整个的处理过程中，这些直方图数据都必须得到维护。每列直方图累积了2r+1个垂直方向上相邻像素的信息，初始的时候，这2r+1个像素是分别以第一行的每个像素为中心的。核的直方图通过累积2r+1个相邻的列直方图数据获取。其实，我们所做的就是将核直方图分解成他对应的列直方图的集合，在整个滤波的过程中，这些直方图数据在两个步骤内用恒定的时间保持最新。
　　考虑从某个像素向右移动一个像素的情况。对于当前行，核最右侧的列直方图首先需要更新，而此时该列的列直方图中的数据还是以上一行对应位置那个像素为中心计算的。因此需要减去最上一个像素对应的直方图然后加上其下面一像素的直方图信息。这样做的效果就是将列直方图数据降低一行。这一步很明显是个0(1)操作，只有一次加法和一次减法，而于半径r无关。
      第二步更新核直方图，其是2r+1个列直方图之和。这是通过减去最左侧的列直方图数据，然后再加上第一步所处理的那一列的列直方图数据获得的。这一步也是个O(1)操作，如图2所示。如前所述，加法、减法以及计算直方图的中值的耗时都是一些依赖于图像位深的计算，而于滤波半径无关。
　　　　这个算法较之前的Huang算法的差别在于，充分利用了之前的数据，这有点类似于动态规划，也是用空间换取时间的策略。
 参考博客：http://www.cnblogs.com/Imageshop/archive/2013/04/26/3045672.html