用于图像去雾的优化对比度增强算法.docx
《用于图像去雾的优化对比度增强算法.docx》由会员分享，可在线阅读，更多相关《用于图像去雾的优化对比度增强算法.docx（17页珍藏版）》请在冰豆网上搜索。
用于图像去雾的优化对比度增强算法
用于图像去雾的优化对比度增强算法
图像去雾哪家强？
之前我们已经讨论过了著名的基于暗通道先验的图像去雾（KaimingHe,2009）算法，如果你用兴趣可以参考：
此外，网上也有很多同道推荐了一篇由韩国学者所发表的研究论文《Optimizedcontrastenhancementforreal-timeimageandvideodehazing》（你也可以从文末参考文献【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通常用图像中最明亮的颜色来作为估计。
因为大量的灰霾通常会导致一个发亮（发白）的颜色。
然而，在这个框架下，那些颜色比大气光更加明亮的物体通常会被选中，因而便会导致一个本来不应该作为大气光参考值的结果被用作大气光的估计。
为了更加可靠的对大气光进行估计，算法的作者利用了这样一个事实：
通常，那些灰蒙蒙的区域（也就是天空）中像素的方差（或者变动）总体来说就比较小。
基于这个认识，算法的作者提出了一个基于四叉树子空间划分的层次搜索方法。
如下图所示，我们首先把输入图像划分成四个矩形区域。
然后，为每个子区域进行评分，这个评分的计算方法是“用区域内像素的平均值减去这些像素的标准差”（theaveragepixelvaluesubtractedbythestandarddeviation
ofthepixelvalueswithintheregion）。
记下来，选择具有最高得分的区域，并将其继续划分为更小的四个子矩形。
我们重复这个过程直到被选中的区域小于某个提前指定的阈值。
例如下图中的红色部分就是最终被选定的区域。
在这被选定的区域里，我们选择使得距离||（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++版）
voiddehazing:
:
AirlightEstimation（IplImage*imInput）
{
intnMinDistance=65536;
intnDistance;
intnX,nY;
intnMaxIndex;
doubledpScore[3];
doubledpMean[3];
doubledpStds[3];
floatafMean[4]={0};
floatafScore[4]={0};
floatnMaxScore=0;
intnWid=imInput->width;
intnHei=imInput->height;
intnStep=imInput->widthStep;
//4sub-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）;
//comparetothreshold（200）-->biggerthanthreshold,dividetheblock
if（nHei*nWid>200）
{
//computethemeanandstd-devinthesub-block
//upperleftsub-ock
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;
//upperrightsub-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;
}
//lowerleftsub-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;
}
//lowerrightsub-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;
}
//selectthesub-block,whichhasmaximumscore
switch（nMaxIndex）
{
case0:
AirlightEstimation（iplUpperLeft）;break;
case1:
AirlightEstimation（iplUpperRight）;break;
case2:
AirlightEstimation（iplLowerLeft）;break;
case3:
AirlightEstimation（iplLowerRight）;break;
}
}
else
{
//selecttheatmosphericlightvalueinthesub-block
for（nY=0;nY{
for（nX=0;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的实现，显而易见代码更加简单（源码由图像算法研究学习群里“薛定谔的猫”提供，略有修改）：
functionairlight=est_airlight（img）
%computeatmosphericlightAthroughhierarchical
%searchingmethodbasedonthequad-treesubdivision
globalbest;
[w,h,z]=size（img）;
img=double（img）;
ifw*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];
ifmax（x）==score0
est_airlight（lu）;
elseifmax（x）==score1
est_airlight（ru）;
elseifmax（x）==score2
est_airlight（lb）;
elseifmax（x）==score3
est_airlight（rb）;
end
else
fori=1:
w
forj=1:
h
nMinDistance=65536;
distance=sqrt（（255-img（i,j,1））.^2+（255-img（i,j,2））.^2+（255-img（i,j,3））.^2）;
ifnMinDistance>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是一个权重参数用于控制信息损失和对比度之间的一个相对重要性。
WithasmallvalueλL,therestoredimageshavesignificantlyincreasedcontrast,buttheyloseinformationandcontainunnaturallydarkpixelsduetothetruncationofpixelvalues.Onthecontrary,withalargevalueofλL,wecanpreventtheinformationlossbutcannotremovehazefully.总的来说,λL=5strikesabalancebetweentheinformationlosspreventionandthehazeremovaleffectively.
下面是基于OpenCV实现的示例代码：
floatdehazing:
:
NFTrsEstimationColor（int*pnImageR,int*pnImageG,int*pnImageB,intnStartX,intnStartY,intnWid,intnHei）
{
intnCounter;
intnX,nY;
intnEndX;
intnEndY;
intnOutR,nOutG,nOutB;
intnSquaredOut;
intnSumofOuts;
intnSumofSquaredOuts;
floatfTrans,fOptTrs;
intnTrans;
intnSumofSLoss;
floatfCost,fMinCost,fMean;
intnNumberofPixels,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++）
升级会员 