子块部分重叠直方图均衡算法（POSHE算法）原理解读及c++实现(一)

本文链接：https://blog.csdn.net/u012294613/article/details/140999668

写在前面
Matlab版代码：https://blog.csdn.net/weixin_40647819/article/details/89603660

最近刚好写到了直方图均衡算法，因为之前用到过图像增强，就大致地再多了解了一下，看到了POSHE算法，这个算法也算是比较经典的吧，有一些它的优势。其实这篇论文是很老的论文了，是韩国人提出的，目前看来可能意义不大，但是我看了一下引用情况，中文论文有很多基于这个算法各种改进，增强，去雾等等，英文文献的引用也有大概800吧，说明它的原理还是比较重要或者说有借鉴之处的。

论文全名：An advanced contrast enhancement using partially overlapped sub-block histogram equalization System

论文地址：http://medialab.sjtu.edu.cn/teaching/DIP/Projects/chapter_enh/BlockHE.pdf

因为刚实现的时候，文章的表述有些地方我理解的不太清楚，就在网上搜了一下，发现不知道为啥很少有人写这个算法，代码除了一个Matlab版的（我大致看了一下，写的不完整，不知道有没有用），几乎没有其他的版本，原理也只是有人简单的介绍一下，问问题的和要程序的倒是有一些。所以决定写一下吧，然后用c++实现一下，有些地方也可能理解的不到位。

下面对子块部分重叠算法(POSHE算法)作一下简介吧。总的来说，直方图均衡算法可以分为两类，一类是全局直方图均衡，也就是我们熟悉的经典直方图均衡算法，全局的缺点在于某些局部会损失细节；另一类是局部直方图均衡算法，它能够较好得保持局部细节的对比度。

局部直方图均衡算法，又称为子块直方图均衡算法，可按照所均衡子块的重叠程度来分类，可以分为子块不重叠、子块重叠与子块部分重叠三种。

子块不重叠、子块重叠算法都有较大的局限性，子块不重叠算法会导致不可避免的块效应，子块重叠算法虽然克服了块效应，但子块是逐像素移动的，效率很低，从应用来看这两种基本都不会被采用的。而POSHE算法刚好能克服这两点，所以我就看了这篇论文，并实现了（只能说视觉效果还过的去）。

子块部分重叠的均衡算法（POSHE）

该方法（POSHE）与其他子块重叠方法的不同之处在于：

（1）子块不是逐像素移动，而是将移动步长约取为子块尺寸的几分之一。

（2）子块均衡的灰度转换函数不仅用于映射子块中心像素灰度值，而且用于映射子块所有像素的灰度值。

（3）对多次被均衡的像素，将均衡结果取平均作为该像素在输出图像中的灰度值

子块部分重叠算法的优势在于：

（1）由于子块部分重叠方式减少了相邻子块间的均衡函数形状差异，使块效应基本得以消除对于子块边界可能出现的少量块效应，论文中设计了块效应消除滤波器（BERF），可以视觉上克服这种边界上的少量块效应。

（2）由于子块均衡总次数比子块重叠方式少得多，计算效率大幅度提高。

（3）图像细节的增强能力与子块重叠算法相近。

原理
它的原理分为两个部分，先执行POSHE，然后执行BERF。执行POSHE倒是很简单的，后面执行BERF，论文写的含糊不清（是我理解能力有限），搞了我几个小时。

执行POSHE
1.基本思想

为了既能够尽可能削减子块不重叠算法的块效应，还要能够克服子块重叠算法的效率问题，作者采取了一个折中的办法，基本思想就是：将子块移动的步长取子块尺寸的一部分（几分之一），子块在图像上移动，然后对每个子块进行直方图均衡，因为子块是部分重叠，所以有很多像素会被均衡很多次，称为均衡频次，且不同区域的像素均衡频次不同，作者就对这些被均衡的像素进行累加，最后再除以其均衡频次。比如某个像素总共被处理了64次，最后再除以64，并以这个结果作为最终结果。

2.子块和掩码

