1.原理介绍
1.1 大气散射模型
为了表示雾对图像成像的影响，研究者提出了一个物理模型--大气散射模型，用其来表示雾霾等恶劣天气条件对图像造成的影响。该模型由McCartney首先提出，目前已经成为计算机视觉领域中朦胧图像生成过程的经典描述，该模型的数学表达式如下：
$$ I(x)=J(x)t(x)+A(1-t(x)) $$
其中:
- $I(x)$、$J(x)$分别表示有雾图像和对应的的无雾清晰图像，$A$表示全球大气光，$t(x)$是透射矩阵，描述的是光通过传输介质后没有被散射的部分，$t(x)$的数学表达式为：(其中，β指的是光的散射系数,d(x)指的是目标与摄象机之间的距离。)
$$ t(x)=e^{-\beta d(x)} $$
- $J(x)t(x)$叫做直接衰减项，描述的是正常的清晰图像在透射媒介后发生衰减后保留下来的部分。可以看出，清晰图像会随着成像设备与物体的距离即$d(x)$的增加发生指数性衰减。
- $A(1－t(x))$部分称为大气遮罩层，表示的是全局大气光对成像的影响。
因此，根据大气散射模型可知，要想恢复出清晰图像，需要解决两个问题：
- （1）准确的对全局大气光A进行估计求解，
- （2）准确的对透射矩阵t(x)进行求解。
由上面的大气散射模型我们可以很容易得出除雾公式如下：
$$ J(x)=\frac{I(x)-A(1-t(x))}{t(x)}=\frac{I(x)+t(x)}{t(x)}+A $$
1.2 暗通道先验理论介绍
暗通道先验理论为何凯明博士2009年在CVPR(IEEE Conference On Computer Version and Pattern Recogintion/IEEE计算视觉和模式识别会议)上所提出。该理论基于对大量无雾的户外图像的观察和统计，得出了大部分的户外无雾图像的每个局部区域都存在至少一个颜色通道的强度值很低，换言之，这个区域的各个像素点的颜色通道强度的最小值是个很小的数，其值约等于0，即在其暗通道图中表现为黑色。
何凯明博士在其发表的暗通道先验除雾的论文中提出，对于任意的输入图像J，其暗通道Jdark可以用数学公式表达如下所示：
$$ J^{dark}(x)=\min_{c\in{r,g,b}}[\min_{y\in{\Omega(x)}}(J^c(y))] $$
式中$Ω(x)$则是一个正方形区域，该区域的中心是像素$x$，$J^c(x)$代表图像$J$的在某像素点的颜色强度。
由公式可知，暗通道实际上图像的最小灰度图经过最小值滤波得到的。分别对无雾图像和有雾图像进行暗通道图的求解效果如图所示:
1.3通过暗通道先验对大气光A和透射率t(x)进行估计
由大气散射模型可知，若要对图像进行除雾，等价于现在的已经知道的$I(x)$，要求解出$J(x)$，显然这个方程有无数个解。因此就需要用到暗通道先验理论对透射率$t(x)$和大气光$A$进行估计了。首先完成的是对大气光$A$的估计，根据对何凯明博士发表的论文的理解，本文对大气光求解的具体步骤如下：
- 将暗通道图和原图转化为[图片像素数量*通道数]的向量。
- 求出暗通道图对于的向量中亮度前1/1000的向量索引。
- 根据索引在原图转化得到的向量中寻找具有最高亮度的点的值，这就是A的估计值了。
在完成大气$A$的求解后，便可以对$t(x)$的求解进行推导了。参考相关资料对$t(x)$的求解过程进行推导如下:
首先将大气散射模型稍作变形得到如下式子：
$$ \frac{I^c(x)}{A^c}=t(x)\frac{J^c(x)}{A^c}+1-t(x) $$
为了对$t(x)$进行求解，先假设在每一个窗口内透射率$t(x)$为常数$\hat t(x)$，因为上面已经得到了$A$的估计值，所以在这里只需将$A$当做一个常量即可，对上式两边进行两次最小值运算，可以得到下式(i)：
$$ \min_{y \in \Omega(x)}(\min_{c \in \{r,g,b\}}(\frac{I^c(x)}{A^c})) =\hat t(x)\min_{y \in \Omega(x)}(\min_{c \in \{r,g,b\}}(\frac{J^c(x)}{A^c})) +1-\hat t(x) $$
根据暗通道先验理论有：
$$ J^{dark}(x)=min_{c\in{r,g,b}}[min_{y\in{\Omega(x)}}(J^c(y))] =0 $$
因此，可以得出：
$$ \min_{y \in \Omega(x)}(\min_{c \in \{r,g,b\}}(\frac{J^c(x)}{A^c})) =0 $$
代入式(i)可得式(ii)：
$$ \hat t(x) = 1 - \min_{y \in \Omega(x)}(\min_{c \in \{r,g,b\}}(\frac{J^c(x)}{A^c})) $$
到这就得到了透射图$t(x)$的估计结果。但是如果直接使用这一结果进行除雾效果有时产生的结果并不理想。因为空中总是会存在某些漂浮的小颗粒的，这使得我们看远处的东西会或多或少的受到雾的影响，此外，雾的存在还能让人更好的感受到物体的远近，因此，在去雾时需要考虑对一部分的雾进行保留。对雾的保留可以通过在式(ii)中引入一个在[0,1]之间的因子来实现，变化后得到的新的公式如下所示：
$$ \hat t(x) = 1 - \omega\min_{y \in \Omega(x)}(min_{c \in \{r,g,b\}}(\frac{J^c(x)}{A^c})) $$
根据何凯明博士的论文及进行了几次测试，最终本文选定的$\omega$的值为0.95。本文对复现实验中对$t(x)$进行求解过程及效果图如下图所示：
1.4 根据大气光A和透射率t(x)的估计结果进行图像去雾
到这里，我们就可以根据雾天图片退化模型进行清晰图像的恢复了。又因为当像素点x对应的透射图的$t(x)$很小时，会导致根据公式
$$ J(x)=\frac{I(x)+t(x)}{t(x)}+A $$
求解出的$J(x)$的值偏大，从而导致最终的得到的图片某些地方过曝，所以需要对其进行限制，本文选择限制的最大取值为$t_0=0.1$。从而的到最终使用的图像去雾公式为：
$$ J(x)=\frac{I(x)+t(x)}{max(t(x),t_0)}+A $$
2.代码实现
2.0绘图辅助函数
import matplotlib.pyplot as plt
#显示单个图
def show_img_by_plt(img,title):
plt.figure(figsize=(20,10)) #初始化画布
plt.axis('off') # 关掉坐标轴
plt.imshow(img) #显示图片
plt.title(title,y=-0.1) # 设置图像title
plt.show() #show
#通过subplot同时显示2个图
def show_double_img_by_subplot(img1,title1,img2,title2):
plt.figure(figsize=(20,10)) #初始化画布
#第1个位置创建一个小图
plt.subplot(1,2,1)#表示将整个图像窗口分为1行2列, 当前位置为1
plt.axis('off') # 关掉坐标轴
plt.title(title1,y=-0.1)
plt.imshow(img1)
#第2个位置创建一个小图
plt.subplot(1,2,2)#表示将整个图像窗口分为1行2列, 当前位置为2
plt.axis('off') # 关掉坐标轴
plt.title(title2,y=-0.1)
plt.imshow(img2)
plt.show() #show
#通过subplot同时显示2个图
#同时显示三个图
def show_three_img_by_subplot(img1,title1,img2,title2,img3,title3):
plt.figure(figsize=(16,8)) #初始化画布
#第1个位置创建一个小图
plt.subplot(1,3,1)#表示将整个图像窗口分为1行2列, 当前位置为1
plt.axis('off') # 关掉坐标轴
plt.title(title1,y=-0.1)
plt.imshow(img1)
#第2个位置创建一个小图
plt.subplot(1,3,2)#表示将整个图像窗口分为1行2列, 当前位置为2
plt.axis('off') # 关掉坐标轴
plt.title(title2,y=-0.1)
plt.imshow(img2)
#第3个位置创建一个小图
plt.subplot(1,3,3)#表示将整个图像窗口分为1行2列, 当前位置为2
plt.axis('off') # 关掉坐标轴
plt.title(title3,y=-0.1)
plt.imshow(img3)
plt.show() #show2.1 获取图片暗通道图
# 核心函数
#暗通道实际上是在rgb三个通道中取最小值组成灰度图，然后再进行一个最小值滤波得到的。
def get_dark_channel_img(original_img, r=15):
#原图rgb三个通道中取最小值组成灰度图
temp_img = np.min(original_img,2)
#再进行一个最小值滤波
#最小值滤波用腐蚀来替代了，其实腐蚀就是最小值滤波，最大值滤波是膨胀
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (r,r))
"""
cv2.getStructuringElement( ) 返回指定形状和尺寸的结构元素。
这个函数的第一个参数表示内核的形状，有三种形状可以选择。
矩形：MORPH_RECT;
交叉形：MORPH_CROSS;
椭圆形：MORPH_ELLIPSE;
(r,r)表示kernel的size
"""
dst_img = cv2.erode(temp_img, kernel)
"""
dst = cv.dilate(temp_img, kerne) 对图片进行腐蚀
参数说明：
dst_img 输出与输入相同大小和类型的图像.
kernel 用于侵蚀的结构元素可以使用getStructuringElement来创建结构元素。
"""
return dst_img #返回最终暗通道图#调用测试
import imageio
import numpy as np
import cv2
img_src="./image/get_dark_channel_img_clear.jpg"
original_img = imageio.imread(img_src)
dark_channel_img=get_dark_channel_img(original_img)
show_double_img_by_subplot(original_img,"original_img",dark_channel_img,"dark_channel_img")
2.2根据暗通道图和原图估计大气光
#核心函数
"""
1.从暗通道图按照亮度的大小取前0.1%的像素。
2.在这些位置中，在原始有雾图像I中寻找对应的具有最高亮度的点的值，作为A值。
"""
def estimate_A(original_img,dark_channel_img):
#计算图片的像素总数
img_h,img_w,img_c_num = original_img.shape
img_size = img_h*img_w
#计算暗通道图的像素值较大的0.1%的像素的数量
num_px = int(max(math.floor(img_size/1000),1))
#将暗通道图变为列向量
dark_vec = dark_channel_img.reshape(img_size,1);
#将图片变为列向量
img_vec = original_img.reshape(img_size,3);
#将暗通道图对应的列向量进行排序并返回从小到大的数组索引--argsort函数返回的是数组值从小到大的索引值
indices = np.argsort(dark_vec, axis=0)
#取暗通道图像素较大的0.1%的像素的索引
indices = indices[img_size-num_px::]
#将原图对应的 暗通道图的像素值排前0.1% 的像素点的像素值进行求平均
A = np.max(img_vec[indices])
return A#调用测试
import imageio
import numpy as np
import cv2
import math
# img_src=str(input("请输入原图的地址："))
img_src="./image/haze_image.jpg"
original_img = imageio.imread(img_src)/255
dark_channel_img=get_dark_channel_img(original_img)
A = estimate_A(original_img,dark_channel_img)2.3根据原图和大气光A进行投射图t(x)的估计
#核心函数
#根据论文中透射图的估计公式进行书写的
def estimate_transmission(original_img,A, omega = 0.95):
min_min_Ic_div_Ac = np.zeros(original_img.shape,"float64");
for index in range(0,3):
min_min_Ic_div_Ac[:,:,index] = original_img[:,:,index]/A
transmission = 1 - omega*get_dark_channel_img(min_min_Ic_div_Ac);
return transmission#调用测试
import imageio
import numpy as np
import cv2
# img_src=str(input("请输入原图的地址："))
img_src="./image/haze_image.jpg"
original_img = imageio.imread(img_src)
dark_channel_img=get_dark_channel_img(original_img)
A = estimate_A(original_img,dark_channel_img)
transmission = estimate_transmission(original_img,A)
show_double_img_by_subplot(original_img,"original_img",transmission,"transmission")
2.4根据估计好的大气光和透射图进行图片恢复
#核心函数
def recover_haze_by_dark_channel_prior(haze_img,t_max = 0.1):
dark_channel_img=get_dark_channel_img(haze_img)
A = estimate_A(haze_img,dark_channel_img)
transmission = estimate_transmission(haze_img,A)
clear_img = np.empty(haze_img.shape,haze_img.dtype);
for index in range(0,3):
clear_img[:,:,index] = (haze_img[:,:,index]-A[index])/cv2.max(transmission,t_max) + A[index]
return clear_img#调用测试
import imageio
import numpy as np
# haze_img_src=str(input("请输入有雾图片的地址："))
haze_img_src="./image/haze_image1.jpg"
haze_img = imageio.imread(haze_img_src)
clear_img = recover_haze_by_dark_channel_prior(haze_img)
show_double_img_by_subplot(haze_img,"haze_img",clear_img,"clear_img")
transmission = estimate_transmission(haze_img,A)
show_three_img_by_subplot(haze_img,"haze_img",clear_img,"clear_img",transmission,"transmission")
3.暗通道先验图片去雾汇总
import imageio
import numpy as np
import cv2
import math
def get_dark_channel_img(original_img, r=6):
#原图rgb三个通道中取最小值组成灰度图
temp_img = np.min(original_img,2)
#再进行一个最小值滤波
#最小值滤波用腐蚀来替代了，其实腐蚀就是最小值滤波，最大值滤波是膨胀
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (r,r))
dst_img = cv2.erode(temp_img, kernel)
return dst_img #返回最终暗通道图
def estimate_A(original_img,dark_channel_img):
#计算图片的像素总数
img_h,img_w,img_c_num = original_img.shape
img_size = img_h*img_w
#计算暗通道图的像素值较大的0.1%的像素的数量
num_px = int(max(math.floor(img_size/1000),1))
#将暗通道图变为列向量
dark_vec = dark_channel_img.reshape(img_size,1);
#将图片变为列向量
img_vec = original_img.reshape(img_size,3);
#将暗通道图对应的列向量进行排序并返回从小到大的数组索引--argsort函数返回的是数组值从小到大的索引值
indices = np.argsort(dark_vec, axis=0)
#取暗通道图像素较大的0.1%的像素的索引
indices = indices[img_size-num_px::]
#将原图对应的 暗通道图的像素值排前0.1% 的像素点的像素值进行求平均
atm_sum = np.zeros([1,3])
for index in range(1,num_px):
atm_sum = atm_sum + img_vec[indices[index]]
A = atm_sum / num_px;
return A[0]
def estimate_transmission(original_img,A, omega = 0.95):
min_min_Ic_div_Ac = np.zeros(original_img.shape,"float64");
for index in range(0,3):
min_min_Ic_div_Ac[:,:,index] = original_img[:,:,index]/A[index]
transmission = 1 - omega*get_dark_channel_img(min_min_Ic_div_Ac);
return transmission
def recover_haze_by_dark_channel_prior(haze_img,t_max = 0.1):
dark_channel_img=get_dark_channel_img(haze_img)
A = estimate_A(haze_img,dark_channel_img)
transmission = estimate_transmission(haze_img,A)
clear_img = np.empty(haze_img.shape,haze_img.dtype);
for index in range(0,3):
clear_img[:,:,index] = (haze_img[:,:,index]-A[index])/cv2.max(transmission,t_max) + A[index]
return clear_img
haze_img_src="./image/haze_image.jpg"
haze_img = imageio.imread(haze_img_src)
clear_img = recover_haze_by_dark_channel_prior(haze_img)
show_double_img_by_subplot(haze_img,"haze_img",clear_img,"clear_img")
- 应用与视频
from moviepy.editor import VideoFileClip
def process_image(image):
return deHaze(image/255.0)*255
output_file = './vedio/clear_dark_channel_short.mp4'
test_clip = VideoFileClip("./vedio/haze_short.mp4")
new_clip = test_clip.fl_image(process_image)
new_clip.write_videofile(output_file, audio=False)参考资料
- 暗通道先验原理——DCP去雾算法
- He K , Sun J , Fellow, et al. Single Image Haze Removal Using Dark Channel Prior[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2011, 33(12):2341-2353.
评论测试