计算视觉——基于暗通道先验的图像去雾算法雾在摄影时，是一种显著降低图像质量、影响观者视觉体验的干扰因素，虽然在特定情况下

问题背景

雾在摄影时，是一种显著降低图像质量、影响观者视觉体验的干扰因素，虽然在特定情况下，雾是艺术摄影中的重点关注对象，但对于用户日常的拍摄以及地理成像等多种图像应用领域而言，我们更希望能够减少雾的影响。因雾而降低质量的图像如下

为了解决这一类因雾霾气候所造成的图像质量降低，图像去雾的研究领域展开了一系列相关的研究。其中不乏有许多精彩的想法与令人惊叹的实现，本文的重点将会聚焦于计算机视觉顶级会议CVPR2009年的会议论文基于暗通道先验的单图像去雾(Single Image Haze Removal Using Dark Channel Prior, 2009 CVPR Best Paper Award, 2009CVPR最佳论文奖)，该论文所提出的方法和思想非常经典且简单，是图像领域中的学习者、研究者都应该了解的。

此外，本文将从易到难的顺序介绍包含暗通道去雾在内的三种去雾方法，帮助读者深入理解图像去雾中的难点与各方法的贡献与突破。

暴力解法——直方图均衡化去雾

其实不难发现，有雾的图像往往会偏暗或偏白，通过图像三通道直方图我们可以很容易发现这一特征。

通过直方图规律上的统计，我们很容易发现，有雾图像三通道值的直方图分布是不均匀的，即分布向一边倾移。基于这一发现，我们可以把去雾问题转化为直方图的分布不均匀问题去解决，这就可以引出我们的第一种"暴力"解法——直方图均衡化。直方图均衡化这是图像领域非常常见的处理手段，因此本文不会过多解释。

% Matlab 2017a
% Histogram_ equalization_enhancement.m

RGB = imread('C:\Users\as\Desktop\juejin\3——\1.png'); % 读取待去雾彩色图像
subplot(121);
imshow(RGB);
title('待去雾图像');

[R, C, K] = size(RGB); % 新增的K表示颜色通道数

% 统计每个像素值出现次数
cnt = zeros(K, 256);
for i = 1 : R
    for j = 1 : C
        for k = 1 : K
            cnt(k, RGB(i, j, k) + 1) = cnt(k, RGB(i, j, k) + 1) + 1;
        end
    end
end

f = zeros(3, 256);
f = double(f); cnt = double(cnt);

% 统计每个像素值出现的概率， 得到概率直方图
for k = 1 : K
    for i = 1 : 256
        f(k, i) = cnt(k, i) / (R * C);
    end
end

% 求累计概率，得到累计直方图
for k = 1 : K
    for i = 2 : 256
        f(k, i) = f(k, i - 1) + f(k, i);
    end
end

% 用f数组实现像素值[0, 255]的映射。 
for k = 1 : K
    for i = 1 : 256
        f(k, i) = f(k, i) * 255;
    end
end

% 完成每个像素点的映射
RGB = double(RGB);
for i = 1 : R
    for j = 1 : C
        for k = 1 : K
            RGB(i, j, k) = f(k, RGB(i, j, k) + 1);
        end
    end
end

% 输出
RGB = uint8(RGB);
subplot(122);
imshow(RGB);
title('直方图均衡化增强去雾图像');

处理结果如下：

我们可以进行简单的结果分析，在增强图像中其实存在问题，如图像6中的建筑虽然变得清晰了，但叶子却因为直方图的均衡化而发黑，树叶部分的细节信息失真严重。

而图像1与图像5都存在这类问题，这是因为景深变化而导致的图像各位置的雾的浓度不同所引起的，因此直方图均衡化在对于整幅图像含雾程度差不多的情况下才能有理想的效果。

基于Retinex理论的图像去雾

Retinex理论是一种表征人对色彩感知的一个模型。Retinex理论的基本内容是物体的颜色是由物体对长波（红），中波（绿）和短波（蓝）光线的反射能力决定的，而不是由反射光强度的绝对值决定，物体的色彩不受光照非均匀性的影响，具有一致性，即Retinex理论是以颜色恒常性为基础的。如下图所示，观察者所看到的物体的图像颜色是由物体表面反射率对入射光反射得到的，反射率R由物体本身决定，不受入射光大小变化。

