基于暗通道先验的单幅图像去雾算法实现
最近闲来无事尝试了一下去雾算法的效果，从随便搜集到的资料来看多数算法是用于白天带雾图像去雾的，少部分是研究夜间带雾图像去雾的。而白天的去雾算法有很多都衍生自何博士的《Single Image Haze Removal Using Dark Channel Prior》，可见何博士算法之经典。我也只是把别人整理过的算法看懂，并写代码验证效果，并无创新，权当记录。
该算法用到了一个暗通道先验“Dark Channel Prior”，它的含义是：对于图片中绝大对数的非天空区域，总有在某些色彩通道具有很小的值，原作者测试了好多图片都证明该先验是正确的，暗通道的表达式如下：
式中Jc表示彩色图像的每个通道的值，Ω(x)表示以像素x为中心的一个窗口。可以看出获取暗通道就是获取一个窗口内所有像素的各个通道的最小值。
式（5）的意义用代码表达也很简单，首先求出每个像素RGB分量中的最小值，存入一副和原始图像大小相同的灰度图中，然后再对这幅灰度图按窗口进行最小值滤波，窗口大小与滤波半径的关系为WindowSize = 2*Radius + 1;
暗通道的先验理论指出暗通道的值是趋近于0的，即：
实际生活中造成暗原色中低通道值主要有三个因素：a)汽车、建筑物和城市中玻璃窗户的阴影，或者是树叶、树与岩石等自然景观的投影；b)色彩鲜艳的物体或表面，在RGB的三个通道中有些通道的值很低（比如绿色的草地／树／植物，红色或黄色的花朵／叶子，或者蓝色的水面）；c)颜色较暗的物体或者表面，例如灰暗色的树干和石头。总之，我们白天看到的清晰的图像中总是有阴影或者彩色的物体，这些景物的暗通道值是很小的。
业界广为使用的雾图形成模型用表达式表示如下：
I(x)=J(x)t(x)+A(1-t(x)) (1)
假设在一个窗口内的透射率是固定值，用表示，则对(7)在一个窗口求最小值
其中，I(x)就是已有的带雾图像，J(x)是期望的无雾图像，A是全球大气光成分， t(x)为透射率。现在的已知条件就是I(x)，要求目标值J(x)。如果知道A和t(x)就能求出J(x)，还好我们有暗通道先验。
由(1)得到
考虑到一个像素点由RGB三个通道组成，各个通道都需要恢复，用c表示RGB各通道，
假设在一个窗口内的透射率是固定值，用表示，则对(7)在一个窗口求最小值
上式中，J(y)是待求的无雾的图像，根据前述的暗原色先验理论有：
因此
由(8)和(10)得到
这样就得到了在一个窗口内透射率t(x)的估计值，但是根据文章中的说法，即便在明亮的场景空气中也是存在微粒的，而且正是这些空气中的小颗粒让我们感受到景深的存在，因此对式(11)加一个修正因子ω，使去雾后的图片更接近真实，式(11)变为如下:
默认取ω=0.95，越小去雾的效果越弱
之前假设全球大气成分光强度是已知的，现在根据其暗通道图像来计算。依据文章https://www.cnblogs.com/Imageshop/p/3281703.html中的方法，这样来计算：
(1) 从暗通道图中按照亮度的大小取最大的前0.1%的像素。
(2)找出这些像素的位置，根据这些位置从原始有雾图像中取这些像素的三个分量的值相加，然后计算所有这些值的平均值，作为A值。
有了t(x)和A，就可以根据公式(1)来恢复无雾图像了，恢复的公式如下：
J(x)=(I(x) - A)/t(x) + A ,考虑到t过小会导致J偏大，因此对t做了一个保护门限
=0.1,当t<0.1时取t=
,因此最终的恢复公式为：
用上述暗通道先验理论去恢复图像效果还是很明显的，下面的图像在做最小值滤波时没有考虑右侧和下方的边缘处的处理，所以能看到右侧和下方有个明显的异常带。以下例子默认使用的滤波半径都是5
原图 暗通道图
透射率图 恢复图
原图 暗通道图
透射率图 恢复图
从上图可以看出在建筑物的边缘有一个色彩模糊的条带，像是边缘的雾没有去掉似的，这是由于采用最小值滤波获取的暗通道图像过于粗糙导致的，后来何博士还提出用导向滤波来获取更精细的透射率图，从别人实现的效果看确实效果更佳，我没有自己编码实现。
写总结时发现外面雾霾很大顺手拍了一张照片，效果如下：
原图 ω=0.95
ω=0.8 ω=0.5
在文中作者提到也可以用双边滤波替代导向滤波同样能得到较为精细的透射率图，我也尝试了一下，
效果如下：
Radius=5,SigmaS=100,SigmaR=30时的透射率图和恢复图
Radius=3,SigmaS=100,SigmaR=30时的透射率图和恢复图
Radius=5,SigmaS=100,SigmaR=30时的透射率图和恢复图
可以看出用双边滤波得到的透射率图较为精细，物体的边缘就没有了模糊的条带。
代码随便写了一下，没经过优化，如下：
/*获取暗通道图像*/
void GetDarkChannel(UINT8* Src,UINT8* DarkChannel,int Width,int Height)
{
UINT8* ImgPt = Src;
UINT8* DarkPt = DarkChannel;
int Min = 255;
for (int Y = 0; Y < Height; Y++)
{
for (int X = 0; X < Width; X++)
{
Min = *ImgPt;
if (Min > *(ImgPt + 1)) Min = *(ImgPt + 1);
if (Min > *(ImgPt + 2)) Min = *(ImgPt + 2);
*DarkPt = Min;
ImgPt += 3;
DarkPt++;
}
}
}
/*计算大气成分光*/
int GetAcomponent(UINT8* Src, UINT8* pDark, int Width, int Height)
{
int Threld = 0;
int Hist[256] = { 0 };
int sum = 0;
double Accumulator = 0.0;
int Num = 0;
for (int i = 0; i < Height; i++)
{
for (int j = 0; j < Width; j++)
{
Hist[pDark[i * Width + j]]++;
}
}
for (int k = 0; k < 256; k++)
{
sum += Hist[k];
if (sum > Width * Height * 999/1000)
{
Threld = k;
break;
}
}
for (int i = 0; i < Height; i++)
{
for (int j = 0; j < Width; j++)
{
if (pDark[i * Width + j] > Threld)
{
int offset = (i * Width + j) * 3;
Accumulator += Src[offset];
Accumulator += Src[offset+1];
Accumulator += Src[offset+2];
Num++;
}
}
}
Accumulator = Accumulator / (Num * 3);
if(Accumulator > 220)
{
Accumulator = 220;
}
return (int)Accumulator;
}
/*最小值方块滤波*/
void MinFilter(UINT8* DarkChannel, int Width, int Height, int Radius)
{
int Min = 255;
UINT8* DarkSrc = DarkChannel;
for(int i = Radius;i < Height;i+=Radius)
{
for(int j = Radius;j < Width;j+=Radius)
{
int Offset = i * Width + j;
Min = DarkSrc[Offset - Width - 1];
for(int k = i-Radius;k < i;k++)
{
for(int g = j-Radius;g < j;g++)
{
Offset = k * Width + g;
if(Min > DarkSrc[Offset])
{
Min = DarkSrc[Offset];
}
}
}
for(int k = i-Radius;k < i;k++)
{
for(int g = j-Radius;g < j;g++)
{
Offset = k * Width + g;
DarkChannel[Offset] = Min;
}
}
}
}
}
/计算透射率图/
void GetTransmissionRate(UINT8* DarkChannel, int Width, int Height,int Ac)
{
for(int i = 0;i < Height;i++)
{
for(int j = 0;j < Width;j++)
{
int Offset = i * Width + j;
int tg = 255 - 95*(DarkChannel[Offset]*255+Ac/2)/(Ac*100);
if(tg < 25)
{
tg = 25;
}
DarkChannel[Offset] = tg;
}
}
}
/*恢复去雾后的图像*/
void RecoverHazeImage(UINT8* Src,UINT8* Dst,UINT8* TransmissionRate,int Width,int Height,int Ac)
{
//J = ( I - A)/t + A
int Index = 0;
for(int Y = 0; Y < Height; Y++)
{
UINT8* ImgPt = Src + Y * Width*3;
UINT8* ImgPd = Src + Y * Width*3;
for (int X = 0; X < Width; X++)
{
int temp = (ImgPt[0]-Ac)*255/TransmissionRate[Index]+Ac;
ImgPd[0] = (unsigned char)((((unsigned short)temp | ((short)(255 - temp) >> 15)) & ~temp >> 15));
temp = (ImgPt[1]-Ac)*255/TransmissionRate[Index]+Ac;
ImgPd[1] = (unsigned char)((((unsigned short)temp | ((short)(255 - temp) >> 15)) & ~temp >> 15));
temp= (ImgPt[2]-Ac)*255/TransmissionRate[Index]+Ac;
ImgPd[2] = (unsigned char)((((unsigned short)temp | ((short)(255 - temp) >> 15)) & ~temp >> 15));
ImgPt += 3;
ImgPd += 3;
Index++;
}
}
}
/*对暗通道图像做双边滤波*/
void BilaterFilter(UINT8* DarkChannel, UINT8* Dst,int Width, int Height, int Radius,int SigmaS,int SigmaR)
{
int maxk = 0;
float gauss_space_coeff = -0.5/(SigmaS*SigmaS);
float gauss_color_coeff = -0.5/(SigmaR*SigmaR);
for(int i = -Radius;i <= Radius;i++)
{
for(int j = -Radius; j <= Radius;j++)
{
float r = sqrt((float)(i*i)+(float)(j*j));
if(r > Radius)
continue;
space_weight[maxk] = exp(r*r*gauss_space_coeff);
OfsTab[maxk++] = i*Width+j;
}
}
float sum = 0;
float wsum = 0;
UINT8* DarkSrc = DarkChannel;
for(int i = Radius;i < Height - Radius;i++)
{
for(int j = Radius;j < Width - Radius;j++)
{
int Offset = i * Width + j;
int Val0 = DarkSrc[Offset];
sum = 0.0;
wsum = 0.0;
for(int k = 0;k < maxk;k++)
{
int Val1 = DarkSrc[Offset + OfsTab[k]];
int differ = Val1 - Val0;
float w = space_weight[k]*exp((float)(differ*differ*gauss_color_coeff));
sum+=Val1*w;
wsum += w;
}
Val0 = sum/wsum;
Dst[Offset] = Val0;
}
}
}从测试的图片看去雾的效果还是有的，但是有时参数选择不合理有天空的地方还是会恶化，没有别人文章中的效果那么好