-
Notifications
You must be signed in to change notification settings - Fork 0
/
non_local_dehazing.m
132 lines (109 loc) · 4.89 KB
/
non_local_dehazing.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
function [img_dehazed, transmission, cluster] = non_local_dehazing(img_hazy, air_light, gamma)
%The core implementation of "Non-Local Image Dehazing", CVPR 2016
%
% The details of the algorithm are described in our paper:
% Non-Local Image Dehazing. Berman, D. and Treibitz, T. and Avidan S., CVPR2016,
% which can be found at:
% www.eng.tau.ac.il/~berman/NonLocalDehazing/NonLocalDehazing_CVPR2016.pdf
% If you use this code, please cite the paper.
%
% Input arguments:
% ----------------
% img_hazy - A hazy image in the range [0,255], type: uint8
% air_light - As estimated by prior methods, normalized to the range [0,1]
% gamma - Radiometric correction. If empty, 1 is assumed
%
% Output arguments:
% ----------------
% img_dehazed - The restored radiance of the scene (uint8)
% transmission - Transmission map of the scene, in the range [0,1]
%
% Author: Dana Berman, 2016.
%
% The software code of the Non-Local Image Dehazing algorithm is provided
% under the attached LICENSE.md
%% Validate input
[h,w,n_colors] = size(img_hazy);
if (n_colors ~= 3) % input verification
error(['Non-Local Dehazing reuires an RGB image, while input ',...
'has only ',num2str(n_colors),' dimensions']);
end
if ~exist('air_light','var') || isempty(air_light) || (numel(air_light)~=3)
error('Dehazing on sphere requires an RGB airlight');
end
if ~exist('gamma','var') || isempty(gamma), gamma = 1; end
img_hazy = im2double(img_hazy);
img_hazy_corrected = img_hazy.^gamma; % radiometric correction
%% Find Haze-lines
% Translate the coordinate system to be air_light-centric (Eq. (3))
dist_from_airlight = double(zeros(h,w,n_colors));
for color_idx=1:n_colors
dist_from_airlight(:,:,color_idx) = img_hazy_corrected(:,:,color_idx) - air_light(:,:,color_idx);
end
% Calculate radius (Eq. (5))
radius = sqrt( dist_from_airlight(:,:,1).^2 + dist_from_airlight(:,:,2).^2 +dist_from_airlight(:,:,3).^2 );
% Cluster the pixels to haze-lines
% Use a KD-tree impementation for fast clustering according to their angles
dist_unit_radius = reshape(dist_from_airlight,[h*w,n_colors]);
dist_norm = sqrt(sum(dist_unit_radius.^2,2));
dist_unit_radius = bsxfun(@rdivide, dist_unit_radius, dist_norm);
n_points = 1000;
% load pre-calculated uniform tesselation of the unit-sphere
fid = fopen(['TR',num2str(n_points),'.txt']);
points = cell2mat(textscan(fid,'%f %f %f'));
fclose(fid);
mdl = KDTreeSearcher(points);
ind = knnsearch(mdl, dist_unit_radius);
cluster = reshape( ind, h, w);
%% Estimating Initial Transmission
% Estimate radius as the maximal radius in each haze-line (Eq. (11))
K = accumarray(ind, radius(:), [n_points,1], @max);
radius_new = reshape( K(ind), h, w);
% Estimate transmission as radii ratio (Eq. (12))
transmission_estimation = radius./radius_new;
% Limit the transmission to the range [trans_min, 1] for numerical stability
trans_min = 0.1;
transmission_estimation = min(max(transmission_estimation, trans_min),1);
%% Regularization
% Apply lower bound from the image (Eqs. (13-14))
trans_lower_bound = 1 - min(bsxfun(@rdivide,img_hazy_corrected,reshape(air_light,1,1,3)) ,[],3);
transmission_estimation = max(transmission_estimation, trans_lower_bound);
% for i = 1:h
% for j = 1:w
% if cluster(i,j) == 182 || 633
% transmission_estimation(i, j) = transmission_estimation(i, j) + 0.7;
% end
% end
% end
% Solve optimization problem (Eq. (15))
% find bin counts for reliability - small bins (#pixels<50) do not comply with
% the model assumptions and should be disregarded
bin_count = accumarray(ind,1,[n_points,1]);
bin_count_map = reshape(bin_count(ind),h,w);
bin_eval_fun = @(x) min(1, x/50);
% Calculate std - this is the data-term weight of Eq. (15)
K_std = accumarray(ind,radius(:),[n_points,1],@std);
radius_std = reshape( K_std(ind), h, w);
radius_eval_fun = @(r) min(1, 3*max(0.001, r-0.1));
radius_reliability = radius_eval_fun(radius_std./max(radius_std(:)));
data_term_weight = bin_eval_fun(bin_count_map).*radius_reliability;
lambda = 0.1;
transmission = wls_optimization(transmission_estimation, data_term_weight, img_hazy, lambda);
%% Dehazing
% (Eq. (16))
img_dehazed = zeros(h,w,n_colors);
leave_haze = 1.06; % leave a bit of haze for a natural look (set to 1 to reduce all haze)
for color_idx = 1:3
img_dehazed(:,:,color_idx) = ( img_hazy_corrected(:,:,color_idx) - ...
(1-leave_haze.*transmission).*air_light(color_idx) )./ max(transmission,trans_min);
end
% Limit each pixel value to the range [0, 1] (avoid numerical problems)
img_dehazed(img_dehazed>1) = 1;
img_dehazed(img_dehazed<0) = 0;
img_dehazed = img_dehazed.^(1/gamma); % radiometric correction
% For display, we perform a global linear contrast stretch on the output,
% clipping 0.5% of the pixel values both in the shadows and in the highlights
adj_percent = [0.005, 0.995];
img_dehazed = adjust(img_dehazed,adj_percent);
img_dehazed = im2uint8(img_dehazed);
end % function non_local_dehazing