Retinex这个词是由视网膜（视网膜Retina）和大脑皮层（皮层Cortex）两个词组合构成的，是为了表明当时尚不明确HVS(人类视觉系统,Human Vision System)特性究竟取决于此两个生理结构中的哪一个有关。Retinex理论是基于人类视觉的亮度和颜色感知的模型所提出了一种颜色恒常知觉的计算理论，即通常情况下，人眼/大脑消除光源对物体表面色彩的影响。基于Retinex理论，我们可以把图像色彩分解为光源与物体反射两部分，且两部分对最终色彩的影响是乘性的。

需要指出，现在对于人眼色彩感知的理论是基于目前的阶段理论，视网膜中的视锥细胞感知色彩、视杆细胞感知明暗：故在黑暗处我们无法有较好的色彩感知，是因为视锥细胞在暗处几乎不发挥作用，而视杆细胞占据主导地位；在亮度足够时，视锥细胞会占视觉的主导作用；最终由二者采集的光信息转化为电/化学信号输送至大脑，形成人的视觉感知。阶段理论既可以解释颜色颉颃，又能够说明三原色光可以混合出所有色彩的可见光，是目前公认的颜色感知理论。）

因此，根据Retinex理论，一副图像 $S(x,y)$ 可以分解为反射图像 $R(x,y)$ 和亮度图像 $L(x,y)$ ，对于图像 $S$ 中的每个像素点 $(x,y)$ ，我们可以用如下公式表示三者间的关系：

实际上，Retinex理论是利用图像 $S$ 来计算得到物体反射率 $R$ ，也就是消除入射光 $L$ 的影响，从而使得图像呈现未受光源影响的"原本"的颜色（这里的"原本"必须打引号，因为物体本身只有光谱反射率，没有颜色，只有当光线照射在物体上，反射光进入人眼，才产生了颜色感知，即人类视觉系统、光、物体缺一就不会形成颜色；需要强调，颜色不是物理性质）。

单尺度的Retinex算法（SSR Algorithm, Single Scale Retinex Algorithm）是最简单最早提出的一种基于Retinex算法，其算法步骤如下：

Step 1:通过对数域将乘性分量 $L$ 和 $R$ 转换为加性分量

Step 2:用高斯模板对原图像做卷积，得到亮度估计图像 $D(x,y)$ （该步骤的物理解释就是通过对原图像 $S$ 中像素点与周围区域中像素的进行加权平均来对图像中照度变化做估计，即用区域像素估计该区域的亮度信息，在代码中，我们采用的是对396x380图像采用了396x380同样大小的卷积窗口，因此是对整幅图像估计一个光源亮度值）

其中 $F(x,y)$ 为高斯函数，表示为

F(x,y)=\lambda e^{\frac{-(x^2+y^2)}{c^2}}

我们在matlab中通过如下代码创建该模板，其中N1、M1为原图像长、宽上的像素数

sigma = 255;
F = fspecial('gaussian', [N1,M1], sigma);

Step 3:在对数域中消去估计光源的影响，得增强图像的对数域形式 $\log{R(x,y)}$ ，并将其转换到实数域 $R(x,y)$ ，如下

由此，SSR算法的步骤已全部介绍完毕，下面贴上代码。

% Matlab 2017a
% Single_Scale_Retinex_Algorithm.m

clear
I = imread('C:\Users\as\Desktop\juejin\201552036503611.png');

R = I(:, :, 1);
[N1, M1] = size(R);
R0 = double(R);
Rlog = log(R0+1);
Rfft2 = fft2(R0);
 
sigma = 255;
F = fspecial('gaussian', [N1,M1], sigma);
Efft = fft2(double(F));

DR0 = Rfft2.* Efft;
DR = ifft2(DR0);

DRlog = log(DR +1);
Rr = Rlog - DRlog;
EXPRr = exp(Rr);
MIN = min(min(EXPRr));
MAX = max(max(EXPRr));
EXPRr = (EXPRr - MIN)/(MAX - MIN);
EXPRr = adapthisteq(EXPRr);

G = I(:, :, 2);

G0 = double(G);
Glog = log(G0+1);
Gfft2 = fft2(G0);
  
