遥感图像处理学习笔记(1)——matlab中的图像去雾

自柴静《穹顶之下》以来，雾霾一直被持续关注和监测。雾霾天气在我国多个地方频繁出现，不仅影响人们日常生活，还会对某些学科领域的研究带来一定的影响，在雾霾形成的区域，空气中的各种微小粒子增多，其带来的散射效应导致了光在传播过程衰减，使得光信号承载的信息丢失。这种效应在遥感领域造成的影响颇为重大，令图像的后续处理变的困难。所以对有雾遥感图像的恢复算法有着重要的价值。

本文使用matlab进行图像除雾处理，核心算法来源于Single image haze removal using dark channel prior. 这篇文章，即暗通道先验理论。下面按理论说明与算法实现两部分展开。

1、除雾理论

暗通道理论的核心是：在有雾图像中，绝大多数非天空的局部区域里，某一些像素总会有至少一个颜色通道具有很低的值。暗通道的表达式为：

J d a r k ( x ) = min ⁡ y ∈ Ω ( x ) ( min ⁡ C ∈ { r , g , b } J C ( y ) ) ⋯ ⋯ 1 J^{dark}(x)=\min_{y \in \Omega (x)} (\min_{C \in {\left \{ r,g,b \right \} }}J^C(y))\cdots \cdots 1 Jdark(x)=y∈Ω(x)min​(C∈{r,g,b}min​JC(y))⋯⋯1

其中J表示图像，Jc为图像的一个颜色通道（例如Jr表示图像红色通道）, Ω(x)表示以像素X为中心的一个窗口。再来表示暗通道理论：
$$J^{dark}\to 0\cdots \cdots 2$$

暗通道理论不难理解，难的是如何在matlab中运用这个理论来除雾。除了暗通道理论，要想除雾还必须掌握另外一个模型——雾图形成模型：
$$I(x)=J(x)t(x)+A(1-t(x))\cdots \cdots 3$$

其中，I(x)表示接收到的图像，J(x)表示原始无雾的图像；A定义为全球大气光成分，t(x)为透射率。接下来要做的就是根据已知的I(x)一步步求出J(x)。

对3式做变换：

$$\frac{I^C(x)}{A}=\frac{J^C(x)}{A}t(x)+1-t(x)\cdots \cdots 4$$ c ∈ { r , g , b }   c\in\left\{r,g,b\right\} \\\ c∈{r,g,b} 

假设在每一个窗口内,透射率t(x)为常数，定义为t ̃(x)，并且A也为常数，再对上式两边同时取最小值：
$$\begin{gathered} \underset{y\in \Omega (x)}{\min }(\underset{C}{\min }\frac{I^C(y)}{A^C})= \\ \tilde{t}(x)\underset{y\in \Omega (x)}{\min }\left(\underset{C}{\min }\frac{J^C(y)}{A^C}\right)+1-\tilde{t}(x)\cdots \cdots 5 \end{gathered}$$

由2式：
$$J^{\text{dark}}=\underset{y\in \Omega (x)}{\min }\left(\underset{c}{\min }{J}^c(y)\right)=0\cdots \cdots 6$$
将6式带入5，最后得到透射率的预估值：
$$\tilde{t}(x)=1-\underset{y\in \Omega (x)}{\min }\left(\underset{c}{\min }\frac{I^c(y)}{A^c}\right)\cdots \cdots 7$$
查阅文献了解到，为了避免过度去雾，通常在去雾的时候需要保留一定程度的雾，通过在透射率表达式中引入一个在[0,1]之间的修正因子ω，查阅到常用的ω有0.95，0.9等。
$$\tilde{t}(x)=1-\omega \underset{y\in \Omega (x)}{\min }\left(\underset{c}{\min }\frac{I^c(y)}{A^c}\right)\cdots \cdots 8$$
最终求解J(x):
$$J(x)=\frac{I(x)-A}{t(x)}+A\cdots \cdots 9$$
又从文献了解到，透射率图中t(x)的值很小时，会导致J的值偏大，从而使得图像整体向白场过度，因此一般设置一阈值T0，当t值小于T0时，令t=T0，后文中均以T0=0.1为标准计算。
$$J(x)=\frac{I(x)-A}{\max \left(t(x),T_0\right)}+A\cdots \cdots 10$$
根据以上表达式，可以整理出从有雾图像I(x)中获得原图像J(x)的步骤如下：
1、由1获得I(x)的暗通道I^dark (x)；
2、根据A的定义，在暗通道中寻找具有最大亮度值的区域，根据该区域回到原图中，对应的最高亮度作为A值；
3、分别求出每个通道的大气光值A^c，根据8求出t(x)；
4、通过10求出J(x).

2、matlab代码实现

首先要获取图像的暗通道，其中有两个最小值：一个是像素点亮度值在三原色通道中的最小值，一个是以像素为中心的邻域中的最小值。对于前者，只用一行代码就可以求解：

minc=min(pic,[],3);

