本文为基于何恺明博士的Single Image Haze Removal Using Dark Channel Prior和Guided Image Filtering两篇论文的去雾算法python代码实现。

1 一些基本的定义

1.1 雾图成像模型

I(x)为原图，J(x)为无雾图像，A是大气光成分，为一常数。t(x)为透光率。
其含义就是图像I(x)为事物反射的光经过雾气衰减后加上雾气反射的大气光的结合所成的像。

1.2 暗通道定义

1.3 计算投射图t(x)

这里是一个作者统计了大量图片得出的结论，即无雾图片的暗通道是接近于0的。


实际上，即使在晴天，大气中也不是完全没有任何粒子。所以当我们看远处的物体时，雾仍然存在。 所以，我们可以通过引入一个常量参数w(0<w<1)论文取值为0.95，为远处的物体保留少量的雾。
这样t(x)会比原来大，大气光的权重会小点，使得图片整体更暗？

暗通道先验对于天空区域来说不是一个好的先验。幸运的是，在朦胧的图像I中，天空的颜色通常与大气光A非常相似，即下式趋近1。所以，在这些位置t(x)趋于0。因为天空是无限遥远的，它的传输确实接近于零。所以我们不需要预先分离天空区域。

1.4 导向滤波

详见我的另一篇文章导向滤波与opencv python实现。

2 代码

2.1 获得暗通道

def get_min_channel(img): return np.min(img,axis=2) 

最小值滤波器，其实与腐蚀一样：

讯享网def min_filter(img,r): kernel = np.ones((2*r-1,2*r-1)) return cv2.erode(img,kernel)#最小值滤波器，可用腐蚀替代 

2.2 大气光照A计算

def get_A(img_haze,dark_channel,bins_l): hist,bins = np.histogram(dark_channel,bins=bins_l)#得到直方图 d = np.cumsum(hist)/float(dark_channel.size)#累加 # print(bins) threshold=0 for i in range(bins_l-1,0,-1): if d[i]<=0.999: threshold=i break A = img_haze[dark_channel>=bins[threshold]].max() #候选区域可视化 show = np.copy(img_haze) show[dark_channel>=bins[threshold]] = 0,0,255 cv2.imwrite('./most_haze_opaque_region.jpg',show*255) return A 

暗通道最亮的区域如下图红**域处：

2.3 计算t(x)

讯享网def get_t(img_haze,A,t0=0.1,w=0.95): out = get_min_channel(img_haze) out = min_filter(out,r=7) t = 1-w*out/A #需要乘上一系数w，为远处的物体保留少量的雾 t = np.clip(t,t0,1)#论文4.4所提到t(x)趋于0容易产生噪声，所以设置一最小值0.1 return t 

这时候得出来的结果如下图：

因为没有经过导向滤波细化传输，所以效果有点差。

2.4 导向滤波

def guided_filter(I,p,win_size,eps): mean_I = cv2.blur(I,(win_size,win_size)) mean_p = cv2.blur(p,(win_size,win_size)) corr_I = cv2.blur(I*I,(win_size,win_size)) corr_Ip = cv2.blur(I*p,(win_size,win_size)) var_I = corr_I-mean_I*mean_I cov_Ip = corr_Ip - mean_I*mean_p a = cov_Ip/(var_I+eps) b = mean_p-a*mean_I mean_a = cv2.blur(a,(win_size,win_size)) mean_b = cv2.blur(b,(win_size,win_size)) q = mean_a*I + mean_b return q 

得到的结果如下图所示：

与原图相比：

可见去雾效果还是不错的。

2.5 评估

PSNR

def PSNR(target,ref): #必须归一化 target=target/255.0 ref=ref/255.0 MSE = np.mean((target-ref)2) if MSE<1e-10: return 100 MAXI=1 PSNR = 20*math.log10(MAXI/math.sqrt(MSE)) return PSNR 

SSIM

from skimage.metrics import structural_similarity as sk_cpt_ssim ssim = sk_cpt_ssim(J,I*255, win_size=11, data_range=255, multichannel=True) 

2.6 完整代码

# -*- coding: utf-8 -*- # @Time : 2022/10/1 23:08 # @Author : shuoshuo # @File : main.py # @Project : 去雾 import cv2 import numpy as np import matplotlib.pyplot as plt import math from skimage.metrics import structural_similarity as sk_cpt_ssim def guided_filter(I,p,win_size,eps): mean_I = cv2.blur(I,(win_size,win_size)) mean_p = cv2.blur(p,(win_size,win_size)) corr_I = cv2.blur(I*I,(win_size,win_size)) corr_Ip = cv2.blur(I*p,(win_size,win_size)) var_I = corr_I-mean_I*mean_I cov_Ip = corr_Ip - mean_I*mean_p a = cov_Ip/(var_I+eps) b = mean_p-a*mean_I mean_a = cv2.blur(a,(win_size,win_size)) mean_b = cv2.blur(b,(win_size,win_size)) q = mean_a*I + mean_b return q def get_min_channel(img): return np.min(img,axis=2) def min_filter(img,r): kernel = np.ones((2*r-1,2*r-1)) return cv2.erode(img,kernel)#最小值滤波器，可用腐蚀替代 def get_A(img_haze,dark_channel,bins_l): hist,bins = np.histogram(dark_channel,bins=bins_l)#得到直方图 d = np.cumsum(hist)/float(dark_channel.size)#累加 # print(bins) threshold=0 for i in range(bins_l-1,0,-1): if d[i]<=0.999: threshold=i break A = img_haze[dark_channel>=bins[threshold]].max() #候选区域可视化 show = np.copy(img_haze) show[dark_channel>=bins[threshold]] = 0,0,255 cv2.imwrite('./most_haze_opaque_region.jpg',show*255) return A def get_t(img_haze,A,t0=0.1,w=0.95): out = get_min_channel(img_haze) out = min_filter(out,r=7) t = 1-w*out/A #需要乘上一系数w，为远处的物体保留少量的雾 t = np.clip(t,t0,1)#论文4.4所提到t(x)趋于0容易产生噪声，所以设置一最小值0.1 return t def PSNR(target,ref): #必须归一化 target=target/255.0 ref=ref/255.0 MSE = np.mean((target-ref)2) if MSE<1e-10: return 100 MAXI=1 PSNR = 20*math.log10(MAXI/math.sqrt(MSE)) return PSNR if __name__ == '__main__': I = cv2.imread('test.jpg')/255.0 dark_channel = get_min_channel(I) dark_channel_1 = min_filter(dark_channel,r=7) # cv2.imwrite("./dark_channel.jpg", dark_channel_1*255) A = get_A(I,dark_channel_1,bins_l=2000) t = get_t(I,A) t = guided_filter(dark_channel,t,81,0.001) t = t[:,:,np.newaxis].repeat(3,axis=2)#升维至(r,w,3) J = (I-A)/t +A J = np.clip(J,0,1) J = J*255 J =np.uint8(J) cv2.imwrite("./result.jpg",J) #评估 PSNR = PSNR(J,I*255) print(f"PSNR:{ 
     PSNR}") ssim = sk_cpt_ssim(J,I*255, win_size=11, data_range=255, multichannel=True) print(f"ssim:{ 
     ssim}") 

3 参考文献

何恺明的两篇论文：
Single Image Haze Removal Using Dark Channel Prior
Guided Image Filtering
论文链接
Dehazing for Image and Video Using Guided Filter