[论文复现]何恺明博士CVPR2009去雾算法(2)

一、前言

终于，简单实现了何博士论文去雾算法的基础部分。由于CVPR2009论文中的优化方法比较麻烦，速度比较慢，何博士于2010ECCV补充了一篇Guided Image Filtering优化求解速度，此处后面的优化我直接调用了何博士的函数。

主要参考文献：

[1] He K, Jian S, Tang X. Single image haze removal using dark channel prior[C]// IEEE Conference on Computer Vision & Pattern Recognition. 2009.

[2] He K, Jian S, Tang X. Guided image filtering[C]// European Conference on Computer Vision. 2010.

[3] 论文原文及相关资料下载 http://kaiminghe.com/

关键词：暗通道先验、去雾、导向滤波

二、实现

论文求解核心公式是：

现有图已知，未知量为待求图，全局大气光和未被散射光。

A全局大气光采用亮度值(hsv-v)最高的前0.1%像素的rgb值去分别估计三通道的全局大气光

t(x)未被散射光结合暗通道先验公式如下去估计。

因此即可求待求图J(x)，为了控制误差，对t(x)做一个最小值约束，计算公式如下：

三、代码

mydehaze.m

%%
% 未被散射系数t估计
% 2018年12月29日
% 使用方法：
% Yuquan Campus, bamboopu
% 待改进功能：

%%
% clean
close all;
clear;
clc;

%%
% 读图显示图片
a_rgb=im2double(imread(‘tiananmen1.bmp‘));
% a_rgb=im2double(imread(‘stone.jpg‘));
a_gray=rgb2gray(a_rgb);
[height,weight,scale]=size(a_rgb);
% 转hsv
a_hsv=rgb2hsv(a_rgb);
% 暗通道,min(x,[],3)取每个像素位置最小值
c_minrgb=min(a_rgb,[],3);% 先取RGB最小
radius=15;
c_dark=ordfilt2(c_minrgb,1,ones(radius,radius),‘symmetric‘);%窗口最小值滤波
% 暗通道灰度对比
imwrite(c_dark,‘c_dark.bmp‘);
fig=figure;
h1=subplot(1,2,1);
imshow(c_dark);
title(‘The dark channel image‘);
h2=subplot(1,2,2);
imshow(a_gray);
title(‘The grayscale image‘);
% 全局大气光估计(采用最高亮度值)
k_find=0.001;
k_pixel=fix(k_find*height*weight);
k_sort=sort(c_dark(:),‘descend‘);
k_threshold=k_sort(k_pixel);
a_v=a_hsv(:,:,1);
a_v(c_dark<k_threshold)=0;
% ar估计
a_rgb_max=a_rgb(:,:,1);
a_rgb_max(a_v==0)=[];
a_rgb_max=a_rgb_max(:);
ar=sum(a_rgb_max)/length(a_rgb_max);
% ag估计
a_rgb_max=a_rgb(:,:,2);
a_rgb_max(a_v==0)=[];
a_rgb_max=a_rgb_max(:);
ag=sum(a_rgb_max)/length(a_rgb_max);
% ab估计
a_rgb_max=a_rgb(:,:,3);
a_rgb_max(a_v==0)=[];
a_rgb_max=a_rgb_max(:);
ab=sum(a_rgb_max)/length(a_rgb_max);
% ar ag ab最大值修正
max_a=255/255;
if ar>max_a
    ar=max_a;
end
if ag>max_a
    ag=max_a;
end
if ab>max_a
    ab=max_a;
end
% t(x)估计
a_rgb_a=zeros(height,weight,scale);
a_rgb_a(:,:,1)=a_rgb(:,:,1)/ar;
a_rgb_a(:,:,2)=a_rgb(:,:,2)/ag;
a_rgb_a(:,:,3)=a_rgb(:,:,3)/ab;
a_rgb_a_min=min(a_rgb_a,[],3);% 先取RGB最小
w=0.95;
tx=1-w*ordfilt2(a_rgb_a_min,1,ones(radius,radius),‘symmetric‘);
imwrite(tx,‘tx.bmp‘);
% 白场过渡
t0=0.1;
tx(tx<0.1)=0.1;
%tx=guidedfilter_color(a_rgb,tx,100,0.04);
tx=guidedfilter(a_gray,tx,60,0.01);
imwrite(tx,‘tx_guided.bmp‘);
figure;imshow(tx);
% 简化版去雾图估计-rgb
dh_rgb=zeros(height,weight,scale);
dh_rgb(:,:,1)=(a_rgb(:,:,1)-ar)./tx+ar;
dh_rgb(:,:,2)=(a_rgb(:,:,2)-ag)./tx+ag;
dh_rgb(:,:,3)=(a_rgb(:,:,3)-ab)./tx+ab;
figure;
h1=subplot(1,2,1);
imshow(dh_rgb);
title(‘The dehaze image‘);
h2=subplot(1,2,2);
imshow(a_rgb);
title(‘The origin image‘);
imwrite(dh_rgb,‘dh_rgb.bmp‘);

guidedfilter.m

function q = guidedfilter(I, p, r, eps)
%   GUIDEDFILTER   O(1) time implementation of guided filter.
%
%   - guidance image: I (should be a gray-scale/single channel image)
%   - filtering input image: p (should be a gray-scale/single channel image)
%   - local window radius: r
%   - regularization parameter: eps

[hei, wid] = size(I);
N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.

mean_I = boxfilter(I, r) ./ N;
mean_p = boxfilter(p, r) ./ N;
mean_Ip = boxfilter(I.*p, r) ./ N;
cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.

mean_II = boxfilter(I.*I, r) ./ N;
var_I = mean_II - mean_I .* mean_I;

a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
b = mean_p - a .* mean_I; % Eqn. (6) in the paper;

mean_a = boxfilter(a, r) ./ N;
mean_b = boxfilter(b, r) ./ N;

q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
end

四、效果