为了削弱块效应，文中首先打了一个比方，提到了3*3的掩码（mask），如下图，这很好理解，就类似于高斯滤波模板之类的，可以起到平滑作用。使用这个3*3的掩码（mask），中心区域的变换函数就可由其自身及其八个相邻区域（权重）的直方图得到。随着这个掩码尺寸的增加，块效应就会被持续削弱，文章认为一般当尺寸在15*15以上时，块效应就会得到充分抑制，达到

与子块重叠算法相媲美的效果。不过我后面的结果发现15*15以上的依然有块效应，而且与子块的尺寸有关，子块尺寸越小，块效应越明显；也与图像背景也有关，细节纹理比较丰富的图像（树林）块效应不太明显，而相反那种简单背景（天空）的更容易产生块效应。

这种掩码可以通过子块的移动产生，比如3*3的mask，就可以通过子块移动它尺寸的1/2得到。对于一个120*120的子块，移动步长是60，则得到3*3的mask。进一步，如果移动子块尺寸的1/4，则可得到7*7的mask；移动子块尺寸的1/8，就可得到15*15的mask.......也就是说可以通过改变子块的移动步长来控制mask的大小。下面是一个例子：

虚线框是一个子块，显然这是产生的3*3的mask，因为子块移动步长是子块尺寸的1/2，当子块移动完成，白色区域的像素的均衡频次是1，灰色区域是2，深色区域是4，想象一下画一下不难理解。这些被多次处理的像素是累加的，最后除以均衡频次得到映射值。

放几张图片吧，是子块移动形成的网格（不同尺寸的子块和不同的步长，具体多少忘了）：

3.理论推导

当然，虽然过程简单，但也是有理论基础的，我就搬运一下论文的分析过程吧，首先给了一个图：

还是3*3的mask，当然根据前面说的，a,b,c,d......i可不是像素，而是子块的子区域。子块移动四次，所以有1,2,3,4个子块。我们关注中间e区域的POSHE函数，因为e区域被均衡了4次（4个子块），所以有：

而每一个子块的直方图均衡变换函数为：

然后经过一系列推导，就得到了：

发现还真神奇呢，它的系数和3*3高斯滤波的模板一样吧，哈哈。

4.代码实现

这个算法实现，有一点不好就是很多地方需要考虑图像的边界问题，执行POSHE还好，执行BERF是时时刻刻要考虑啊。

1：图像边界扩充，因为子块移动嘛，要保证子块不越能界嘛。不扩充也行，但是不均匀，为了计算方便我还是扩充了，用的是对称法，就是边界附近的像素对称到要扩充的区域；

2：设置子块尺寸，输出图像初始化为0；

3：子块移动，对当前子块进行直方图均衡，记录均衡频次，结果在输出图像中累加，直至移动完成整个图像；（这里累加要注意数据类型，8bit的最大值是255，不排除累加输出会超过255的，所以输出图像可以先设置16位的，或者浮点型的）

4：输出图像图像每个像素除以对应的均衡频次(频次可以用一个矩阵表示)；

5：数据类型转换，转到8bit图像；

6：裁剪，因为之前扩充了边界嘛，完成。

感受：最繁琐的就是子块移动吧，要细心，不能越界，再就要找准参考点，我是以子块左顶点为参考点，从0开始，然后下一个子块的左顶点就是当前左顶点+移动步长，最后判断这个左顶点是不是最边界框的左顶点，当然也可以以右顶点为参考。用for循就行了。

结果
先看结果图吧，在附代码吧。

原图仅仅执行POSHE 执行了BERF

基本就是这个结果吧，仅仅执行POSHE ，会发现有块效应，边缘出现了先，鸭子图片还不怎么明显，胸部CT图片就很明显了，执行了BERF之后，块效应的线基本消除了。饿了，没时间了........下一篇再写BERF吧，BERF代码有点多，但都是相似的。

再看看POSHE和经典全局的对比吧：

全局处理 POSHE

好像区别不大吼，仔细看鸭子头部、左上角房子、右边的树干，POSHE算法的效果在细节增强恢复上还是好一些的。执行效率方面下一篇在对比一下吧。

