论文阅读：Fast Single Image Dehazing Using Saturation Based Transmission Map Estimation

代码在文末给出，因为只有一个py文件。该代码是作者给出的。

1. 提出的方法

1.1 大气散射模型

    大气散射模型：
$$\mathbf{H}(x)=t(x)\mathbf{J}(x)+[1-t(x)]\mathbf{A}\text{\hspace{1em}(1)}$$

    移项整理后：
$$\mathbf{J}(x)=D(\mathbf{H}(x),\mathbf{A},t(x))=\frac{\mathbf{H}(x)-\mathbf{A}}{t(x)}+\mathbf{A}\text{\hspace{1em}(2)}$$

    这里雾图使用$\mathbf{H}(x)$而不是$\mathbf{I}(x)$，因为后面会用到$I$的符号。

1.2 传输率图公式

    给定$\mathbf{A}$的情况下，用$\mathbf{A}$归一化$(1)$得：
H c ( x ) A c = t ( x ) J c ( x ) A c + 1 − t ( x ) , c ∈ { r , g , b } (3) \frac{H^{c}(x)}{A^{c}}=t(x) \frac{J^{c}(x)}{A^{c}}+1-t(x), c \in\{r, g, b\} \tag{3} AcHc(x)​=t(x)AcJc(x)​+1−t(x),c∈{r,g,b}(3)

    两边对通道求最小值操作，得到：
$$\underset{c}{\min }{H}_n^c=t(x)\underset{c}{\min }{J}_n^c+1-t(x)\text{\hspace{1em}(4)}$$

    其中，$H_n^c(x)=\frac{H^c(x)}{A^c}$和$J_n^c(x)=\frac{J^c(x)}{A^c}$。移项整理后得：
$$t(x)=\frac{1-m_{{\mathbf{H}}_n}(x)}{1-m_{{\mathbf{J}}_n}(x)}\text{\hspace{1em}(5)}$$

    其中，$m_{{\mathbf{H}}_n}(x)=\underset{c}{\min }{H}_n^c(x)$，$m_{{\mathbf{J}}_n}(x)=\underset{c}{\min }{J}_n^c(x)$。下图展示了本文提出的$m_{{\mathbf{J}}_n}(x)$和暗通道先验中，清晰图像用$\mathbf{A}$归一化的结果。

图1 大气光归一化结果

1.3 传输率图估计

    定义饱和度公式：给定一幅图像$\mathbf{K}(x)$，像素点$x$的饱和度定义为$S_{\mathbf{K}}(x)$。公式如下：
$$S_{\mathbf{K}}(x)=1-\frac{m_{\mathbf{K}}(x)}{I_{\mathbf{K}}(x)}\text{\hspace{1em}(6)}$$

    其中， m K ( x ) = min ⁡ c ∈ { r , g , b } K c ( x ) m_{\mathbf{K}}(x)=\min\limits_{c \in\{r, g, b\}} K^{c}(x) mK​(x)=c∈{r,g,b}min​Kc(x)，$I_{\mathbf{K}}(x)$为$\mathbf{K}(x)$的强度，定义为：
$$I_{\mathbf{K}}(x)=\frac{K^r(x)+K^g(x)+K^b(x)}{3}\text{\hspace{1em}(7)}$$

    $(6)$移项整理得到$m_{\mathbf{K}}(x)=I_{\mathbf{K}}(x)\left(1-S_{\mathbf{K}}(x)\right)$，用该形式重写$(5)$得：
$$t(x)=\frac{1-I_{{\mathbf{H}}_n}(x)\left(1-S_{{\mathbf{H}}_n}(x)\right)}{1-I_{{\mathbf{J}}_n}(x)\left(1-S_{{\mathbf{J}}_n}(x)\right)}\text{\hspace{1em}(8)}$$

    由上可知，传输率图可以通过归一化场景辐射的强度和饱和度来确定。两个未知数（$I_{{\mathbf{J}}_n}(x)$和$S_{{\mathbf{J}}_n}(x)$）用于估计$t(x)$。
    用$\mathbf{A}$归一化$(2)$，加入强度定义$(7)$得到：
$$I_{{\mathbf{J}}_n}(x)=\frac{I_{{\mathbf{H}}_n}(x)-1}{t(x)}+1\text{\hspace{1em}(9)}$$

