图像去雾哪家强？之前我们已经讨论过了著名的基于暗通道先验的图像去雾（Kaiming He, 2009）算法，如果你用兴趣可以参考：
此外，网上也有很多同道推荐了一篇由韩国学者所发表的研究论文《Optimized contrast enhancement for real-time image and video dehazing》（你也可以从文末参考文献【1】给出的链接中下载到这篇经典论文），其中原作者就提出了一个效果相当不错的图像去雾算法。最近有朋友在我们的图像处理算法研究学习群中也提到了该算法，恰巧想到主页君也很久未发图像相关之文章了，今天遂心血来潮，向大家科普一下这个新算法。
欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。
算法核心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-block
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;
// upper right sub-block
cvCvtPixToPlane(iplUpperRight, iplR, iplG, iplB, 0);
cvMean_StdDev(iplR, dpMean, dpStds);
cvMean_StdDev(iplG, dpMean+1, dpStds+1);
cvMean_StdDev(iplB, dpMean+2, dpStds+2);
dpScore[0] = dpMean[0]-dpStds[0];
dpScore[1] = dpMean[1]-dpStds[1];
dpScore[2] = dpMean[2]-dpStds[2];
afScore[1] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);
if(afScore[1] > nMaxScore)
{
nMaxScore = afScore[1];
nMaxIndex = 1;
}
// lower left sub-block
cvCvtPixToPlane(iplLowerLeft, iplR, iplG, iplB, 0);
cvMean_StdDev(iplR, dpMean, dpStds);
cvMean_StdDev(iplG, dpMean+1, dpStds+1);
cvMean_StdDev(iplB, dpMean+2, dpStds+2);
dpScore[0] = dpMean[0]-dpStds[0];
dpScore[1] = dpMean[1]-dpStds[1];
dpScore[2] = dpMean[2]-dpStds[2];
afScore[2] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);
if(afScore[2] > nMaxScore)
{
nMaxScore = afScore[2];
nMaxIndex = 2;
}
// lower right sub-block
cvCvtPixToPlane(iplLowerRight, iplR, iplG, iplB, 0);
cvMean_StdDev(iplR, dpMean, dpStds);
cvMean_StdDev(iplG, dpMean+1, dpStds+1);
cvMean_StdDev(iplB, dpMean+2, dpStds+2);
dpScore[0] = dpMean[0]-dpStds[0];
dpScore[1] = dpMean[1]-dpStds[1];
dpScore[2] = dpMean[2]-dpStds[2];
afScore[3] = (float)(dpScore[0]+dpScore[1]+dpScore[2]);
if(afScore[3] > nMaxScore)
{
nMaxScore = afScore[3];
nMaxIndex = 3;
}
// select the sub-block, which has maximum score
switch (nMaxIndex)
{
case 0:
AirlightEstimation(iplUpperLeft); break;
case 1:
AirlightEstimation(iplUpperRight); break;
case 2:
AirlightEstimation(iplLowerLeft); break;
case 3:
AirlightEstimation(iplLowerRight); break;
}
}
else
{
// select the atmospheric light value in the sub-block
for(nY=0; nY<nHei; nY++)
{
for(nX=0; nX<nWid; nX++)
{
// 255-r, 255-g, 255-b
nDistance = int(sqrt(float(255-(uchar)imInput->imageData[nY*nStep+nX*3])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3])
+float(255-(uchar)imInput->imageData[nY*nStep+nX*3+1])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3+1])
+float(255-(uchar)imInput->imageData[nY*nStep+nX*3+2])*float(255-(uchar)imInput->imageData[nY*nStep+nX*3+2])));
if(nMinDistance > nDistance)
{
nMinDistance = nDistance;
m_anAirlight[0] = (uchar)imInput->imageData[nY*nStep+nX*3];
m_anAirlight[1] = (uchar)imInput->imageData[nY*nStep+nX*3+1];
m_anAirlight[2] = (uchar)imInput->imageData[nY*nStep+nX*3+2];
}
}
}
}
cvReleaseImage(&iplUpperLeft);
cvReleaseImage(&iplUpperRight);
cvReleaseImage(&iplLowerLeft);
cvReleaseImage(&iplLowerRight);
cvReleaseImage(&iplR);
cvReleaseImage(&iplG);
cvReleaseImage(&iplB);
}
或者你也可以使用下面这个Matlab的实现，显而易见代码更加简单（源码由图像算法研究学习群里“薛定谔的猫”提供，略有修改）：
function airlight = est_airlight(img)
%compute atmospheric light A through hierarchical
%searching method based on the quad-tree subdivision
global best;
[w,h,z] = size(img);
img = double(img);
if w*h > 200
lu = img(1:floor(w/2),1:floor(h/2),:);
ru = img(1:floor(w/2),floor(h/2):h,:);
lb = img(floor(w/2):w,1:floor(h/2),:);
rb = img(floor(w/2):w,floor(h/2):h,:);
lu_m_r = mean(mean(lu(:,:,1)));
lu_m_g = mean(mean(lu(:,:,2)));
lu_m_b = mean(mean(lu(:,:,3)));
ru_m_r = mean(mean(ru(:,:,1)));
ru_m_g = mean(mean(ru(:,:,2)));
ru_m_b = mean(mean(ru(:,:,3)));
lb_m_r = mean(mean(lb(:,:,1)));
lb_m_g = mean(mean(lb(:,:,2)));
lb_m_b = mean(mean(lb(:,:,3)));
rb_m_r = mean(mean(rb(:,:,1)));
rb_m_g = mean(mean(rb(:,:,2)));
rb_m_b = mean(mean(rb(:,:,3)));
lu_s_r = std2(lu(:,:,1));
lu_s_g = std2(lu(:,:,2));
lu_s_b = std2(lu(:,:,3));
ru_s_r = std2(ru(:,:,1));
ru_s_g = std2(ru(:,:,2));
ru_s_b = std2(ru(:,:,3));
lb_s_r = std2(lb(:,:,1));
lb_s_g = std2(lb(:,:,2));
lb_s_b = std2(lb(:,:,3));
rb_s_r = std2(rb(:,:,1));
rb_s_g = std2(rb(:,:,2));
rb_s_b = std2(rb(:,:,3));
score0 = lu_m_r + lu_m_g + lu_m_b - lu_s_r - lu_s_g - lu_s_b;
score1 = ru_m_r + ru_m_g + ru_m_b - ru_s_r - ru_s_g - ru_s_b;
score2 = lb_m_r + lb_m_g + lb_m_b - lb_s_r - lb_s_g - lb_s_b;
score3 = rb_m_r + rb_m_g + rb_m_b - rb_s_r - rb_s_g - rb_s_b;
x = [score0,score1,score2,score3];
if max(x) == score0
est_airlight(lu);
elseif max(x) == score1
est_airlight(ru);
elseif max(x) == score2
est_airlight(lb);
elseif max(x) == score3
est_airlight(rb);
end
else
for i = 1:w
for j = 1:h
nMinDistance = 65536;
distance = sqrt((255 - img(i,j,1)).^2 + (255 - img(i,j,2)).^2 + (255 - img(i,j,3)).^2);
if nMinDistance > distance
nMinDistance = distance;
best = img(i,j,:);
end
end
end
end
airlight =best;
end
算法核心2：透射率的计算
我们首先给出图像对比度度量的方法（论文中，原作者给出了三个对比度定义式，我们只讨论其中第一个）：
CMSE=∑p=1N=(Jc(p)?Jˉc)2N
其中 c∈{r,g,b} 是颜色通道的索引标签，Jˉc 是Jc(p)的平均值，并且p=1,?,N，N 是块中像素的数量。
根据之前给出的有雾图像与原始（没有雾的）图像之间的关系模型
J(p)=1t(I(p)?A)+A
我们可以把上述对比度定义式重新为
CMSE=∑p=1N=(Ic(p)?Iˉc)2t2N
其中 Iˉc 是Ic(p)的平均值，而且你会发现上述式子也告诉我们对比度是关于 t 的递减函数。
既然我们希望通过增强对比度的方法来去雾，那么不妨将一个区块B内三个颜色通道上的MSE对比度加总，然后再取负，如下
Econtrast=?∑c∈{r,g,b}∑p∈B(Jc(p)?Jˉc)2NB=?∑c∈{r,g,b}∑p∈B(Ic(p)?Iˉc)2t2NB
由于加了负号，所以取对比度最大就等同于取上式最小。
另外一方面，因为对比度得到增强，可能会导致部分像素的调整值超出了0和255的范围，这样就会造成信息的损失以及视觉上的瑕疵。所以算法作者又提出了一个信息量损失的计算公式：
于是我们把所有问题都统一到了求下面这个式子的最小值问题上
E=Econtrast+λLEloss
其中，λL 是一个权重参数用于控制信息损失和对比度之间的一个相对重要性。With a small value λL, the restored images have significantly increased contrast, but they lose information and contain unnaturally dark pixels due to the truncation of pixel values. On the contrary, with a large value of λL, we can prevent the information loss but cannot remove haze fully. 总的来说, λL= 5 strikes a balance between the information loss prevention and the haze removal effectively.
下面是基于OpenCV实现的示例代码：
float dehazing::NFTrsEstimationColor(int *pnImageR, int *pnImageG, int *pnImageB, int nStartX, int nStartY, int nWid, int nHei)
{
int nCounter;
int nX, nY;
int nEndX;
int nEndY;
int nOutR, nOutG, nOutB;
int nSquaredOut;
int nSumofOuts;
int nSumofSquaredOuts;
float fTrans, fOptTrs;
int nTrans;
int nSumofSLoss;
float fCost, fMinCost, fMean;
int nNumberofPixels, nLossCount;
nEndX = __min(nStartX+m_nTBlockSize, nWid);
nEndY = __min(nStartY+m_nTBlockSize, nHei);
nNumberofPixels = (nEndY-nStartY)*(nEndX-nStartX) * 3;
fTrans = 0.3f;
nTrans = 427;
for(nCounter=0; nCounter<7; nCounter++)
{
nSumofSLoss = 0;
nLossCount = 0;
nSumofSquaredOuts = 0;
nSumofOuts = 0;
for(nY=nStartY; nY<nEndY; nY++)
{
for(nX=nStartX; nX<nEndX; nX++)
{
nOutB = ((pnImageB[nY*nWid+nX] - m_anAirlight[0])*nTrans + 128*m_anAirlight[0])>>7; // (I-A)/t + A --> ((I-A)*k*128 + A*128)/128
nOutG = ((pnImageG[nY*nWid+nX] - m_anAirlight[1])*nTrans + 128*m_anAirlight[1])>>7;
nOutR = ((pnImageR[nY*nWid+nX] - m_anAirlight[2])*nTrans + 128*m_anAirlight[2])>>7;
if(nOutR>255)
{
nSumofSLoss += (nOutR - 255)*(nOutR - 255);
nLossCount++;
}
else if(nOutR < 0)
{
nSumofSLoss += nOutR * nOutR;
nLossCount++;
}
if(nOutG>255)
{
nSumofSLoss += (nOutG - 255)*(nOutG - 255);
nLossCount++;
}
else if(nOutG < 0)
{
nSumofSLoss += nOutG * nOutG;
nLossCount++;
}
if(nOutB>255)
{
nSumofSLoss += (nOutB - 255)*(nOutB - 255);
nLossCount++;
}
else if(nOutB < 0)
{
nSumofSLoss += nOutB * nOutB;
nLossCount++;
}
nSumofSquaredOuts += nOutB * nOutB + nOutR * nOutR + nOutG * nOutG;;
nSumofOuts += nOutR + nOutG + nOutB;
}
}
fMean = (float)(nSumofOuts)/(float)(nNumberofPixels);
fCost = m_fLambda1 * (float)nSumofSLoss/(float)(nNumberofPixels)
- ((float)nSumofSquaredOuts/(float)nNumberofPixels - fMean*fMean);
if(nCounter==0 || fMinCost > fCost)
{
fMinCost = fCost;
fOptTrs = fTrans;
}
fTrans += 0.1f;
nTrans = (int)(1.0f/fTrans*128.0f);
}
return fOptTrs;
}
在原文中，作者还提供了并行实现的算法（如果你参考作者的源码，你会发现他使用了OpenMP）用于实时对视频进行去雾，这已经超出本文的研究范围，我们不做深究。
实验效果
作者的项目主页（文献[2]）中提供有一个基于OpenCV的实现，不过由于OpenCV代码配置起来太麻烦，所以我并未采用。非常感谢图像算法研究群里的同学分享给我MATLAB版的代码（我略有修改）。其实这个算法用MATLAB来写真的非常方便，大约不超过150行即可搞定，特别是MATLAB2014之后集成了导向滤波算法，所以一个函数调用直接搞定，省去很多麻烦。下面是我基于这个MATLAB版程序给出一些测试结果，可供参考（由于我和原作者所使用的参数可能存在差异，所以效果并不保证完全相同）。
主要结论：
- 算法对于天空部分的处理相当到位，更优于Kaiming He的暗通道先验算法；
- 对比度过大时，图像很容易发暗，可以后期将图片稍微调亮一些。
- 算法本身是从对比增强的角度来进行去雾操作的，所以你可以看出结果自动带对比度增强加成，所以这个算法所取得的结果通常更加鲜亮。
===========
无冥冥之志者，无昭昭之明，无惛惛之事者，无赫赫之功。
===========
如果你是图像处理的同道中人，欢迎加入图像处理学习群（529549320）。为保证本群质量，入群前请先阅读群规（即本博客置顶文章http://blog.csdn.Net/baimafujinji/article/details/50570976），并在博客置顶帖中留言，否则将不予入群，请不要做无谓的尝试。
参考文献与推荐阅读材料
[1] Kim, J.H., Jang, W.D., Sim, J.Y. and Kim, C.S., 2013. Optimized contrast enhancement for real-time image and video dehazing. Journal of Visual Communication and Image Representation, 24(3), pp.410-425. （百度云下载链接）
[2] 原作者的项目主页——可以下载到基于OpenCV实现的程序源代码
[3] laviewpbt更早之前发布的一篇研究该算法的文章——该博客上同时包含有很多其他图像去雾方面的文章，推荐有兴趣的读者参考