《Single Image Haze Removal Using Dark Channel Prior》去雾代码实现分析

 本文主要为了通过对代码进行分析，帮助更进一步了解《Single Image Haze Removal Using Dark Channel Prior》的操作步骤。在看本篇文章前，需要对《Single Image Haze Removal Using Dark Channel Prior》有大致的了解。提供一个文章解析参考的连接：https://www.cnblogs.com/Imageshop/p/3281703.html。本文的代码分析也是基于上述文章中的代码资源进行分析的，本文主要在于结合代码梳理出更加清晰的思路，并且对代码进行了较为详细的注释。

一、主程序
首先在这给出主程序及其相关注释

clear
clc
close all

kenlRatio = 0.01;
minAtomsLight = 240;            %设定的大气光的上限值
% image_name =  'test images21.bmp';
image_name =  'D:1.png';       %原图像
img=imread(image_name);         %读取原始图片文件中的数据
figure,imshow(uint8(img)), title('src');   %显示原始图像  并且显示标题 

sz=size(img);       %返回一个行向量sz，该行向量的第一个元素时矩阵的行数，第二个元素是矩阵的列数。
w=sz(2);            %读取所提取矩阵的列数
h=sz(1);             %读取所提取矩阵的行数
dc = zeros(h,w);      %返回一个h×w的0矩阵给dc，实际上这个矩阵的size=原输入图像size

for y=1:h               %行遍历          
    for x=1:w           %列遍历
        dc(y,x) = min(img(y,x,:));  %对原图像的每一个像素求RGB3个通道中的最小亮度值,然后赋值给dc矩阵的y行x列
    end
end

figure,imshow(uint8(dc)), title('Min(R,G,B)');  %显示原图像的灰度图

krnlsz = floor(max([3, w*kenlRatio, h*kenlRatio])) %krnlsz为最小滤波窗口的边长尺寸大小
%kenlRatio是一个预设的值，求尺寸时拿kenlRatio和灰度图矩阵长和宽（即w和h）乘上这个预设值和3进行比较，选择最大数作为窗口尺寸
dc2 = minfilt2(dc, [krnlsz,krnlsz]);        %对得到的灰度图进行最小值滤波获得暗通道图像

