写在前面

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                                      执行了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));	
}