代码
里面的BERF(newsrc, dst, step_r, step_c, step_levels, large_dir)，是执行BERF。下次写博客再贴出来吧。

对了，全局直方图均衡的c++实现:https://blog.csdn.net/weixin_40647819/article/details/88383606

/***************************************
POSHE
src - 输入图像
dst - 输出图像
s  -  子块尺寸=原图尺寸/s ，2的倍数
k  -  子块移动步长=子块尺寸/k, 2的倍数
****************************************/
void Poshe(cv::Mat& src, cv::Mat& dst, cv::Mat& newsrc_draw,float s, float k){
    int nr = src.rows;
    int nc = src.cols;
 
    //边界扩充
    int newnr = ceil(nr / s / k)*s*k;
    int newnc = ceil(nc / s / k)*s*k;
    cv::Mat newsrc;
    cv::copyMakeBorder(src, newsrc, 0 , newnr - nr, 0, newnc - nc, cv::BORDER_REFLECT);
    dst=cv::Mat::zeros(newnr, newnc, CV_16U); //因为8位无符号整型的范围是0~255，而在累加过程中会超过这个范围，所以采用16位无符号整型
 
 
    //创建子块，确定子块移动步长
    int sub_block_r = newnr / s; //子块高度
    int sub_block_c = newnc / s; //子块宽度
    int step_r = sub_block_r / k; //垂直移动步长
    int step_c = sub_block_c / k; //水平移动步长
 
    //对当前子块进行直方图均衡
    int sub_block_x ; //子块左顶点坐标(行)
    int sub_block_y; //子块左顶点坐标(列)
    cv::Mat HE_frequency = cv::Mat::zeros(newnr, newnc, CV_8U); //均衡频率计数矩阵
    cv::Mat sub_block_HE;
    newsrc.copyTo(newsrc_draw);
    for (sub_block_x=0; sub_block_x <= (newnr - sub_block_r); ){
        for (sub_block_y = 0; sub_block_y <= (newnc - sub_block_c); ){
            cv::Mat sub_block = newsrc(cv::Rect(sub_block_y, sub_block_x, sub_block_c, sub_block_r));
            Histogram_equalization(sub_block, sub_block_HE);
            cv::rectangle(newsrc_draw, cv::Rect(sub_block_y, sub_block_x, sub_block_c, sub_block_r), cv::Scalar(0, 0, 0), 1.8, 1, 0);
            
            //将直方图均衡后子块的像素值映射至输出图像    
            int sub_block_HE_i = 0;
            for (int i = sub_block_x; i < sub_block_x + sub_block_r; i++){    
                 int sub_block_HE_j = 0;
                for (int j = sub_block_y; j < sub_block_y + sub_block_c; j++){     
                     dst.at<ushort>(i, j) = dst.at<ushort>(i, j) + sub_block_HE.at<uchar>(sub_block_HE_i, sub_block_HE_j);
                     HE_frequency.at<uchar>(i, j)++;
                     sub_block_HE_j++;
                }
                sub_block_HE_i++;    
            }
 
            sub_block_y = sub_block_y + step_c;
        } 
        sub_block_x = sub_block_x + step_r;
    }
 
    for (int i = 0; i < dst.rows; ++i){
        for (int j = 0; j < dst.cols; ++j){
            dst.at<ushort>(i, j) = (dst.at<ushort>(i, j)) / (HE_frequency.at<uchar>(i, j));
        }
    }
    dst.convertTo(dst, CV_8U, 1, 0); //数据类型转换
 
    //BERF
    int step_levels = 3; 
    int large_dir = 40;
    BERF(newsrc, dst, step_r, step_c, step_levels, large_dir);
 
    dst = dst(cv::Rect(0, 0, src.cols, src.rows));    
}

————————————————

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

原文链接：https://blog.csdn.net/weixin_40647819/article/details/88416512

接着上次的博客写，上次已经介绍完了POSHE部分，也为这次写BERF做了一些铺垫。在这次正式开始之前，还是先放两张图片，以说明执行BERF的必要性：

仅执行POSHE