DG0 = Gfft2.* Efft;
DG = ifft2(DG0);

DGlog = log(DG +1);
Gg = Glog - DGlog;
EXPGg = exp(Gg);
MIN = min(min(EXPGg));
MAX = max(max(EXPGg));
EXPGg = (EXPGg - MIN)/(MAX - MIN);
EXPGg = adapthisteq(EXPGg);

B = I(:, :, 3);

B0 = double(B);
Blog = log(B0+1);
Bfft2 = fft2(B0);
  
DB0 = Bfft2.* Efft;
DB = ifft2(DB0);

DBlog = log(DB+1);
Bb = Blog - DBlog;
EXPBb = exp(Bb);
MIN = min(min(EXPBb));
MAX = max(max(EXPBb));
EXPBb = (EXPBb - MIN)/(MAX - MIN);
EXPBb = adapthisteq(EXPBb);

result = cat(3, EXPRr, EXPGg, EXPBb);
subplot(121), imshow(I);
subplot(122), imshow(result);

上图是通过上述代码（参数默认未修改）所运行得到的结果，值得一提的是，SSR算法不仅适用于图像去雾增强，在白平衡中使用同样是可行的，见上图书本色彩。

基于暗通道先验的单图像去雾

前面介绍的两种方法其实都存在一个问题，就是因距离变化而导致图像局部区域有雾程度不同所引起的去雾残留（即图像中的雾并没有去除完全），而何凯明（暗通道去雾的提出者）更是给出了一个绝妙的想法解决了这一领域长期存在的问题。不多吹，看完自会明白提出该方法的论文为何能成为计算机视觉顶会的最佳论文，下面我们进入正题。

算法基本思想是：在绝大多数非天空的局部区域里，某些像素至少会有一个很低的颜色通道的值，这一点是作为先验知识进行图像去雾的有力佐证。对于任意一副图像 $J$ 。该图像的暗通道(Dark Channel Prior)可以用如下用数学语言表示

J^{Dark}(x)=\mathop{\min(\min(J^c(y)))}\limits_{c\in \{r, g, b \} y\inΩ(x)}

其中， $c$ 代表通道，可以为 ${r,g,b}$ 任意一者， $Ω(x)$ 为以像素x为中心的一个窗口，窗口大小未预先定义； $y$ 是 $Ω(x)$ 区域中的其中一个点。

而基于作者大量观察没有雾图像的暗通道作为先验，结论指出

到这里为止仍然有些晦涩难懂，用通俗易懂的语言来说，作者发现，在非天空区域的无雾图像中的局部区域(比如20x20、30x30、50x50)，在该区域内R、G、B三通道中存在的最小值趋于0（即在区域窗口内搜索r、g、b的最小值，并非三个都需要求最小，而是区域内三者中的最小值，一个区域一个最小值，若该图像没有雾，那么最小值应当趋于0。）

给出该部分的代码实现。

% Matlab 2017a
% DarkPrior.m
clc
clear
close all
imageRGB = imread('C:\Users\as\Desktop\juejin\testimage1.png');
figure;
imshow(imageRGB);
imageRGB = double(imageRGB);
imageRGB = imageRGB./255;

r = imageRGB(:,:,1);
g = imageRGB(:,:,2);
b = imageRGB(:,:,3);
[m n] = size(r);
a = zeros(m,n);
for i = 1:m;
    for j = 1:n;
        a(i,j) = min(r(i,j),g(i,j));

        a(i,j) = min(r(i,j),b(i,j));
    end
end
d = ones(15,15); % 15x15 window
fun = @(block_struct)min(min(block_struct.data))*d;
dark = blockproc(a, [15 15], fun);   % 15x15 window
dark = dark(1:m, 1:n);

figure;
imshow(dark);

有雾图像的暗通道如下：

直方图去雾图像以及自然无雾图像的暗通道图如下：

其中的显著特点是，有雾图像或去雾不完全的图像，在暗通道的图像仍然偏白，而自然图像与已去雾区域的暗通道是非常黑、非常暗的。这也是该通道被称为暗通道的原因。

由这些图像我们可以很明显看到暗通道先验理论的普遍性。在作者何凯明的论文中，统计了5000多副图像的特征，也都基本符合这个先验，因此，我们可以认为这是一条定理。作者采用的窗口大小相对比博主使用的更大一些，结果如下。

