# 用于图像去雾的优化对比度增强算法

## 算法核心1：计算大气光值

I(p)=t(p)J(p)+(1t(p))A

of the pixel values within the region）。记下来，选择具有最高得分的区域，并将其继续划分为更小的四个子矩形。我们重复这个过程直到被选中的区域小于某个提前指定的阈值。例如下图中的红色部分就是最终被选定的区域。在这被选定的区域里，我们选择使得距离 ||(Ir(p),Ig(p),Ib(p))(255,255,255)||$||(I_r(p),I_g(p),I_b(p))-(255,255,255)||$ 最小化的颜色（包含 r,g,b$r,g,b$ 三个分量）来作为大气光的参考值。注意，这样做的意义在于我们希望选择那个离纯白色最近的颜色（也就是最亮的颜色）来作为大气光的参考值。

J(p)=1t(I(p)A)+A

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);
}




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

J(p)=1t(I(p)A)+A

CMSE=p=1N=(Ic(p)I¯c)2t2N

Econtrast=c{r,g,b}pB(Jc(p)J¯c)2NB=c{r,g,b}pB(Ic(p)I¯c)2t2NB

E=Econtrast+λLEloss

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;
}


## 实验效果

• 算法对于天空部分的处理相当到位，更优于Kaiming He的暗通道先验算法；
• 对比度过大时，图像很容易发暗，可以后期将图片稍微调亮一些。
• 算法本身是从对比增强的角度来进行去雾操作的，所以你可以看出结果自动带对比度增强加成，所以这个算法所取得的结果通常更加鲜亮。

===========

## 参考文献与推荐阅读材料

[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更早之前发布的一篇研究该算法的文章——该博客上同时包含有很多其他图像去雾方面的文章，推荐有兴趣的读者参考

• 广告
• 抄袭
• 版权
• 政治
• 色情
• 无意义
• 其他

120