新年没看见什么新气象，可能雾霾太大，真个望长城内外，惟余莽莽，于是想起写一篇有关去雾的博客吧。图像去雾的最大挑战是：去雾的同时，能够保持图像不偏色，并且速度要快。本文算法来自《 Optimized contrast enhancement for real-time image and video dehazing》这篇文章，在实际测试中效果还可以。算法原理同样基于大气散射模型：
算法流程如下图所示：
1. 估算大气环境光：算法采用图像分块策略，递归查找大气环境光所在的图像块，定位到图像块后，在遍历图像块内像素，根据下述公式：
计算大气环境光数值。
void EstimateAirlight(BMPINFO *pSrcBitmap, uchar *airlight)
{
int blockMaxValIndex = 0;
float blockCurVal = 0.0f;
float blockMaxVal = 0.0f;
float mean[3] = { 0.0f };
float stddev[3] = { 0.0f };
float channelVal[3] = { 0.0f };
int width = pSrcBitmap->lWidth;
int height = pSrcBitmap->lHeight;
int size = width * height;
int halfWidth = width / 2;
int halfHeight = height / 2;
uchar *pSrcData = NULL;
uchar *pBlockData = NULL;
BMPINFO upperLeftBitmap = { 0 }, upperRightBitmap = { 0 }, lowerLeftBitmap = { 0 }, lowerRightBitmap = { 0 };
if (size > 400)
{
// 分块并填充数据
// 左上
upperLeftBitmap.dwPixelFormat = BMPFORMAT_RGB32_R8G8B8A8;
upperLeftBitmap.lWidth = halfWidth;
upperLeftBitmap.lHeight = halfHeight;
upperLeftBitmap.lPitch[0] = upperLeftBitmap.lWidth * 4;
upperLeftBitmap.pPlane[0] = (uchar *)malloc(upperLeftBitmap.lPitch[0] * upperLeftBitmap.lHeight);
pSrcData = pSrcBitmap->pPlane[0];
pBlockData = upperLeftBitmap.pPlane[0];
for (int i = 0; i < upperLeftBitmap.lHeight; i++)
{
memcpy(pBlockData, pSrcData, upperLeftBitmap.lPitch[0]);
pBlockData += upperLeftBitmap.lPitch[0];
pSrcData += pSrcBitmap->lPitch[0];
}
// 右上
upperRightBitmap.dwPixelFormat = BMPFORMAT_RGB32_R8G8B8A8;
upperRightBitmap.lWidth = width - halfWidth;
upperRightBitmap.lHeight = halfHeight;
upperRightBitmap.lPitch[0] = upperRightBitmap.lWidth * 4;
upperRightBitmap.pPlane[0] = (uchar *)malloc(upperRightBitmap.lPitch[0] * upperRightBitmap.lHeight);
pSrcData = pSrcBitmap->pPlane[0] + halfWidth*4;
pBlockData = upperRightBitmap.pPlane[0];
for (int i = 0; i < upperRightBitmap.lHeight; i++)
{
memcpy(pBlockData, pSrcData, upperRightBitmap.lPitch[0]);
pBlockData += upperRightBitmap.lPitch[0];
pSrcData += pSrcBitmap->lPitch[0];
}
// 左下
lowerLeftBitmap.dwPixelFormat = BMPFORMAT_RGB32_R8G8B8A8;
lowerLeftBitmap.lWidth = halfWidth;
lowerLeftBitmap.lHeight = height - halfHeight;
lowerLeftBitmap.lPitch[0] = lowerLeftBitmap.lWidth * 4;
lowerLeftBitmap.pPlane[0] = (uchar *)malloc(lowerLeftBitmap.lPitch[0] * lowerLeftBitmap.lHeight);
pSrcData = pSrcBitmap->pPlane[0] + pSrcBitmap->lPitch[0]*halfHeight;
pBlockData = lowerLeftBitmap.pPlane[0];
for (int i = 0; i < lowerLeftBitmap.lHeight; i++)
{
memcpy(pBlockData, pSrcData, lowerLeftBitmap.lPitch[0]);
pBlockData += lowerLeftBitmap.lPitch[0];
pSrcData += pSrcBitmap->lPitch[0];
}
// 右下
lowerRightBitmap.dwPixelFormat = BMPFORMAT_RGB32_R8G8B8A8;
lowerRightBitmap.lWidth = width - halfWidth;
lowerRightBitmap.lHeight = height - halfHeight;
lowerRightBitmap.lPitch[0] = lowerRightBitmap.lWidth * 4;
lowerRightBitmap.pPlane[0] = (uchar *)malloc(lowerRightBitmap.lPitch[0] * lowerRightBitmap.lHeight);
pSrcData = pSrcBitmap->pPlane[0] + pSrcBitmap->lPitch[0]*halfHeight + halfWidth*4;
pBlockData = lowerRightBitmap.pPlane[0];
for (int i = 0; i < lowerRightBitmap.lHeight; i++)
{
memcpy(pBlockData, pSrcData, lowerRightBitmap.lPitch[0]);
pBlockData += lowerRightBitmap.lPitch[0];
pSrcData += pSrcBitmap->lPitch[0];
}
// 计算最优块值
// 左上
CompMean_StdDev(&upperLeftBitmap, mean, stddev);
channelVal[0] = mean[0] - stddev[0];
channelVal[1] = mean[1] - stddev[1];
channelVal[2] = mean[2] - stddev[2];
blockCurVal = channelVal[0] + channelVal[1] + channelVal[2];
blockMaxVal = blockCurVal;
blockMaxValIndex = 0;
// 右上
CompMean_StdDev(&upperRightBitmap, mean, stddev);
channelVal[0] = mean[0] - stddev[0];
channelVal[1] = mean[1] - stddev[1];
channelVal[2] = mean[2] - stddev[2];
blockCurVal = channelVal[0] + channelVal[1] + channelVal[2];
if (blockCurVal > blockMaxVal)
{
blockMaxVal = blockCurVal;
blockMaxValIndex = 1;
}
// 左下
CompMean_StdDev(&lowerLeftBitmap, mean, stddev);
channelVal[0] = mean[0] - stddev[0];
channelVal[1] = mean[1] - stddev[1];
channelVal[2] = mean[2] - stddev[2];
blockCurVal = channelVal[0] + channelVal[1] + channelVal[2];
if (blockCurVal > blockMaxVal)
{
blockMaxVal = blockCurVal;
blockMaxValIndex = 2;
}
// 右下
CompMean_StdDev(&lowerRightBitmap, mean, stddev);
channelVal[0] = mean[0] - stddev[0];
channelVal[1] = mean[1] - stddev[1];
channelVal[2] = mean[2] - stddev[2];
blockCurVal = channelVal[0] + channelVal[1] + channelVal[2];
if (blockCurVal > blockMaxVal)
{
blockMaxVal = blockCurVal;
blockMaxValIndex = 3;
}
// 递归查找大气环境光所在的块
switch (blockMaxValIndex)
{
case 0:
free(upperRightBitmap.pPlane[0]);
upperRightBitmap.pPlane[0] = NULL;
free(lowerLeftBitmap.pPlane[0]);
lowerLeftBitmap.pPlane[0] = NULL;
free(lowerRightBitmap.pPlane[0]);
lowerRightBitmap.pPlane[0] = NULL;
EstimateAirlight(&upperLeftBitmap, airlight);
break;
case 1:
free(upperLeftBitmap.pPlane[0]);
upperLeftBitmap.pPlane[0] = NULL;
free(lowerLeftBitmap.pPlane[0]);
lowerLeftBitmap.pPlane[0] = NULL;
free(lowerRightBitmap.pPlane[0]);
lowerRightBitmap.pPlane[0] = NULL;
EstimateAirlight(&upperRightBitmap, airlight);
break;
case 2:
free(upperLeftBitmap.pPlane[0]);
upperLeftBitmap.pPlane[0] = NULL;
free(upperRightBitmap.pPlane[0]);
upperRightBitmap.pPlane[0] = NULL;
free(lowerRightBitmap.pPlane[0]);
lowerRightBitmap.pPlane[0] = NULL;
EstimateAirlight(&lowerLeftBitmap, airlight);
break;
case 3:
free(upperLeftBitmap.pPlane[0]);
upperLeftBitmap.pPlane[0] = NULL;
free(upperRightBitmap.pPlane[0]);
upperRightBitmap.pPlane[0] = NULL;
free(lowerLeftBitmap.pPlane[0]);
lowerLeftBitmap.pPlane[0] = NULL;
EstimateAirlight(&lowerRightBitmap, airlight);
break;
default:
break;
}
}
else
{
int minDistance = 99999;
int eucDistance = 0;
pSrcData = pSrcBitmap->pPlane[0];
// 估算大气环境光
for (int i = 0; i < size; i++, pSrcData += 4)
{
eucDistance = int(sqrt((255.0f - pSrcData[AXJ_BLUE])*(255.0f - pSrcData[AXJ_BLUE]) + (255.0f - pSrcData[AXJ_GREEN])*(255.0f - pSrcData[AXJ_GREEN]) + (255.0f - pSrcData[AXJ_RED])*(255.0f - pSrcData[AXJ_RED])));
if (eucDistance < minDistance)
{
minDistance = eucDistance;
airlight[AXJ_BLUE] = pSrcData[AXJ_BLUE];
airlight[AXJ_GREEN] = pSrcData[AXJ_GREEN];
airlight[AXJ_RED] = pSrcData[AXJ_RED];
}
}
}
// 释放资源
if (upperLeftBitmap.pPlane[0] != NULL)
{
free(upperLeftBitmap.pPlane[0]);
upperLeftBitmap.pPlane[0] = NULL;
}
if (upperRightBitmap.pPlane[0] != NULL)
{
free(upperRightBitmap.pPlane[0]);
upperRightBitmap.pPlane[0] = NULL;
}
if (lowerLeftBitmap.pPlane[0] != NULL)
{
free(lowerLeftBitmap.pPlane[0]);
lowerLeftBitmap.pPlane[0] = NULL;
}
if (lowerRightBitmap.pPlane[0] != NULL)
{
free(lowerRightBitmap.pPlane[0]);
lowerRightBitmap.pPlane[0] = NULL;
}
}
2. 估算大气透射率：这应该是本文的核心部分了，所谓优化对比度去雾，即尽量提高有雾图像的对比度，同时尽量减少图像信息损失。优化对比度计算公式：
信息量损失计算公式：
最后选取合适的参数，以获取最优透射率图：
3. 修正透射率图：本文在计算透射率图时，是以图像块而非像素点为单位计算，因此计算得到的透射率图类似马赛克效果，于是作者采用导向滤波的方式修正透射率图，即以原图的灰度图为导向图，以透射率图为滤波输入图，导向滤波结果即为修正的透射率图。导向滤波也是个挺神奇的算法，以后再讲。另外需要注意的是图像块大小也将影响最后的去雾效果。
4. 复原图像：大气环境光和透射率图，根据前面提到的大气散射模型就可以复原图像了。复原过程中可以适当做一下gamma校正，避免复原后的图像过于昏暗。下面是去雾主代码示例，自己在实现算法时，在效果和速度上肯定需要做些权衡，所以算法流程和原算法相比还是有些差异的。
void ImageDeHazing(BMPINFO *pSrcBitmap)
{
// 原图下采样
int width = 1000, height = 750;
if (pSrcBitmap->lWidth > pSrcBitmap->lHeight)
{
width = (pSrcBitmap->lWidth > 1000) ? 1000 : 500;
height = width*pSrcBitmap->lHeight / pSrcBitmap->lWidth;
}
else
{
height = (pSrcBitmap->lHeight > 1000) ? 1000 : 500;
width = height*pSrcBitmap->lWidth / pSrcBitmap->lHeight;
}
int size = width * height;
BMPINFO dSamplingBmp = { 0 };
dSamplingBmp.dwPixelFormat = BMPFORMAT_RGB32_R8G8B8A8;
dSamplingBmp.lWidth = pSrcBitmap->lWidth;
dSamplingBmp.lHeight = pSrcBitmap->lHeight;
dSamplingBmp.lPitch[0] = dSamplingBmp.lWidth * 4;
dSamplingBmp.pPlane[0] = (uchar*)malloc(dSamplingBmp.lPitch[0] * dSamplingBmp.lHeight);
memcpy(dSamplingBmp.pPlane[0], pSrcBitmap->pPlane[0], dSamplingBmp.lPitch[0]*pSrcBitmap->lHeight);
ImageResampling(&dSamplingBmp, width, height);
float *pGrayData = (float *)malloc(size * sizeof(float));
uchar *pSrcData = NULL, *pDstData = NULL;
pSrcData = dSamplingBmp.pPlane[0];
for (int i = 0; i < size; i++, pSrcData += 4)
{
pGrayData[i] = (pSrcData[AXJ_BLUE] + pSrcData[AXJ_GREEN] + pSrcData[AXJ_RED]) / 3.0f;
}
// 估算大气环境光
uchar airlight[3] = { 0 };
EstimateAirlight(&dSamplingBmp, airlight);
// 估算大气透射率
const int blockSize = 5;
float *pAirTransmission = (float *)malloc(width * height * sizeof(float));
EstimateTransmission(&dSamplingBmp, airlight, pAirTransmission, blockSize, blockSize);
// 修正透射率图
BMPINFO transbmp = { 0 };
transbmp.dwPixelFormat = BMPFORMAT_GRAY8;
transbmp.lWidth = width;
transbmp.lHeight = height;
transbmp.lPitch[0] = transbmp.lWidth;
transbmp.pPlane[0] = (uchar *)malloc(transbmp.lPitch[0] * transbmp.lHeight);
FastGuidedFilter(pGrayData, pAirTransmission, transbmp.pPlane[0], width, height, 40, 0.01f);
// 释放资源
free(pGrayData);
pGrayData = NULL;
free(pAirTransmission);
pAirTransmission = NULL;
free(dSamplingBmp.pPlane[0]);
dSamplingBmp.pPlane[0] = NULL;
// 输出校正
uchar ret_gammaLut[256];
for(int i = 0; i < 256; i++)
{
ret_gammaLut[i] = (uchar)(pow((float)i / 255.0f, 0.65f) * 255.0f);
}
// 复原图像
float *horLookupTable = (float *)malloc(pSrcBitmap->lWidth * sizeof(float));
float *verLookupTable = (float *)malloc(pSrcBitmap->lHeight * sizeof(float));
for (int i = 0; i < pSrcBitmap->lWidth; i++)
{
horLookupTable[i] = ((float)i / pSrcBitmap->lWidth) * width;
}
for (int i = 0; i < pSrcBitmap->lHeight; i++)
{
verLookupTable[i] = ((float)i / pSrcBitmap->lHeight) * height;
}
uchar transmission = 0;
pSrcData = pSrcBitmap->pPlane[0];
for (int i = 0; i < pSrcBitmap->lHeight; i++)
{
for (int j = 0; j < pSrcBitmap->lWidth; j++, pSrcData += 4)
{
BilinearInterGray(transbmp.pPlane[0], horLookupTable[j], verLookupTable[i], transbmp.lWidth, transbmp.lHeight, &transmission);
pSrcData[AXJ_BLUE] = ret_gammaLut[(uchar)CLAMP0255((float)(pSrcData[AXJ_BLUE] - airlight[AXJ_BLUE])*255 / transmission + airlight[AXJ_BLUE])];
pSrcData[AXJ_GREEN] = ret_gammaLut[(uchar)CLAMP0255((float)(pSrcData[AXJ_GREEN] - airlight[AXJ_GREEN])*255 / transmission + airlight[AXJ_GREEN])];
pSrcData[AXJ_RED] = ret_gammaLut[(uchar)CLAMP0255((float)(pSrcData[AXJ_RED] - airlight[AXJ_RED])*255 / transmission + airlight[AXJ_RED])];
}
}
// 释放资源
free(horLookupTable);
horLookupTable = NULL;
free(verLookupTable);
verLookupTable = NULL;
free(transbmp.pPlane[0]);
transbmp.pPlane[0] = NULL;
}下面是一些去雾效果图：
原文献代码下载链接：http://download.csdn.net/detail/u013085897/9776286
参考资料：
9792
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
