用于图像去雾的优化对比度增强算法
图像去雾哪家强？之前我们已经讨论过了著名的基于暗通道先验的图像去雾（Kaiming He, 2009）算法，如果你用兴趣可以参考：
此外，网上也有很多同道推荐了一篇由韩国学者所发表的研究论文《Optimized contrast enhancement for real-time image and video dehazing》（你也可以从文末参考文献【1】给出的链接中下载到这篇经典论文），其中原作者就提出了一个效果相当不错的图像去雾算法。最近有朋友在我们的图像处理算法研究学习群中也提到了该算法，恰巧想到主页君也很久未发图像相关之文章了，今天遂心血来潮，向大家科普一下这个新算法。
为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。
算法核心1：计算大气光值
通常，图像去雾问题的基本模型可以用下面这个公式来表示（这一点在基于暗通道先验的图像去雾中我们也使用过）：
I(p)=t(p)J(p)+(1?t(p))A
其中，J(p)=(Jr(p),Jg(p),Jb(p))T 表示原始图像（也就是没有雾的图像）；I(p)=(Ir(p),Ig(p),Ib(p))T 表示我们观察到的图像（也就是有雾的图像）。r、g、b 表示位置 p 处的像素的三个分量。A=(Ar,Ag,Ab)T 是全球大气光，它表示周围环境中的大气光。
此外，t(p)∈[0,1] 是反射光的透射率, 由场景点到照相机镜头之间的距离所决定。因为光传播的距离越远，那么通常光就约分散而且越发被削弱。所以上面这个公式的意思就是，本来没有被雾所笼罩的图像 J 与大气光 A 按一定比例进行混合后就得到我们最终所观察到的有雾图像。
大气光 A 通常用图像中最明亮的颜色来作为估计。因为大量的灰霾通常会导致一个发亮（发白）的颜色。然而，在这个框架下，那些颜色比大气光更加明亮的物体通常会被选中，因而便会导致一个本来不应该作为大气光参考值的结果被用作大气光的估计。为了更加可靠的对大气光进行估计，算法的作者利用了这样一个事实：通常，那些灰蒙蒙的区域（也就是天空）中像素的方差（或者变动）总体来说就比较小。
基于这个认识，算法的作者提出了一个基于四叉树子空间划分的层次搜索方法。如下图所示，我们首先把输入图像划分成四个矩形区域。然后，为每个子区域进行评分，这个评分的计算方法是“用区域内像素的平均值减去这些像素的标准差”（the average pixel value subtracted by the standard deviation
of the pixel values within the region）。记下来，选择具有最高得分的区域，并将其继续划分为更小的四个子矩形。我们重复这个过程直到被选中的区域小于某个提前指定的阈值。例如下图中的红色部分就是最终被选定的区域。在这被选定的区域里，我们选择使得距离 ||(Ir(p),Ig(p),Ib(p))?(255,255,255)|| 最小化的颜色（包含 r,g,b 三个分量）来作为大气光的参考值。注意，这样做的意义在于我们希望选择那个离纯白色最近的颜色（也就是最亮的颜色）
来作为大气光的参考值。
我们假设在一个局部的小范围内，场景深度是相同的（也就是场景内的各点到相机镜头的距离相同），所以在一个小块内（例如32×32）我们就可以使用一个固定的透射率 t，所以前面给出的有雾图像与原始（没有雾的）图像之间的关系模型就可以改写为
J(p)=1t(I(p)?A)+A 可见，在求得大气光 A 的估计值之后，我们希望复原得到的原始（没有雾的）图像 J(p) 将依赖于散射率 t。
总的来说，一个有雾的块内，对比度都是比较低的，而被恢复的块内的对比度则随着 t 的估计值的变小而增大，我们将设法来估计一个最优的 t 值，从而使得去雾后的块能够得到最大的对比度。
下面是原作者给出的大气光计算函数代码（C/C++版）
void dehazing::AirlightEstimation(IplImage* imInput) {
int nMinDistance = 65536; int nDistance;
int nX, nY;
int nMaxIndex; double dpScore[3]; double dpMean[3]; double dpStds[3];
float afMean[4] = {0}; float afScore[4] = {0}; float nMaxScore = 0;
int nWid = imInput->width; int nHei = imInput->height;
int nStep = imInput->widthStep;
// 4 sub-block
IplImage *iplUpperLeft = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3); IplImage *iplUpperRight = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3);
IplImage *iplLowerLeft = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3); IplImage *iplLowerRight = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 3);
IplImage *iplR = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1); IplImage *iplG = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1); IplImage *iplB = cvCreateImage(cvSize(nWid/2, nHei/2),IPL_DEPTH_8U, 1);
// divide
cvSetImageROI(imInput, cvRect(0, 0, nWid/2, nHei/2)); cvCopyImage(imInput, iplUpperLeft);
cvSetImageROI(imInput, cvRect(nWid/2+nWid%2, 0, nWid, nHei/2)); cvCopyImage(imInput, iplUpperRight);
cvSetImageROI(imInput, cvRect(0, nHei/2+nHei%2, nWid/2, nHei)); cvCopyImage(imInput, iplLowerLeft);
cvSetImageROI(imInput, cvRect(nWid/2+nWid%2, nHei/2+nHei%2, nWid, nHei)); cvCopyImage(imInput, iplLowerRight);
// compare to threshold(200) --> bigger than threshold, divide the block if(nHei*nWid > 200) {
// compute the mean and std-dev in the sub-block // upper left sub-blwww.shanxiwang.netock
cvCvtPixToPlane(iplUpperLeft, iplR, iplG, iplB, 0);
cvMean_StdDev(iplR, dpMean, dpStds);
cvMean_StdDev(iplG, dpMean+1, dpStds+1); cvMean_StdDev(iplB, dpMean+2, dpStds+2); // dpScore: mean - std-dev
dpScore[0] = dpMean[0]-dpStds[0]; dpScore[1] = dpMean[1]-dpStds[1]; dpScore[2] = dpMean[2]-dpStds[2];
afScore[0] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);
nMaxScore = afScore[0]; nMaxIndex = 0;