暗通道先验去雾算法和引导滤波学习笔记

关于暗通道先验图像去雾的一些理论概述

本文是基于何恺明的暗通道先验去雾算法的论文Single Image Haze Removal Using Dark Channel Prior，大家可以去谷歌学术下载看一看。
https://ieeexplore.ieee.org/document/5206515

大气散射模型

在Narasimhan的Vision and the Atmosphere论文中，他假设环境光波长与大气散射系数无关，提出了一个简化的大气散射模型，这也是图像去雾领域中广泛使用的物理模型之一，其公式如下所示:
$$I(x)=J(x)t(x)+A[1-t(x)]\text{\hspace{1em}(1)}$$

上式中，$x$表示图像的每个像素点的坐标，$I(x)$表示在有雾条件下的观测图像，$J(x)$表示的是要恢复对的无雾的图像，$A$表示整个图像的全局大气光，$t(x)$表示介质透射率，代表在坐标$x$位置处所有光线穿过雾的比例，在均匀介质中，$t(x)$的表达式为
$$t(x)=e^{-\beta d(x)}\text{\hspace{1em}(2)}$$
其中$d(x)$表示场景深度，$\beta$表示大气散射系数。

图像去雾的最终目的就是在$I(x)$的基础上恢复出$J(x)$，方程(1)是一个不定方程，无法对其直接求解，所以我们必须要有先验知识。

暗通道先验理论

暗通道先验理论认为在大多数非天空的图片区域中，RGB三个颜色通道中至少有一个颜色通道在某些像素处具有非常低的强度，根据这样的先验理论，就有如下定义式：

J dark  ( x ) = min ⁡ c ∈ { r , g , b } { min ⁡ y ∈ Ω ( x ) { J c ( y ) } } (3) J^{\text {dark }}(x)=\min _{c \in\{r, g, b\}}\left\{\min _{y \in \Omega(x)}\left\{J^c(y)\right\}\right\}\tag3 Jdark (x)=c∈{r,g,b}min​{y∈Ω(x)min​{Jc(y)}}(3)

其中，$\Omega (x)$是以$x$为中心的局部像素，$J^c$是$J$的颜色通道。

利用matlab求得的有雾图像的暗通道如图1所示

透射率估计