$$\mathbf{J}(x)=\frac{\mathbf{H}(x)-\mathbf{A}}{t(x)}+\mathbf{A}\text{\hspace{1em}(2)}$$

    用$\mathbf{A}$归一化：
$$\frac{J^c(x)}{A^c}=\frac{\frac{H^c(x)}{A^c}-1}{t(x)}+1$$

    加入定义$J_n^c(x)=\frac{J^c(x)}{A^c}$，$H_n^c(x)=\frac{H^c(x)}{A^c}$得：
$$J_n^c(x)=\frac{H_n^c(x)-1}{t(x)}+1$$

    三个通道：
$$J_n^r(x)=\frac{H_n^r(x)-1}{t(x)}+1$$

$$J_n^g(x)=\frac{H_n^g(x)-1}{t(x)}+1$$

$$J_n^b(x)=\frac{H_n^b(x)-1}{t(x)}+1$$

    上面那三个式子相加，加入$(7)$的定义得：
$$3I_{{\mathbf{J}}_n}(x)=\frac{3I_{{\mathbf{H}}_n}(x)-3}{t(x)}+3$$

    化简得$(9)$。

    将$(9)$带入$(8)$中，得：
$$t(x)=\frac{1-I_{{\mathbf{H}}_n}(x)\left(1-S_{{\mathbf{H}}_n}(x)\right)}{1-(\frac{I_{{\mathbf{H}}_n}(x)-1}{t(x)}+1)\left(1-S_{{\mathbf{J}}_n}(x)\right)}\text{\hspace{1em}(10)}$$

    化简得：
$$t(x)=1-I_{{\mathbf{H}}_n}(x)\left(1-\frac{S_{{\mathbf{H}}_n}(x)}{S_{{\mathbf{J}}_n}(x)}\right)\text{\hspace{1em}(11)}$$

    $(10)$移项：
$$t(x)-(I_{{\mathbf{H}}_n}(x)-1+t(x))(1-S_{{\mathbf{J}}_n}(x))=1-I_{{\mathbf{H}}_n}(x)(1-S_{{\mathbf{H}}_n}(x))$$

$$t(x)(1-(1-S_{{\mathbf{J}}_n}(x)))-(I_{{\mathbf{H}}_n}(x)-1)(1-S_{{\mathbf{J}}_n}(x))=1-I_{{\mathbf{H}}_n}(x)(1-S_{{\mathbf{H}}_n}(x))$$

$$\begin{array}{rl}t(x)S_{{\mathbf{J}}_n}(x)=& 1-I_{{\mathbf{H}}_n}(x)(1-S_{{\mathbf{H}}_n}(x))+I_{{\mathbf{H}}_n}(x)(1-S_{{\mathbf{J}}_n}(x))-(1-S_{{\mathbf{J}}_n}(x))\\ =& S_{{\mathbf{J}}_n}(x)-I_{{\mathbf{H}}_n}(S_{{\mathbf{J}}_n}(x)-S_{{\mathbf{H}}_n}(x))\end{array}$$

    两边乘以$S_{{\mathbf{J}}_n}(x)$得$(11)$。

    有了$(11)$，未知单元就只剩下$S_{{\mathbf{J}}_n}(x)$。估计$S_{{\mathbf{J}}_n}(x)$比估计$m_{\mathbf{J}}(x)$要简单得多。

1.4 估计场景辐射得饱和度

    $(11)$移项变形后得：
$$S_{{\mathbf{H}}_n}(x)=t(x)\frac{I_{{\mathbf{J}}_n}(x)}{I_{{\mathbf{H}}_n}(x)}{S}_{{\mathbf{J}}_n}(x)\text{\hspace{1em}(12)}$$

    $I_{{\mathbf{H}}_n}(x)$和$I_{{\mathbf{J}}_n}(x)$的关系通过$(9)$可得：
$$I_{{\mathbf{H}}_n}(x)=t(x)I_{{\mathbf{J}}_n}(x)+1-t(x)\text{\hspace{1em}(13)}$$

    因为$0\le t(x)\le 1$，从$(12)(13)$可得$I_{{\mathbf{J}}_n}(x)\le I_{{\mathbf{H}}_n}(x)$和$S_{{\mathbf{J}}_n}(x)\ge S_{{\mathbf{H}}_n}(x)$。

    由$(13)$和$0\le t(x)\le 1$得：