这条表达式的含义是图像矩阵pic在第三个维度上的最小元素。通过这行代码，我们把三维图像矩阵化成了二维的，其中每个像素点的值是原图中三原色通道的最小值。接下来我们需要求出，在以每个像素为中心的窗口中的最小值。仔细一想，图像卷积处理的原理是用一个n×n的矩阵，依次遍历图像每一个像素点做矩阵相乘。这给我们提供了一个很棒的思路：我们使用：
$$\left(\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\\ 1& 1& 1\end{array}\right)$$这一矩阵对上一步求得的二维矩阵进行卷积处理，并取每次运算得到的矩阵中的最小值，就可以解出暗通道了。查询matlab文档了解到可以使用ordfilt2函数，它的文档描述如下：

B = ordfilt2(A,order,domain) replaces each element in A by the orderth element in the sorted set of neighbors specified by the nonzero elements in domain.

也就是说，用domain作为卷积模板，取其中的第order个值代替中心值。我们这里A就是待处理的I(x)，order取1表示最小值，domain设置一个15×15的、值全部为1的矩阵求解：

domain = ones(15)；
I_dark=ordfilt2(minc, 1, domain, ‘symmetric’);

至此便完成了去雾的第一步——获取暗通道。接下来需要求出“两个”大气光值：全图的大气光值A，和每个通道中的大气光值A^c。这里我使用了一种迭代算法求全图的A：设想一张有雾图片中，总有一部分包括天空（或雾气浓度最大），因此首先将一张图片四等分，求出4个部分的平均亮度值，选择最大的那个部分，再次四等分，重复迭代直到区域大小低于限定值，记录最终得到的区域，最后回到原图中，求解该区域的亮度平均值作为该大气光值。联想一下计算机科学中的四叉树，好像原理也类似。

代码实现

function [imax, jmax] = QuadMat(mat, minLength, i, j)
% edited by TerryYuan, 2020
% 迭代法求最大值
% 输入i, j为数组索引
% min为最小长度
% 输出为最大子树索引
h = length(i);
w = length(j);
if h>minLength && w>minLength
    imin = min(i);
    imax = max(i);
    jmin = min(j);
    jmax = max(j);
    imid = floor(h/2) + imin - 1;
    jmid = floor(w/2) + jmin - 1;
    child0 = mat(imin : imid, jmin : jmid);
    child1 = mat(imin : imid, jmid + 1 : jmax);
    child2 = mat(imid + 1 : imax, jmin : jmid);
    child3 = mat(imid + 1 : imax, jmid + 1 : jmax);
    w0 = mean(child0(:)) - power(std2(child0),2);
    w1 = mean(child1(:)) - power(std2(child1),2);
    w2 = mean(child2(:)) - power(std2(child2),2);
    w3 = mean(child3(:)) - power(std2(child3),2);
    weight=[w0,w1,w2,w3];
    maxchild = find(weight==max(weight));
    if maxchild < 3
        newi = [imin : imid];
    else
        newi = [imid + 1 : imax];
    end
    if (maxchild == 3) | (maxchild ==1)
        newj = [jmin : jmid];
    else
        newj = [jmid+1 : jmax];
    end
%     disp("子树："+ maxchild(1));
    [imax, jmax] = QuadMat(mat, minLength, newi, newj);

else
    imax = i;
    jmax = j;
end

end

使用上面所示的QuadMat函数可求解A。

接下来要根据式8求解t，仔细观察发现式8和暗通道表达式（式1）在形式上是一样的，所以我们可以对$$I^c(y)/A^c$$使用ordfilt2处理。代码实现：

function [Tx] = Trans(pic, filtersize)
% edited by TerryYuan, 2020
% 计算透射率
[H,W,~] = size(pic);
omega = 0.9;
r = pic(:,:,1);
g = pic(:,:,2);
b = pic(:,:,3);
% 四叉树迭代找大气光图像区域
% 求解各通道的Ac
[arx,ary] = QuadMat(r, 10, [1:H], [1: W]);
Ar = mean(mean(r(arx, ary)));
[agx,agy] = QuadMat(g, 10, [1:H], [1: W]);
Ag = mean(mean(g(agx, agy)));
[abx,aby] = QuadMat(b, 10, [1:H], [1: W]);
Ab = mean(mean(b(abx, aby)));
IA = cat(3, r/Ar, g/Ag, b/Ab);
Tx = 1 - omega * autoFilter(IA, filtersize);
end

到此为止我们成功求出了A和t，接下来便可根据式10计算导出图像矩阵：

J = (pic - A)./max(Tx, t0) + A;

对文章开头左边的原图按上述步骤处理，得到以下图像

发现图中轮廓四周存在白色包围，推测应该是卷积处理时的精度问题，所以将卷积模板换成3×3的矩阵，再次对图像进行去雾变换，得到以下图像

不难发现，3×3处理的图像更为精细，然而整体图像存在较为严重的偏色，可以判断为除雾过度导致的结果。