执行BERF后
从上面两幅图中就可以看出，仅执行POSHE会产生不可避免的快效应，尤其是边界产生了线。为了抑制这种干扰，论文提出来一种名为BERF(blocking effect reduction filter)的滤波器，中文叫块效应抑制滤波器。

原理
步骤：
本质上就是一种空间域的滤波，原理比较简单，我就直接写步骤：

1、判断是否产生块效应。论文没有给出具体的判断方法，按照我的理解，就是计算经过POSHE处理后图像的子块边界相邻两侧两个像素的灰度差异f1，同时计算原图的子块边界相邻两侧两个像素的灰度差异f2，如果f1-f2大于s，则产生了块效应，s可取3~4。比如，边界某像素在图像矩阵第5行，则边界相邻两侧两个像素就分别是第4和第6行。

2、若产生块效应，则计算第一步中相邻两侧两个像素的灰度平均值，并将该边界像素的灰度值替换为平均值。

3、在垂直于边界的方向上，以递增/递减顺序的像素值被缓慢变化的值替换，而且而是从上步计算的平均值开始。简单地说，就是在垂直于边界的方向上，从平均值开始，将本来正在处于递增或递减的像素值，代替为缓慢递增或递减的像素值。比如，边界为第5行，平均灰度值为100，那么这个边界像素的灰度值被替换为100，将另外两侧相邻像素中较大的像素替换为104，较小的替换为96，然后沿着垂直于边界的方向，如果不满足终止条件，就会一直替换下去，递增的一侧就会是104、108、112、116....,递减的一侧就会是96、92、88、84.........，直至满足终止条件。

终止条件：
论文用文字对3个终止条件进行了描述，但是我觉得用文字描述很不直观：

我用公式描述吧，以行边界（垂直边界方向就是列方向咯），且像素向上递增为例：

用x1表示当前已经执行过BERF的像素的灰度值（论文写的pixels that have previously have previously...）,x2表示下一个准备执行BERF的像素的灰度值（论文写的the present pixel）。

条件1：x2-x1<4（step size），step size可取4及以下的值；

条件2：x2-x1> 40（large enough），large enough可取40及以上的值；

条件3：x2>X1。

以上就是以行边界（垂直边界方向就是列方向咯），且像素向上递增情况下的终止条件，我们编程时当然用的是执行条件咯：

以上3个终止条件改写成执行条件就是：

4<=x2-x1<=40;

很简单了有木有。其他情况的条件可以依次类推，或者看我的程序。

代码
结果图开篇已经放了，直接把代码放上吧，代码有点多，但都是相似的，搞懂一种情况即可：

PS:这里只放了BERF的代码，整个可以运行的代码在：github