而暗通道在求解雾图的透射率中非常有用，下面会讲到雾图像的物理模型。

在计算机视觉和计算机图形学的研究中（是的你没有看错，计算机图形学，原因在于电子游戏中有一些非常逼真的浓雾场景，也同样是采用下面这个公式去生成的），下述方程所描述的雾形成模型被广泛使用：

其中， $I^c(x)$ 代表待去雾图像， $J^c(x)$ 代表无雾图像， $A^c$ 表示大气光成分， $t(x)$ 为透射率， $c \in \{r, g, b\}$ 。

这个模型的可解释性很强，如果大气中没有雾，那么透射率就是 $t=100\%$ ，则 $I^c(x)=J^c(x)$ 成立，图像本身就没有雾。反之，如果大气中全是雾(实际上不可能)，那么透射率 $t=0\%$ ，有雾图像 $I^c(x)$ 和无雾图像 $J^c(x)$ 一点关系都没有，此时 $I^c(x)=A^c$ ，图像全部颜色直接与大气中的雾颜色完全相同。

通过上面这些解释，大家可以理解这个雾的模型了。那么我们现在需要做的，其实就是已知 $I(x)$ ，求解 $J(x)$ 。

而如果要单用这一个模型，通过已有的待去雾图像 $I(x)$ 计算得到去雾图像 $J(x)$ ，这本质上是个病态问题，已知参数显然不足，我们既不知道透射率 $t(x)$ ，也不知道大气光成分 $A^c$ 。

但是，大气光成分可以大致进行推测，何凯明给出的一种方法是，取暗通道中亮度前 $1\%$ 的像素位置对应的像素的三通道分别求平均，这样所得的三个值可以大致作为大气光成分的三通道值（图像中纯粹是雾时的三通道值）。

那么，现在这个模型中只剩下透射率 $t(x)$ 无法得知了。接下来就是劈头盖脸一通朴实无华的公式推导，建立起暗通道与透射率的关系。

首先，将雾模型等式两边同时除以 $A^c$ ，如下

\frac{I^c(x)}{A^c}=\widetilde t(x)\frac{J^c(x)}{A^c}+1-\widetilde t(x)

我们不妨先假设每一个窗口内的透射率 $t(x)$ 为常数 $\widetilde t(x)$ 。对于上式，同时对两侧进行最小化运算，如下：

\mathop{\min}\limits_{ y\inΩ(x)} (\mathop{\min} \limits_c \frac{I^c(x)}{A^c})=\widetilde t(x) \mathop{\min}\limits_{ y\inΩ(x)} (\mathop{\min} \limits_c\frac{J^c(x)}{A^c})+1-\widetilde t(x)

还记得之前暗通道的定义吗，暗通道有如下定义，且在无雾图像中趋向于0，如下：

J^{Dark}(x)=\mathop{\min(\min(J^c(y)))}\limits_{c\in \{r, g, b \} y\inΩ(x)} \rightarrow 0

将上式定义代入前一个个式子，则有

\widetilde t(x)=1-\mathop{\min}\limits_{y\inΩ(x)}(\mathop{\min} \limits_c\frac{I^c(x)}{A^c})

而有雾的图像的暗通道是我们可以求得的，那么透射率 $\widetilde t(x)$ 我们也已经求得了。

不过，为了保留部分图像景深，作者采用引入了一个常数保留极少数的雾，如下：

\widetilde t(x)=1-0.95 \times \mathop{\min}\limits_{y\inΩ(x)}(\mathop{\min} \limits_c\frac{I^c(x)}{A^c})

通过雾模型，我们可以得到下式：

J^c(x)=\frac{I^c(x)-A^c}{\max(t(x), t_0)+A}

其中 $\max(t(x), t_0)$ 是为了避免分母中出现 $t(x)=0$ 引起的问题， $t_0=0.1$ 。即求解透射率小于0.1时，不使透射率小于0.1（局部区域雾影响太严重时，图像已严重失真，去雾只会放大噪声，进一步降低图像质量）。

由此，已对基于暗通道去雾的思想和方法进行了详细的解读，下面将进入代码环节。