$$t(x)=\frac{I_{{\mathbf{H}}_n}(x)-1}{I_{{\mathbf{J}}_n}(x)-1}\le 1$$

    因为$I_{{\mathbf{H}}_n}(x)$和$I_{{\mathbf{J}}_n}(x)$都已归一化，即$0\le I_{{\mathbf{H}}_n}(x),I_{{\mathbf{J}}_n}(x)\le 1$。所以$I_{{\mathbf{J}}_n}(x)\le I_{{\mathbf{H}}_n}(x)$，结合$(12)$得：$S_{{\mathbf{J}}_n}(x)\ge S_{{\mathbf{H}}_n}(x)$。

    通过以上推导，可知雾图的饱和度小于清晰图像的饱和度。本文通过将雾图的饱和度进行拉伸得到清晰图像的饱和度。文中给出了三种饱和度拉伸的方法，分别如下：

1. 拉伸函数1
   $$\tau (z,p)=\{\begin{array}{ll}0.5\left(1-{\left(1-\frac{z}{0.5}\right)}^p\right),& z\le 0.5\\ 0.5+0.5{\left(\frac{z-0.5}{0.5}\right)}^2,& 0.5<z\le 1\end{array}\text{\hspace{1em}(14)}$$

    其中，$\tau (z,p)$是关于$z$自适应的拉伸函数，$p(p>0)$为功率常数。

2. 拉伸函数2
   $$q(z)=z(2-z)\text{\hspace{1em}(15)}$$

    在$(15)$中，$q(z)$是$z$的二次函数，在$z=1$时最大值为$1$。

3. 拉伸函数3
   $$\eta (z,\gamma )=\frac{\left(z^{\frac{1}{\gamma }}+\left(1-(1-z{)}^{\frac{1}{\gamma }}\right)\right)}{2}\text{\hspace{1em}(16)}$$

    $\gamma$为决定$\eta (z,\gamma )$形状的常数。
    图2显示了三个拉伸函数的示例。图2(a)展示了这些拉伸函数的曲线图。在图2(b)中，展示了一幅雾图及其灰度图，以及三种对灰度图不同拉伸函数的结果图。

图2 三种拉伸函数

    图3展示了使用各种拉伸函数的去雾图像。如图3所示，去雾结果在很大程度上不取决于拉伸功能。

图3 三种拉伸函数的去雾结果

1.5 处理颜色偏色

     通常，假定三个通道大气光值都是相似的，所以雾几乎没有颜色成分，假定所有大气光都是相似的。但是，黄色的粉尘或某些照明条件下可能会使雾带有颜色偏色。图5展示了颜色偏色的情况。颜色偏色的情况下，去雾效果并不理想。

图4 大气光偏色情况

    本文提出了一种使用白平衡方法的简单的偏色去除算法。 $\mathbf{H}(x)$白平衡处理后的结果${\mathbf{H}}_{WB}(x)$通过下式获得：
$$H_{WB}^c(x)=W\left(H^c(x)\right)=\frac{\mu \left(H^g(x)\right)}{\mu \left(H^c(x)\right)}{H}^c(x)\text{\hspace{1em}(17)}$$

    其中，$\mu (H^c(x))$为$H^c(x)$的均值，$W(\cdot )$为白平衡函数。
    经过白平衡处理后，得到白平衡雾图${\mathbf{H}}_{WB}(x)$。该图的大气光值为${\mathbf{A}}_{WB}$。如果$\mathbf{A}$和${\mathbf{A}}_{WB}$有较大偏差，说明雾是有偏色的。因此，比较$\mathbf{A}$和${\mathbf{A}}_{WB}$的偏差，决定是否使用白平衡后的雾图来去雾。
    定义$\Delta (\mathbf{y})$为$\mathbf{y}$的最大值和最小值的差值。分两种情况：

1. $\Delta (\mathbf{A})\le \Delta \left({\mathbf{A}}_{WB}\right)$
       此时，上述白平衡处理会增加$\mathbf{A}$最大值和最小值的差值。所以不使用白平衡处理，直接使用$(2)$得到去雾图像，即$\mathbf{J}(\mathrm{x})=D(\mathbf{H}(\mathrm{x}),\mathbf{A},\mathrm{t}(\mathrm{x}))$。

2. $\Delta (\mathbf{A})>\Delta \left({\mathbf{A}}_{WB}\right)$
       此时，由于灰度世界假设，上述白平衡处理会减少$\mathbf{A}$最大值和最小值的差值。此时使用如下公式去雾：$\mathbf{J}(\mathrm{x})=W\left[D\left(\mathbf{H}(\mathrm{x}),{\mathbf{A}}_{WB},t(\mathrm{x})\right)\right]$。

    即，对于一幅雾图$\mathbf{H}(x)$，其去雾结果通过下面式子获得：
$$\mathbf{J}(x)=\{\begin{array}{ll}D(\mathbf{H}(x),\mathbf{A},t(x)),& \Delta (\mathbf{A})\le \Delta \left({\mathbf{A}}_{WB}\right)\\ W\left[D\left(\mathbf{H}(x),{\mathbf{A}}_{WB},t(x)\right)\right],& \text{ otherwise }\end{array}\text{\hspace{1em}(18)}$$

    如下图5所示。对于(a)，$\Delta (\mathbf{A})<\Delta \left({\mathbf{A}}_{WB}\right)$，比较白平衡处理与否的结果，显然，未处理的结果更好。反之(b)，$\Delta (\mathbf{A})>\Delta \left({\mathbf{A}}_{WB}\right)$，白平衡处理的结果更好。从而验证了$(18)$的正确性。

图5 白平衡处理效果

1.6 大气光估计

    本文采用了三种不同的方法估计大气光，分别为He，Tang和quad-tree subdivision-based方法。实验结果表明，使用何种大气光估计方法，对去雾结果影响不大。如下图6所示。

图6 三种估计大气光方法对去雾结果的影响

1.7 算法总结

    由于博客的公式编号和原文不同，所以上述总结中公式编号还请查看原文。

2. 实验结果

    参考原文。

3. 读后感

    本文提出了一种传统方法，该方法简单有效。
    一句话总结：将大气散射模型通过作者定义的饱和度变换成如下形式$t(x)=1-I_{{\mathbf{H}}_n}(x)\left(1-\frac{S_{{\mathbf{H}}_n}(x)}{S_{{\mathbf{J}}_n}(x)}\right)$，再通过对$S_{{\mathbf{H}}_n}(x)$的对比度拉伸获得唯一未知量$S_{{\mathbf{J}}_n}(x)$，从而得到传输率图。

4. 代码

# -*- coding: utf-8 -*-
"""
Created on Thu Sep 19 15:09:58 2019

@author: Se Eun Kim, Tae Hee Park, and Il Kyu
Paper: Fast Single Image Dehazing Using Saturation Based Transmission Map Estimation, 
To be appear in IEEE Transaction on Image Processing, 2020
""" 

import numpy as np
import cv2
from skimage import morphology
import math
from tkinter import filedialog
from tkinter import *
from tkinter.filedialog import askopenfilename
import time

'''
   pixel value: i2f: 0-255 to 0-1, f2i: 0-1 to 0-255 
'''
def i2f(i_image):
    f_image = np.float32(i_image)/255.0
    return f_image

def f2i(f_image):
    i_image = np.uint8(f_image*255.0)
    return i_image 

'''
    Compute 'A' as described by Tang et al. (CVPR 2014)
'''
def Compute_A_Tang(im):

    # Parameters
    erosion_window = 15
    n_bins = 200

    R = im[:, :, 2]
    G = im[:, :, 1]
    B = im[:, :, 0]
    
    # compute the dark channel
    dark = morphology.erosion(np.min(im, 2), morphology.square(erosion_window))

    [h, edges] = np.histogram(dark, n_bins)
    numpixel = im.shape[0]*im.shape[1]
    thr_frac = numpixel*0.99
    csum = np.cumsum(h)
    nz_idx = np.nonzero(csum > thr_frac)[0][0]
    dc_thr = edges[nz_idx]
    mask = dark >= dc_thr
    # similar to DCP till this step
    # next, median of these top 0.1% pixels
    # median of the RGB values of the pixels in the mask
    rs = R[mask]
    gs = G[mask]
    bs = B[mask]
    
    A = np.zeros((1,3))

    A[0, 2] = np.median(rs)
    A[0, 1] = np.median(gs)
    A[0, 0] = np.median(bs)

    return A

'''
    Compute intensity: GetIntensity, and Saturation: GetSauration
'''

def GetIntensity(fi):
    return cv2.divide(fi[:, :, 0] + fi[:, :, 1] + fi[:, :, 2], 3)
    
def GetSaturation(fi, intensity):
    min_rgb = cv2.min(cv2.min( fi[:, :, 0], fi[:, :, 1]), fi[:, :, 2])
    me = np.finfo(np.float32).eps
    S = 1.0 - min_rgb/(intensity+me)
    return S 

'''
    Estimate saturation of scene radiance: 3 memthods
    
'''
def EstimateSaturation(h_saturation, p1):
    p2 = 2.0
    k1 =0.5*(1.0-cv2.pow(1.0-2.0*h_saturation, p1))
    k2 = 0.5+0.5*cv2.pow((h_saturation-0.5)/0.5, p2)
    j_saturation = np.where(h_saturation<=0.5, k1, k2)
    j_saturation = np.maximum(j_saturation, h_saturation)     

    return j_saturation      

def EstimateSaturation_Quadratic(h_saturation):
    return h_saturation*(2.0-h_saturation)   

def EstimateSaturation_Gamma(h_saturation, g):
    j_saturation = (np.power(h_saturation, 1.0/g) +1.0- np.power(1.0-h_saturation,1.0/g))/2.0
    j_saturation = np.maximum(j_saturation, h_saturation)     
    
    return j_saturation  

'''
    Estimate Transmission Map
'''
def EstimateTransimission(h_intensity, h_saturation, j_saturation):
      
    Td = h_intensity*(j_saturation-h_saturation)
    Tmn =j_saturation
    Tmap = 1.0 - (Td/Tmn)
    me = np.finfo(np.float32).eps
    Tmap = np.clip(Tmap,me, 1.0)
    cv2.imshow('Transmission Map ', f2i(Tmap))
    cv2.waitKey()
    a=1
    return Tmap 

'''
    Recover dehazed image
'''
def Recover(im, tmap, A):
 
    res = np.empty(im.shape,im.dtype) 
    
    for ind in range(0,3):
        res[:,:,ind] = (im[:,:,ind]-A[0,ind])/tmap+A[0,ind]
        res[:,:,ind] = np.clip(res[:,:,ind], 0.0, 1.0)
   
    return res

'''
    Adjust image range
'''
def Adjust(im, perh, perl):
    aim = np.empty(im.shape,im.dtype)
    temp = np.empty(im.shape,im.dtype)
    im_h = np.percentile(im, perh)
    im_l = np.percentile(im, perl)

    for ind in range(0,3):
 
        aim[:,:,ind] = (im[:,:,ind]-im_l)/(im_h-im_l)
        temp[:,:,ind] = np.clip(aim[:,:,ind], 0.0, 1.0)
        a=1
    return aim


'''                
    Normalize image 0 between 1
'''
def Normalize(im):
    aim = np.empty(im.shape,im.dtype) 

    for ind in range(0,3):
        im_h = np.max(im[:,:,ind])
        im_l = np.min(im[:,:,ind])     
        aim[:, :, ind] = (im[:, :, ind]-im_l)/(im_h-im_l)
        aim[:,:,ind] = np.clip(aim[:,:,ind], 0.0, 1.0)     
        
    return aim

'''
   White balance using grayworld assumption
'''
def gray_world(im):
    aim = np.empty(im.shape,im.dtype) 

    mu_r = np.average(im[:,:,2])
    mu_g = np.average(im[:,:,1])
    mu_b = np.average(im[:,:,0])
    aim[:,:,0] = np.minimum(im[:,:,0]*(mu_g/mu_b),1.0)
    aim[:,:,2] = np.minimum(im[:,:,2]*(mu_g/mu_r),1.0)
    aim[:,:,1] = im[:,:,1]
 
    return  aim

'''
  CLAHE
'''
def Clahe(im, clip):
  HSV = cv2.cvtColor(f2i(im), cv2.COLOR_BGR2HSV)
  clahe = cv2.createCLAHE(clipLimit=clip, tileGridSize=(8,8))
  HSV[:, :, 2] = clahe.apply(HSV[:, :, 2])
  result_im = i2f(cv2.cvtColor(HSV, cv2.COLOR_HSV2BGR))  

  return result_im
   
'''
 Main
'''
Tk().withdraw()
filename = askopenfilename(title = 'Select input hazy image')
hazy_image = i2f(cv2.imread(filename, cv2.IMREAD_COLOR))

start_time = time.time()

A = Compute_A_Tang(hazy_image)
S_A = np.max(A)-np.min(A)

'''
  Compute white balanced A
'''
hazy_imageWB = gray_world(hazy_image)
A_WB = Compute_A_Tang(hazy_imageWB)
S_AWB = np.max(A_WB)-np.min(A_WB)

'''
  Parameter set
      parameters for Adjust: perh = 99.9, perl = 0.5 
     (for full-reference image: perh = 100, perl = 0
      parameter for selecting two branch: epsilon = 0.00 - 0.1
      parameter for CLAHE: clip size  = 1 (0-2)
'''
perh = 99.9
perl = 0.5
epsilon = 0.02
cl  = 1

if S_A < S_AWB+epsilon:
   print('Phase I - Normal')    
   hazy_imagen = np.empty(hazy_image.shape,hazy_image.dtype) 

   for ind in range(0,3):
       hazy_imagen[:,:,ind]=hazy_image[:,:,ind]/A[0,ind]   

   hazy_imagen = Normalize(hazy_imagen)       
   hazy_I = GetIntensity(hazy_imagen)
   hazy_S = GetSaturation(hazy_imagen, hazy_I)         

   '''
   Stretch function I: 0.2 (heavy haze) - 0.4 (low haze)
   Stretch function II: no parameter
   Stretch function III: 4.0 (heavy haze) - 2.0 (low haze) 
   ''' 

   est_S = EstimateSaturation_Gamma(hazy_S, 0.2)
   #est_S = EstimateSaturation_Quadratic(hazy_S)
   #est_S = EstimateSaturation(hazy_S, 2.0)   
   Transmap = EstimateTransimission(hazy_I, hazy_S, est_S)
   r_image = Recover(hazy_image,Transmap, A)
   r_image = Adjust(r_image, perh, perl)
  
else: 
   print('Phase II -White Balance')    
   hazy_imagen = np.empty(hazy_image.shape,hazy_image.dtype) 
   
   for ind in range(0,3):
       hazy_imagen[:,:,ind]=hazy_image[:,:,ind]/A_WB[0,ind]
   
   hazy_imagen = Normalize(hazy_imagen)       
   hazy_I = GetIntensity(hazy_imagen)
   hazy_S = GetSaturation(hazy_imagen, hazy_I)   

   est_S = EstimateSaturation_Gamma(hazy_S, 0.2)
   #est_S = EstimateSaturation_Quadratic(hazy_S)
   #est_S = EstimateSaturation(hazy_S, 2.0)
    
   Transmap = EstimateTransimission(hazy_I, hazy_S, est_S)
   r_image = Recover(hazy_image,Transmap, A_WB)
   r_image = Adjust(r_image, perh, perl)
   r_image = gray_world(r_image)
   
result_ce = Clahe(r_image, cl)

end_time = time.time()

print("--- %s seconds ---" %(end_time - start_time))

cv2.imwrite('result.png', f2i(r_image))
cv2.imwrite('result_ce.png', f2i(result_ce))
cv2.imwrite('trmap.png', f2i(Transmap))

(h,w) = hazy_image.shape[:2]
max_size = 800
if h>=w:
    if h>max_size:
        ns = h/max_size
        nh = int(h/ns)
        nw = int(w/ns)
    else:
        nh = h
        nw = w
        
else:
    if w>max_size:
        ns = w/max_size
        nh = int(h/ns)
        nw = int(w/ns)
    else:
        nh = h
        nw = w      

hazy_r = cv2.resize(hazy_image, (nw,nh))
trans_r = cv2.resize(Transmap, (nw,nh))
r_r = cv2.resize(r_image, (nw,nh))
ce_r = cv2.resize(result_ce, (nw,nh))

cv2.imshow('Input hazy image', f2i(hazy_r))
cv2.imshow('Transmission Map ', f2i(trans_r))
cv2.imshow('first dehazed image', f2i(r_r))
cv2.imshow('enhanced dehazed image',f2i(ce_r))
cv2.waitKey()
cv2.destroyAllWindows()