我们对图像的局部像素应用公式(1)，假设大气光$A$是已知的，则有等式
$$\underset{y\in \Omega (x)}{\min }\left(I^c(y)\right)=\tilde{t}(x)\underset{y\in \Omega (x)}{\min }(J^c(x)+(1-\tilde{t}(x)A^c\text{\hspace{1em}(4)}$$

其中，$\tilde{t}(x)$为透射率估计值。
之后，我们对式(4)变形，再对图像的RGB三个颜色通道取最小值，可以得到：
$$\underset{c\in R,G,B}{\min }(\underset{y\in \Omega (x)}{\min }(\frac{I^c(y)}{A^c}))=\tilde{t}(x)\underset{c\in R,G,B}{\min }(\underset{y\in \Omega (x)}{\min }(\frac{J^c(y)}{A^c})+(1-\tilde{t}(x))\text{\hspace{1em}(5)}$$

显然，根据暗通道先验理论，对于无雾霾图像的$J^{d}ark$是趋向于0的，即
$$\underset{c\in R,G,B}{\min }(\underset{y\in \Omega (x)}{\min }(\frac{J^c(y)}{A^c})=0\text{\hspace{1em}(6)}$$

又因为大气光$A^c$是已知的且大于0，综合式(5)和式(6)，可以得到估计透射率的具体公式：

$$\tilde{t}(x)=1-\underset{c\in R,G,B}{\min }(\underset{y\in \Omega (x)}{\min }(\frac{I^c(y)}{A^c}))\text{\hspace{1em}(7)}$$

此外，我们还需要考虑如果彻底的除去雾霾，图像可能会变得不够真实，丢失深度信息，所以还需要引入一个常数$\alpha (0<\alpha <1)$,对远处的物体保持较少的雾气。修正后的透射率估计式如下
$$\tilde{t}(x)=1-\alpha \underset{c\in R,G,B}{\min }(\underset{y\in \Omega (x)}{\min }(\frac{I^c(y)}{A^c}))\text{\hspace{1em}(8)}$$

为了防止恢复出来的无雾图像的$t(x)$趋向于0，导致整体图像有噪声，降低图像质量，我们设$t_0$为$t(x)$的下界，本文中，代码里的$t_0$取0.1，最终，我们可以得到恢复的图像表达式:
$$J(x)=\frac{I(x)-A}{\max (t(x),t_0)}+A\text{\hspace{1em}(9)}$$

可以看出来，利用暗通道先验理论估计出的透射率在图像的边缘部分较为模糊，不能很好地展现原图像的边缘特征，所以我们还需要对透射率进行细化。

导向滤波

这里我们依然是介绍何恺明的论文Guided Image Filtering中的方法
https://link.springer.com/chapter/10.1007/978-3-642-15549-9_1

引导滤波方法也属于保持边缘的一种滤波方法，在进行滤波时，需要一幅引导图像，当引导图像为图像本身时，引导滤波就是一个保持边缘的滤波操作。

设需要处理的原图像为$P$，导向图为$I$，滤波输出为$Q$，则有

$$Q_i=\sum _j{W}_{ij}(I)P_j$$

其中，$Q_i$表示以$i$为中心的窗口进行滤波得到的输出图像,$W_{i}j$是滤波核，它是导向图$I$的函数，我们认为输出图像$Q$是引导图像$I$的线性变换\cite{he2012guided},即有
$$q_i=a_k{I}_i+b_k,\forall i\in w_k$$

这里的$a_k$和$b_k$是常数项系数，$w_k$是以$k$为中心的核函数窗口。引导滤波主要是寻求线性变换的系数$a_k$和$b_k$的最佳值，使得$P$和$Q$的差异最小化，故我们可以得到目标函数如下：
$$\min E\left(a_k,b_k\right)=\sum _{i\in w_k}\left({\left(a_k{I}_i+b_k-p_i\right)}^2+ϵa_k^2\right)$$

其中，$ϵ$是正则化参数，避免系数$a_k$过大。利用最小二乘法，我们可以求得最优的$a_k$和$b_k$如下
$$\begin{gathered} a_k=\frac{\frac{1}{|\omega |}{\sum }_{i\in {\omega }_k}{I}_i{p}_i-u_k{\bar{p}}_k}{{\sigma }_k^2+\epsilon } \\ {b}_k={\bar{p}}_k-a_k{u}_k \end{gathered}$$

其中，$|\omega |$是$w_k$内的像素总数量，$u_k$和${\sigma }_k^2$是窗口$w_k$内的像素的均值和方差。最后，我们对所有包含像素点$i$的窗口$w_k$的邻域的输出求平均值，可以得到引导滤波的结果式如下：
$$q_i=\frac{1}{|\omega |}\sum _{k,j\in w_k}(a_k{I}_i+b_k)=\overline{a_i}{I}_i+\overline{b_i}$$

去雾结果

matlab代码

本次实验的代码参考的是GitHub上的引导滤波优化去雾算法
https://github.com/sjtrny/Dark-Channel-Haze-Removal

主要函数如下

• 获取图像暗通道

function dark_channel = get_dark_channel(image, win_size)
[m, n, ~] = size(image);
pad_size = floor(win_size/2);
padded_image = padarray(image, [pad_size pad_size], Inf);
dark_channel = zeros(m, n); 
for j = 1 : m
    for i = 1 : n
        patch = padded_image(j : j + (win_size-1), i : i + (win_size-1), :);
        dark_channel(j,i) = min(patch(:));
     end
end
end

• 滤波窗

function sum_img = window_sum_filter(image, r)

% sum_img(x, y) = = sum(sum(image(x-r:x+r, y-r:y+r)));

[h, w] = size(image);
sum_img = zeros(size(image));

% Y axis
im_cum = cumsum(image, 1);

sum_img(1:r+1, :) = im_cum(1+r:2*r+1, :);
sum_img(r+2:h-r, :) = im_cum(2*r+2:h, :) - im_cum(1:h-2*r-1, :);
sum_img(h-r+1:h, :) = repmat(im_cum(h, :), [r, 1]) - im_cum(h-2*r:h-r-1, :);

% X axis
im_cum = cumsum(sum_img, 2);

sum_img(:, 1:r+1) = im_cum(:, 1+r:2*r+1);
sum_img(:, r+2:w-r) = im_cum(:, 2*r+2:w) - im_cum(:, 1:w-2*r-1);
sum_img(:, w-r+1:w) = repmat(im_cum(:, w), [1, r]) - im_cum(:, w-2*r:w-r-1);

end

• 透射率估计

function trans_est = get_transmission_estimate(image, atmosphere, omega, win_size)
[m, n, ~] = size(image);
rep_atmosphere = repmat(reshape(atmosphere, [1, 1, 3]), m, n); 
trans_est = 1 - omega * get_dark_channel( image ./ rep_atmosphere, win_size);
end

• 获得大气光
  这里的估计大气光的方法是我重新对原本的方法进行了改进，原本的方法是取暗通道中亮度最亮的10%的像素的均值作为大气光，但是考虑到有些图片中会有一些干扰的白色背景部分，这样就会造成大气光估计不够准确。

function new_atmosphere=get_new_atmosphere(image)
%RGB变换到HSV
rgb=im2double(image);
r=rgb(:,:,1);
g=rgb(:,:,2);
b=rgb(:,:,3);
V=max(max(r,g),b);
tmp1=min(min(r,g),b);
tmp2=r+g+b;
tmp2(tmp2==0)=eps;
S=1-3.*tmp1./tmp2;
%求亮色区域
[m1,n1]=size(S);
C=zeros(m1,n1);
for i=1:m1
    for j=1:n1
        if S(i,j)<=0.2
            C(i,j)=1;
        else
            C(i,j)=0;
        end
    end
end

l=zeros(m1,n1);
P_sky=zeros(m1,n1);
L=m1;
for i=1:m1
    for j=1:n1
        l(i,j)=i;
        P_sky(i,j)=1-l(i,j)./L;
    end
end
SKY=P_sky(1:fix(m1*n1*0.1));
N=size(SKY);
atmosphere=mean(V(1:N));
new_atmosphere=[atmosphere,atmosphere,atmosphere];

• 获得辐射值

function radiance = get_radiance(image, transmission, atmosphere)
[m, n, ~] = size(image);
rep_atmosphere = repmat(reshape(atmosphere, [1, 1, 3]), m, n);
max_transmission = repmat(max(transmission, 0.1), [1, 1, 3]);
radiance = ((image - rep_atmosphere) ./ max_transmission) + rep_atmosphere;
end

• 引导滤波

function q = guided_filter(guide, target, radius, eps)
[h, w] = size(guide);
  
avg_denom = window_sum_filter(ones(h, w), radius);

mean_g = window_sum_filter(guide, radius) ./ avg_denom;
mean_t = window_sum_filter(target, radius) ./ avg_denom;

corr_gg = window_sum_filter(guide .* guide, radius) ./ avg_denom;
corr_gt = window_sum_filter(guide .* target, radius) ./ avg_denom;

var_g = corr_gg - mean_g .* mean_g;
cov_gt = corr_gt - mean_g .* mean_t;

a = cov_gt ./ (var_g + eps);
b = mean_t - a .* mean_g;

mean_a = window_sum_filter(a, radius) ./ avg_denom;
mean_b = window_sum_filter(b, radius) ./ avg_denom;

q = mean_a .* guide + mean_b;

end

• 总的图像去雾的函数

function [ radiance ] = dehaze_fast( image, omega, win_size )
%DEHZE Summary of this function goes here
%   Detailed explanation goes here

if ~exist('omega', 'var')
    omega = 0.95;
end

if ~exist('win_size', 'var')
    win_size = 15;
end

r = 15;
res = 0.001;

[m, n, ~] = size(image);
dark_channel = get_dark_channel(image, win_size);
atmosphere = get_atmosphere(image, dark_channel);
trans_est = get_transmission_estimate(image, atmosphere, omega, win_size);
x = guided_filter(rgb2gray(image), trans_est, r, res);
transmission = reshape(x, m, n);
radiance = get_radiance(image, transmission, atmosphere);
end

• 运行demo

tic;
image = double(imread('forest.jpg'))/255;

image = imresize(image, 0.6);

result = dehaze_fast(image, 0.95, 5);
toc;

figure, imshow(image)
figure, imshow(result)