代码实战部分

% Matlab 2017a 
% dehaze_without_guide_filter.m
clear
I=double(imread('C:\Users\as\Desktop\juejin\testimage123.png')) / 255;
w0 = 0.95 ;%去雾系数选0.95
[h,w,c]=size(I) ;
dark_channel=zeros(h,w);
for i= 1 : h
    for j= 1 : w
        dark_channel(i,j)=min(I(i,j,:)) ;
    end
end
%求到暗通道，并对它进行最小值滤波
dark_channel=ordfilt2(dark_channel,1,ones(9,9),'symmetric');
%求大气亮度值A
A = max(max(dark_channel));
[a,b]=find(dark_channel == A);
a=a(1);
b=b(1);
A=mean(I(a,b,:));
%求初始透射率
transmission=1 - w0 * dark_channel ./ A;
t0=0.1;
t=max(transmission,t0);
dehaze=zeros(h,w,c);
for l = 1:c
    dehaze(:,:,l)=(I(:,:,l)-A)./t + A;
end
figure
imshow(dehaze);

运行结果分析

实战部分的代码经过非常高的简化，基本是目前能给出的matlab中的最优代码，虽然看起来很短，但是很多部分如果靠自己编程去实现还是有很多难点。下面给出结果对比。

孰优孰劣一看便知，但暗通道仍然存在一定问题，由于我们是通过区域计算暗通道，所以投射率图、暗通道图都存在一定程度的像素化（马赛克化），这导致了边缘位置的不连续性，如(d)的边缘位置残留的雾区域。但这可以通过导向滤波进行优化，代码如下：

% Matlab 2017a 
% dehaze_with_guide_filter.m
clear
I=double(imread('C:\Users\as\Desktop\juejin\testimage123.png')) / 255;
w0 = 0.95 ;%去雾系数选0.95
[h,w,c]=size(I) ;
dark_channel=zeros(h,w);
for i= 1 : h
    for j= 1 : w
        dark_channel(i,j)=min(I(i,j,:)) ;
    end
end
%求到暗通道，并对它进行最小值滤波
dark_channel=ordfilt2(dark_channel,1,ones(9,9),'symmetric');
%求大气亮度值A
A = max(max(dark_channel));
[a,b]=find(dark_channel == A);
a=a(1);
b=b(1);
A=mean(I(a,b,:));
%求初始透射率
transmission=1 - w0 * dark_channel ./ A;
guided_image=I(:,:,1);
%直接调用引导滤波，对初始透射率进行引导滤波
transmission2=imguidedfilter(transmission,guided_image,'NeighborhoodSize',[30,30]);
t0=0.1;
% t=max(transmission,t0);
t=max(transmission2,t0);
dehaze=zeros(h,w,c);
for l = 1:c
    dehaze(:,:,l)=(I(:,:,l)-A)./t + A;
end
figure
imshow(dehaze);

优化前后结果对比：

其中，softmax是作者给出的一种软抠图方法，效果上来说确实比导向滤波更好，但计算效率远远不如导向滤波，有兴趣的读者朋友可以自行了解softmax。

下面再给出两组测试图像的结果。

在此之中，有效果相当卓越的存在，也有一些不尽人意的结果。虽然暗通道具有一定的参数依赖（受窗口大小等参数影响）以及在天空、柏油马路等区域(这两种情况下虽然没有雾，但局部区域RGB最小值仍然会很大，导致计算所得的暗通道值较大，被错误地判别为有雾区域进行去雾，这会导致非常影响视觉体验的图像噪声出现)有非常差的表现，但这一创造性的思想是非常值得了解学习的。通过几个公式推导和两个假设就能实现的暗通道先验去雾，是去雾领域独树一帜的方法。

在此之后有非常多基于暗通道先验进行创新的去雾算法，有兴趣的朋友可以自己再深入研究。

写在最后

Paper on "Single Image Haze Removal Using Dark Channel Prior" :

www.csie.ntu.edu.tw/~b98902025/…

如果你喜欢我的这篇博客的话，请点个赞吧！~如果你喜欢我的风格并且想了解更多相关的内容，可以关注我~我会不定期更新有趣的图像与视觉方面的知识并加入个人的独特理解~

您的支持就是对我最大的肯定。