本文为基于何恺明博士的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计算

通常雾最浓地方的颜色被作为大气光的估计值，可以通过暗通道来确定雾气最浓的地方，即暗通道最亮的区域为雾气最浓的地方，此时，大气光为唯一的光源。
首先选择暗通道中最亮的0.1%像素。这些像素通常是最不透明的。在这些像素中，选择输入图像I中强度最高的像素作为大气光。

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 导向滤波

以还未进行最小值滤波的暗通道图（用暗通道图效果会很差）作为导向图对t(x)进行导向滤波，细化传输。
输入图像必须得先进行归一化后再处理，否则滤波后会有超过255的值

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