基于记忆性的中值滤波O(r)与O(1)复杂度的算法实现

　　本文参考博客：https://www.cnblogs.com/Imageshop/archive/2013/04/26/3045672.html

　　原生的中值滤波是基于排序算法的，这样的算法复杂度基本在O(r²)左右，当滤波半径较大时，排序算法就显得很慢。对此有多种改进算法，这里介绍经典

的Huang算法与O(1)算法，两者都是基于记忆性的算法，只是后者记性更强。

　　排序算法明显的一个不足之处就是无记忆性。当核向右移动一列后，只是核的最左和最右列数据发生了变化，中间不变的数据应当被存储起来，而排序算法

并没有做到这点。Huang算法的思想是建立一个核直方图，来统计核内的各灰度的像素数。当核向右移动时，就将新的一列所有数据加入到直方图中，同时将最

左列的旧数据从直方图中删除，如下图所示。这样做使得大部分数据能够被记忆，减少重复操作。当直方图更新完毕后，就可以通过从左到右累计像素来找到中值。

　　　基于记忆性的中值滤波O(r)与O(1)复杂度的算法实现　　　　　　　　　　　

下面是算法具体实现步骤与代码：

　　1.每行开始都将直方图、像素计数、中值变量清零，将核覆盖的所有像素加到直方图中。

　　2.计算中值，sumcnt为小于中值灰度的像素数和。如果当前sumcnt大于等于阈值，则表明，实际中值比当前median小，则直方图向左减去像素数，同时median

　　　也减小，直到sumcnt小于阈值，则此时median为中值（sumcnt加上中值像素数将大于或等于阈值）；如果sumcnt加上median像素数仍小于阈值，表明，实

　　　际中值比当前median大，则直方图向右累加像素数，同时median增加，直到sumcnt加上median像素数大于等于阈值，则此时median为中值（sumcnt小于阈值）。

　　3.核每向右移动一列，更新核直方图，比较每一个新加入的灰度值与median，小于则median加1，比较每一个减去的灰度值与median，小于则median减1。执行

　　　一次步骤2后继续移动核，直到整幅图都被滤波完毕。

%中值滤波黄法
%输入：8位深度图像、滤波半径（3*3的半径为1,中值百分比）
function grayimg=medianfilterhuang(grayimg,radius,percentage)
    [sizex,sizey]=size(grayimg);
    thres=int16(((2*radius+1)^2*percentage));
    %扩展边界
    extendimg=int16(horzcat(repmat(grayimg(:,1),1,radius),grayimg,repmat(grayimg(:,sizey),1,radius)));
    extendimg=int16(vertcat(repmat(extendimg(1,:),radius,1),extendimg,repmat(extendimg(sizex,:),radius,1)));
    for i=radius+1:sizex+radius
        %初始化直方图，每行第一个将所有像素加入直方图
        histogram=int16(zeros(256,1));
        sumcnt=int16(0);
        median=int16(0);
        for j=radius+1:sizey+radius
            %遍历卷积覆盖行
            for k=i-radius:i+radius
                if (j>radius+1)
                    %删去左列，加入新右列
                    histogram(extendimg(k,j-radius-1)+1)=histogram(extendimg(k,j-radius-1)+1)-1;
                    histogram(extendimg(k,j+radius)+1)=histogram(extendimg(k,j+radius)+1)+1;
                    %只关心中值灰度以左的像素数量
                    if (extendimg(k,j-radius-1)<median) 
                        sumcnt=sumcnt-1;
                    end
                    if (extendimg(k,j+radius)<median)
                        sumcnt=sumcnt+1;
                    end
                else
                    %行首
                    for h=j-radius:j+radius
                        histogram(extendimg(k,h)+1)=histogram(extendimg(k,h)+1)+1;
                    end
                end
            end
            %sumcnt不将中值像素数累加进去
            %旧中值偏大
            while(sumcnt>=thres)
                median=median-1;
                sumcnt=sumcnt-histogram(median+1);
            end
            %旧中值偏小
            while(sumcnt+histogram(median+1)<thres)
                sumcnt=sumcnt+histogram(median+1);
                median=median+1;
            end
            grayimg(i-radius,j-radius)=median;
        end
    end
end

Huang算法