dc2(h,w)=0;             %dc2矩阵的h行w列赋值0
figure,imshow(uint8(dc2)), title('After filter ');  %显示滤波后的图像
t = 255 - dc2;          %下面4行对应公式11，即求出传输特性t（x)的估计值
figure,imshow(uint8(t)),title('t');     %显示传输特性t的图像(原始的方法求t）
t_d=double(t)/255;
sum(sum(t_d))/(h*w)     %t(x)表示的是x像素出的t的估计值，而此处将t_d的矩阵进行平均处理求出一个平均的t值


A = min([minAtomsLight, max(max(dc2))]) %根据论文提示的方法求A值，不过为了防止A过高，设置了一个上限值minAtomsLight
                                         %max(max(dc2))表示从暗通道图中按照亮度的大小取前0.1%的像素。在这些位置中，在原始有雾图像I中寻找对应的具有最高亮度的点的值，作为A值
J = zeros(h,w,3);       %生成一个0矩阵J 大小为h×w 此处的3表示rgb3个通道
img_d = double(img);    %原始图像类型转换成double

J(:,:,1) = (img_d(:,:,1) - (1-t_d)*A)./t_d; %根据公式22求出图像矩阵所有元素的R通道的值，下面两行为G和B的值
J(:,:,2) = (img_d(:,:,2) - (1-t_d)*A)./t_d;
J(:,:,3) = (img_d(:,:,3) - (1-t_d)*A)./t_d;

figure,imshow(uint8(J)), title('J');    %显示去雾图像J


%----------------------------------
%由于透射率图过于粗糙，而通过导向滤波的方式来获得较好的透射率图，下面是导向滤波计算t的方式
r = krnlsz*4
eps = 10^-6;

% filtered = guidedfilter_color(double(img)/255, t_d, r, eps);
filtered = guidedfilter(double(rgb2gray(img))/255, t_d, r, eps);%导向滤波计算t

t_d = filtered;%导向滤波计算的值赋值给t_d矩阵
figure,imshow(t_d,[]),title('filtered t');  %显示导向滤波后求出的t值

J(:,:,1) = (img_d(:,:,1) - (1-t_d)*A)./t_d; %根据公式22求出图像矩阵所有元素的R通道的值，下面两行为G和B的值
J(:,:,2) = (img_d(:,:,2) - (1-t_d)*A)./t_d;
J(:,:,3) = (img_d(:,:,3) - (1-t_d)*A)./t_d;

imwrite(uint8(J),'D:11.bmp');%将去雾后的图像数据写入到文件中
figure,imshow(uint8(J)), title('J_guild_filter');%显示导向滤波后求出的去雾图像J

 

二、主程序拆解分析
1.首先，根据文章我们可以得知：基于暗通道先验的思想，我们可以先求出t的预估值：
在暗通道先验的原理下，求t预估值的公式转化为：

 

对于一张确定的图，我们认为Ac大气光是一个定值。
我们首先计算出原输入图像的灰度图，即括号内的部分minc Ic(y)。大概思路是对原图像的所有像素x求RGB通道中亮度最小的值，从而得到灰度图。

image_name =  'D:1.png';       %原图像
img=imread(image_name);         %读取原始图片文件中的数据
figure,imshow(uint8(img)), title('src');   %显示原始图像  并且显示标题 
sz=size(img);       %返回一个行向量sz，该行向量的第一个元素时矩阵的行数，第二个元素是矩阵的列数。
w=sz(2);            %读取所提取矩阵的列数
h=sz(1);             %读取所提取矩阵的行数
dc = zeros(h,w);      %返回一个h×w的0矩阵给dc，实际上这个矩阵的size=原输入图像size
for y=1:h               %行遍历          
    for x=1:w           %列遍历
        dc(y,x) = min(img(y,x,:));  %对原图像的每一个像素求RGB3个通道中的最小亮度值,然后赋值给dc矩阵的y行x列
    end
end
figure,imshow(uint8(dc)), title('Min(R,G,B)');  %显示原图像的灰度图

对于公式11减号后面内容的实现：对上面得到的灰度图像进行最小值滤波，minfilt2是一个最小值滤波的函数。

kenlRatio = 0.01;
krnlsz = floor(max([3, w*kenlRatio, h*kenlRatio])) %krnlsz为最小滤波窗口的边长尺寸大小
%kenlRatio是一个预设的值，求尺寸时拿kenlRatio和灰度图矩阵长和宽（即w和h）乘上这个预设值和3进行比较，选择最大数作为窗口尺寸
dc2 = minfilt2(dc, [krnlsz,krnlsz]);        %对得到的灰度图进行最小值滤波获得暗通道图像

最后进行归一化、均值处理实现公式11的效果：

dc2(h,w)=0;             %dc2矩阵的h行w列赋值0
figure,imshow(uint8(dc2)), title('After filter ');  %滤波后的图像结果
t = 255 - dc2;          %下面4行对应公式11，即求出传输特性t（x)的估计值
figure,imshow(uint8(t)),title('t');     %显示传输特性t的图像(原始的方法求t）
t_d=double(t)/255;
sum(sum(t_d))/(h*w)     %t(x)表示的是x像素出的t的估计值，而此处将t_d的矩阵进行平均处理求出一个平均的t值

2.预估了传输特性t值后，现在根据文章求出大气光A的估值。

minAtomsLight = 240;            %设定的大气光的上限值
A = min([minAtomsLight, max(max(dc2))]) %根据论文提示的方法求A值，不过为了防止A过高，设置了一个上限值minAtomsLight

3.根据最终的求去雾图像J的公式来还原去雾图像，根据的公式如下：

代码中的平均估值t_d替换的是公式中的max(t(x),t0)部分。

J = zeros(h,w,3);       %生成一个0矩阵J 大小为h×w 此处的3表示rgb3个通道
img_d = double(img);    %原始图像类型转换成double
J(:,:,1) = (img_d(:,:,1) - (1-t_d)*A)./t_d; %根据公式22求出图像矩阵所有元素的R通道的值，下面两行为G和B的值
J(:,:,2) = (img_d(:,:,2) - (1-t_d)*A)./t_d;
J(:,:,3) = (img_d(:,:,3) - (1-t_d)*A)./t_d;
figure,imshow(uint8(J)), title('J');    %显示去雾图像J

4.由于透射率图过于粗糙，而通过导向滤波的方式来获得较好的透射率图，下面是导向滤波计算t的方式。

r = krnlsz*4
eps = 10^-6;
filtered = guidedfilter(double(rgb2gray(img))/255, t_d, r, eps);%导向滤波计算t
t_d = filtered;%导向滤波计算的值赋值给t_d矩阵
figure,imshow(t_d,[]),title('filtered t');  %显示导向滤波后求出的t值

5.利用导向滤波计算的t_d求出最后的去雾图像

J(:,:,1) = (img_d(:,:,1) - (1-t_d)*A)./t_d; %根据公式22求出图像矩阵所有元素的R通道的值，下面两行为G和B的值
J(:,:,2) = (img_d(:,:,2) - (1-t_d)*A)./t_d;
J(:,:,3) = (img_d(:,:,3) - (1-t_d)*A)./t_d;

imwrite(uint8(J),'D:11.bmp');%将去雾后的图像数据写入到文件中
figure,imshow(uint8(J)), title('J_guild_filter');%显示导向滤波后求出的去雾图像J

三、程序运行结果

1.显示原始输入图像

2.显示原始图像的灰度图

3.显示对原始图像灰度图进行最小滤波后的结果

4.求出利用暗通道先验计算的t

5.根据求出的t和A计算出去雾图像J

6.为了细化t，此处求导向滤波后的t

7.利用导向滤波后求出的去雾图像

四、说明
这里值分析了主程序的思路，对于最小滤波和导向滤波的程序未分析，具体的代码见文章开头给出的链接。