/***************************************
BERF滤波器
src - 输入图像
dst - 输出图像
step_r  - 子块垂直移动步长
step_c  - 子块水平移动步长
step_levels - 灰度级缓慢变化的步长
large_dir  - 最大灰度差异
****************************************/
void BERF(cv::Mat& src, cv::Mat& dst, int step_r, int step_c , int step_levels, int large_dir){
 
    //处理子块的行边界（横向边界）
    for (int i = 0; i < dst.rows;){
        for (int j = 0; j < dst.cols; ++j){
            if (i - 1 >= 0 && i + 1 <= dst.rows - 1){ //图像边界判断，先排除第一行和最后一行
                //块效应检查
                int dfacter_newsrc = abs(src.at<uchar>(i, j) - src.at<uchar>(i + 1, j)) + abs(src.at<uchar>(i, j) - src.at<uchar>(i - 1, j)); //原图子块边界与相邻两侧像素的灰度差异
                int dfacter_dst = abs(dst.at<uchar>(i, j) - dst.at<uchar>(i + 1, j)) + abs(dst.at<uchar>(i, j) - dst.at<uchar>(i - 1, j)); //POSHEed图子块边界与相邻两侧像素的灰度差异
                if (dfacter_dst - dfacter_newsrc > step_levels){
                    //存在块效应执行BERF
                    int b = 0;
                    int ave_bound = (int)(dst.at<uchar>(i - 1, j) + dst.at<uchar>(i + 1, j)) / 2;
                    dst.at<uchar>(i, j) = ave_bound;
 
                    if (dst.at<uchar>(i - 1, j) > dst.at<uchar>(i + 1, j)){ //垂直于子块边界，向上递增方向（ increasing rule）
                        if (i - 2 - b >= 0){ //图像边界判断(下一位置像素)
                            int pixel_present_add = dst.at<uchar>(i - 1 - b, j);
                            int pixel_next_add = dst.at<uchar>(i - 2 - b, j);
                            pixel_present_add = ave_bound + step_levels;
                            while (i - 2 - b >= 0 && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir){
                                pixel_next_add = pixel_present_add + step_levels;
                                b += 1;
                                if (i - 2 - b >= 0){ //图像边界判断(下一位置像素)
                                    pixel_present_add = dst.at<uchar>(i - 1 - b, j);
                                    pixel_next_add = dst.at<uchar>(i - 2 - b, j);
                                }
                            }
                        }
                        //垂直于子块边界，向下递减方向（ decreasing rule）
                        if (i + 2 + b <dst.rows){ //图像边界判断
                            int pixel_present_dec = dst.at<uchar>(i + 1 + b, j);
                            int pixel_next_dec = dst.at<uchar>(i + 2 + b, j);
                            pixel_present_dec = ave_bound - step_levels;
                            while (i + 2 + b >= 0 && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir){
                                pixel_next_dec = pixel_present_dec - step_levels;
                                b += 1;
                                if (i + 2 + b >= 0){ //图像边界判断
                                    pixel_present_dec = dst.at<uchar>(i + 1 + b, j);
                                    pixel_next_dec = dst.at<uchar>(i + 2 + b, j);
                                }
                            }
                        }
 
                    }
                    else if (dst.at<uchar>(i - 1, j) < dst.at<uchar>(i + 1, j)){ 垂直于子块边界，向下递增方向（ increasing rule）
                        if (i + 2 + b < dst.rows){ //图像边界判断
                            int pixel_present_add = dst.at<uchar>(i + 1 + b, j);
                            int pixel_next_add = dst.at<uchar>(i + 2 + b, j);
                            pixel_present_add = ave_bound + step_levels;
                            while (i + 2 + b >= 0 && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir){
                                pixel_next_add = pixel_present_add + step_levels;
                                b += 1;
                                if (i + 2 + b >= 0){ //图像边界判断
                                    pixel_present_add = dst.at<uchar>(i + 1 + b, j);
                                    pixel_next_add = dst.at<uchar>(i + 2 + b, j);
                                }
                            }
                        }
                        //垂直于子块边界，向上递减方向（ decreasing rule）
                        if (i - 2 - b >= 0){ //图像边界判断
                            int pixel_present_dec = dst.at<uchar>(i - 1 - b, j);
                            int pixel_next_dec = dst.at<uchar>(i - 2 - b, j);
                            pixel_present_dec = ave_bound - step_levels;
                            while (i - 2 - b >= 0 && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir){
                                pixel_next_dec = pixel_present_dec - step_levels;
                                b += 1;
                                if (i - 2 - b >= 0){ //图像边界判断
                                    pixel_present_dec = dst.at<uchar>(i - 1 - b, j);
                                    pixel_next_dec = dst.at<uchar>(i - 2 - b, j);
                                }
                            }
                        }
 
                    }
                }
            }
        }
        i += step_r;
    }
 
    //处理子块的列边界（纵向向边界）
    for (int j = 0; j <dst.cols;){
        for (int i = 0; i < dst.rows; ++i){
            if (j - 1 >= 0 && j + 1 <= dst.cols - 1){ //图像边界判断，先排除第一列和最后一列
                //块效应检查
                int dfacter_newsrc = abs(src.at<uchar>(i, j) - src.at<uchar>(i, j + 1)) + abs(src.at<uchar>(i, j) - src.at<uchar>(i, j - 1));
                int dfacter_dst = abs(dst.at<uchar>(i, j) - dst.at<uchar>(i, j + 1)) + abs(dst.at<uchar>(i, j) - dst.at<uchar>(i, j + 1));
                if (dfacter_dst - dfacter_newsrc > step_levels){
                    //存在块效应执行BERF
                    int b = 0;
                    int ave_bound = (int)(dst.at<uchar>(i, j - 1) + dst.at<uchar>(i, j + 1)) / 2;
                    dst.at<uchar>(i, j) = ave_bound;
 
                    if (dst.at<uchar>(i, j - 1) > dst.at<uchar>(i, j + 1)){ //垂直于子块边界，向左递增方向（ increasing rule）
                        if (j - 2 - b >= 0){ //图像边界判断
                            int pixel_present_add = dst.at<uchar>(i, j - 1 - b);
                            int pixel_next_add = dst.at<uchar>(i, j - 2 - b);
                            pixel_present_add = ave_bound + step_levels;
                            while (j - 2 - b >= 0 && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir){
                                pixel_next_add = pixel_present_add + step_levels;
                                b += 1;
                                if (j - 2 - b >= 0){ //图像边界判断
                                    pixel_present_add = dst.at<uchar>(i, j - 1 - b);
                                    pixel_next_add = dst.at<uchar>(i, j - 2 - b);
                                }
                            }
                        }
                        //垂直于子块边界，向右递减方向（ decreasing rule）
                        if (j + 2 + b <dst.cols){ //图像边界判断
                            int pixel_present_dec = dst.at<uchar>(i, j + 1 + b);
                            int pixel_next_dec = dst.at<uchar>(i, j + 2 + b);
                            pixel_present_dec = ave_bound - step_levels;
                            while (j + 2 + b >= 0 && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir){
                                pixel_next_dec = pixel_present_dec - step_levels;
                                b += 1;
                                if (j + 2 + b >= 0){ //图像边界判断
                                    pixel_present_dec = dst.at<uchar>(i, j + 1 + b);
                                    pixel_next_dec = dst.at<uchar>(i, j + 2 + b);
                                }
                            }
                        }
 
                    }
                    else if (dst.at<uchar>(i, j - 1) < dst.at<uchar>(i, j + 1)){ 垂直于子块边界，向右递增方向（ increasing rule）
                        if (j + 2 + b < dst.cols){ //图像边界判断
                            int pixel_present_add = dst.at<uchar>(i, j + 1 + b);
                            int pixel_next_add = dst.at<uchar>(i, j + 1 + b);
                            pixel_present_add = ave_bound + step_levels;
                            while (j + 2 + b >= 0 && pixel_next_add - pixel_present_add >= step_levels && pixel_next_add - pixel_present_add <= large_dir){
                                pixel_next_add = pixel_present_add + step_levels;
                                b += 1;
                                if (j + 2 + b >= 0){ //图像边界判断
                                    pixel_present_add = dst.at<uchar>(i, j + 1 + b);
                                    pixel_next_add = dst.at<uchar>(i, j + 2 + b);
                                }
                            }
                        }
                        //垂直于子块边界，向左递减方向（ decreasing rule）
                        if (j - 2 - b >= 0){ //图像边界判断
                            int pixel_present_dec = dst.at<uchar>(i, j - 1 - b);
                            int pixel_next_dec = dst.at<uchar>(i, j - 2 - b);
                            pixel_present_dec = ave_bound - step_levels;
                            while (j - 2 - b >= 0 && pixel_present_dec - pixel_next_dec >= step_levels && pixel_present_dec - pixel_next_dec <= large_dir){
                                pixel_next_dec = pixel_present_dec - step_levels;
                                b += 1;
                                if (j - 2 - b >= 0){ //图像边界判断
                                    pixel_present_dec = dst.at<uchar>(i, j - 1 - b);
                                    pixel_next_dec = dst.at<uchar>(i, j - 2 - b);
                                }
                            }
                        }
                    }
                }
            }
        }
        j += step_c;
    }
}

————————————————

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

原文链接：https://blog.csdn.net/weixin_40647819/article/details/88768435