clear all
close all
clc

imageName='sam_3.bmp';
image=imread(imageName);
sz=size(image);
height=sz(1);
width=sz(2);

image=double(image);
figure,imshow(uint8(image)),title('original');

radius=10;
p=0.95;
minAtmosLight=240;

W=zeros(height,width);

for i=1:height
    for j=1:width
         W(i,j)=min(image(i,j,:));
    end
end

% figure,title('W(x,y)'),imshow(uint8(W));

A=zeros(height,width);

A=medfilt2(W,[radius,radius],'symmetric');

% figure,imshow(uint8(A)),title('A(x,y)');

B=A-medfilt2(abs(W-A),[radius,radius],'symmetric');

V=max(min(p*B,W),0);
figure,imshow(uint8(V)),title('V');

% estimate the atmospheric light,Using He's method in
% ex_darkchannel_guildedfilter.m
ALight = min([minAtmosLight, max(max(255-W))])

image_gray=double(rgb2gray(image));
R_gray=(image_gray-V)./(ones(height,width)-V/ALight);
figure,imshow(uint8(R_gray)),title('R\_gray');

V1=zeros(height,width,3);
for i=1:height
    for j=1:width
        V1(i,j,:)=V(i,j);
    end
end

R=(image-V1)./(1-V1./ALight);

figure,imshow(uint8(R)),title('R(x,y)');
% figure,imhist(rgb2gray(uint8(R))),title('imhist of R');

gamma=1.3;
R=R.^gamma;
figure,imshow(uint8(R)),title('after gamma');