首页 > 代码库 > 用于图像去雾的优化对比度增强算法
用于图像去雾的优化对比度增强算法
图像去雾哪家强?之前我们已经讨论过了著名的基于暗通道先验的图像去雾(Kaiming He, 2009)算法,如果你用兴趣可以参考:
- 暗通道优先的图像去雾算法(上)
- 暗通道优先的图像去雾算法(下)
此外,网上也有很多同道推荐了一篇由韩国学者所发表的研究论文《Optimized contrast enhancement for real-time image and video dehazing》(你也可以从文末参考文献【1】给出的链接中下载到这篇经典论文),其中原作者就提出了一个效果相当不错的图像去雾算法。最近有朋友在我们的图像处理算法研究学习群中也提到了该算法,恰巧想到主页君也很久未发图像相关之文章了,今天遂心血来潮,向大家科普一下这个新算法。
欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji,为保证公式、图表得以正确显示,强烈建议你从该地址上查看原版博文。本博客主要关注方向包括:数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。
算法核心1:计算大气光值
通常,图像去雾问题的基本模型可以用下面这个公式来表示(这一点在基于暗通道先验的图像去雾中我们也使用过):
其中,
此外,
大气光
基于这个认识,算法的作者提出了一个基于四叉树子空间划分的层次搜索方法。如下图所示,我们首先把输入图像划分成四个矩形区域。然后,为每个子区域进行评分,这个评分的计算方法是“用区域内像素的平均值减去这些像素的标准差”(the average pixel value subtracted by the standard deviation
of the pixel values within the region)。记下来,选择具有最高得分的区域,并将其继续划分为更小的四个子矩形。我们重复这个过程直到被选中的区域小于某个提前指定的阈值。例如下图中的红色部分就是最终被选定的区域。在这被选定的区域里,我们选择使得距离
我们假设在一个局部的小范围内,场景深度是相同的(也就是场景内的各点到相机镜头的距离相同),所以在一个小块内(例如
可见,在求得大气光
总的来说,一个有雾的块内,对比度都是比较低的,而被恢复的块内的对比度则随着
下面是原作者给出的大气光计算函数代码(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:透射率的计算
我们首先给出图像对比度度量的方法(论文中,原作者给出了三个对比度定义式,我们只讨论其中第一个):
其中
根据之前给出的有雾图像与原始(没有雾的)图像之间的关系模型
我们可以把上述对比度定义式重新为
其中
既然我们希望通过增强对比度的方法来去雾,那么不妨将一个区块B内三个颜色通道上的MSE对比度加总,然后再取负,如下
由于加了负号,所以取对比度最大就等同于取上式最小。
另外一方面,因为对比度得到增强,可能会导致部分像素的调整值超出了0和255的范围,这样就会造成信息的损失以及视觉上的瑕疵。所以算法作者又提出了一个信息量损失的计算公式:
于是我们把所有问题都统一到了求下面这个式子的最小值问题上
其中,
下面是基于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更早之前发布的一篇研究该算法的文章——该博客上同时包含有很多其他图像去雾方面的文章,推荐有兴趣的读者参考
用于图像去雾的优化对比度增强算法