至此可以总结比较两种不同大小的卷积模板的优劣：使用15×15的模板，可以得到颜色更真实自然的图像，但图中轮廓边缘有不自然的白边；使用3×3的模板可以得到精细的去雾图，但整体色调偏差严重。显然两种卷积模板都不能满足预想的去雾要求，我尝试了各种卷积窗口，变换了修正因子ω，结果都不尽人意。再次打开知网查找解决方案，了解到：引导滤波器可以结合这两种暗通道（精细与粗糙）的优点来优化透射率，实现途径为matlab的imguidedfilter函数。文档说明为：

imguidedfilter(A,G) filters binary, grayscale, or RGB image A using the guided filter, where the filtering process is guided by image G.

其原理是利用较为精细的3×3滤波后得到的Tx3为导向图，对优化后的Tx15进行导向滤波，最终得到较为精细的Tx。代码实现步骤为：
Tx1 = Trans(pic, 15);
Tx2 = Trans(pic, 3);
Tx3 = autoFilter(Tx1, 15);
Tx4 = autoFilter(Tx3, Tx2);

说明：Trans函数可以获取图像的预估透射率，autoFilter是我编写的函数，通过第二个参数类型来选择使用ordfilt2或imguidedfilter：若第二个参数为整数，那么autoFilter对第一个参数矩阵进行卷积处理（即最小值滤波），若第二个参数为矩阵，那么使用引导滤波器处理第一个参数矩阵。代码如下

function pic_filtered = autoFilter(pic, filter)
% edited by TerryYuan, 2020
% 图像卷积处理
% 根据输入filter的类型选用最小值滤波器或引导滤波器
% 首先获取c∈(R,G,B)中最小值
minc=min(pic,[],3);
if numel(filter)==1
	domain = ones(filter);
	pic_filtered=ordfilt2(minc,1,domain, 'symmetric');
else
	pic_filtered=imguidedfilter(minc,filter);
end

end

按此方法处理得到的四个透射率，将其成图比较

利用导向滤波得到的透射率，带入式10，求出优化后的去雾图像

主程序代码

% edited by TerryYuan, 2020

close all;clear all;clc;
filename = "2.jpg";
pic=imread(filename);
[H,W,c]=size(pic);

%   暗通道原理去雾
%   uint8转double
p = im2double(pic);
%   autoFilter获取暗通道
iDark1=autoFilter(p,15) ;
iDark2=autoFilter(p, 3) ;
figure;
imshowpair(iDark1,iDark2, 'montage');
title("dark tunnel");

%   根据暗通道计算透射率
Tx1 = Trans(p, 15);
Tx2 = Trans(p, 3);
Tx3 = autoFilter(Tx1, 15);
Tx4 = autoFilter(Tx3, Tx2);

figure;
subplot(2,2,1); imshow(Tx1); title("15×15 预估透射率");
set(gca,'position',[0 0.5 0.49 0.45]);
subplot(2,2,2); imshow(Tx2); title("3×3 导向图 ");
set(gca,'position',[0.51 0.5 0.49 0.45]);
subplot(2,2,3); imshow(Tx3); title("15×15 优化透射率");
set(gca,'position',[0 0.01 0.49 0.45]);
subplot(2,2,4); imshow(Tx4); title("使用导向滤波后的透射率");
set(gca,'position',[0.51 0.01 0.49 0.45]);

[imax, jmax] = QuadMat(iDark1, 10, [1:H], [1:W]);
A = mean(mean(p(imax, jmax)));
A = min(A, 0.7);
%   求大气光值A
%   多次试验后，为了方式天空区域过曝要限制A不能过大

% J = J + (255-J).*J.*0.1;
% 亮度补偿

t0 = 0.1;
J1 = (p - A)./max(Tx1, t0) + A;
J1 =  im2uint8(J1);
figure;
imshowpair(p, J1, 'montage');
title("Filter 15×15");

J2 = (p - A)./max(Tx2, t0) + A;
J2 =  im2uint8(J2);
figure;
imshowpair(p, J2, 'montage');
title("Filter 3×3");

J = (p - A)./max(Tx4, t0) + A;
J = im2uint8(J);
figure('NumberTitle', 'off', 'name','处理后');
imshowpair(p, J, 'montage');
title({["优化导向滤波"],["ω=0.7"]});

figure;
subplot(2,2,1); imshow(p); title("Origin");
set(gca,'position',[0 0.5 0.49 0.45]);
subplot(2,2,2); imshow(J1); title("Filter 15×15");
set(gca,'position',[0.51 0.5 0.49 0.45]);
subplot(2,2,3); imshow(J2); title("Filter 3×3");
set(gca,'position',[0 0.01 0.49 0.45]);
subplot(2,2,4); imshow(J); title("Filter Optimized");
set(gca,'position',[0.51 0.01 0.49 0.45]);

% 保存图像
% name = split(filename,'.');
% name = strcat(name(1), '-filter.jpg');
% imwrite(J,name);