function dehaze

NS = 64 ;
N = NS * NS ;

A = [0.8 0.8 0.9]' ;

R = [0.8 0.4 0.8]' ;

Rt = R - A * (R' * A) / (norm(A)^2);

true_eta = A' * R / (norm(A) * norm(Rt)) ;

%
% Generating synthetic transmission & shading functions
%

i=0;
for y=1:NS
    for x=1:NS
    i=i+1 ;
    
    xx = x / NS - 0.5 ;
    yy = y / NS - 0.5 ;

    rr = norm([xx yy]) ;
    
    t(i) = 0.3 + 0.2 * sin(rr * 60) + 0.1 * rand ;
    l(i) = 1 + 0.3 * sin((xx+yy) * 80) + 0.1 * rand ;
    
    end
end

t=t';
l=l';
    

%
% Generating synthetic input image
%

i=0;
for y=1:NS
    for x=1:NS
    i=i+1 ;

    I(i,:) = t(i) .* l(i) .* R + (1-t(i)) .* A ;
    
    im(x,y,1) = I(i,1) ; im(x,y,2) = I(i,2) ; im(x,y,3) = I(i,3) ;
    tim(x,y,1) = l(i) .* R(1) ; tim(x,y,2) = l(i) .* R(2) ; tim(x,y,3) = l(i) .* R(3) ; 
    end
end

% estimating eta, trans. (and shading)

[est_t est_l est_eta] = estimate(A,I) ;

est_t = reshape(est_t,NS,NS) ;
est_l = reshape(est_l,NS,NS) ;


t = reshape(t,NS,NS) ;
l = reshape(l,NS,NS) ;


%
% Generating output haze-free image
%

i=0;
for y=1:NS
    for x=1:NS
    i=i+1 ;

    J(i,:) = (I(i,:) - (1-est_t(i)) * A') / est_t(i) ;
    
    oim(x,y,1) = J(i,1) ; oim(x,y,2) = J(i,2) ; oim(x,y,3) = J(i,3) ;
    end
end

figure(1)
clf

ai = ones(3,3,3) ;
ai(2,2,:) = A ;
ai = uint8(ai*255) ; 
subplot(2,4,5), imagesc(ai), title(sprintf('True eta = %f \nEst. eta = %f\nTrue A=%.2f,%.2f,%.2f\n\nAirlight color',true_eta,est_eta,A(1),A(2),A(3))) ;

colormap(gray(256))
subplot(2,4,3), imagesc(est_t), title('Est. trans.') ;
subplot(2,4,4), imagesc(est_l), title('Est. shading') ;
subplot(2,4,7), imagesc(t), title('True trans.') ;
subplot(2,4,8), imagesc(l), title('True shading') ;

im = uint8(im*255) ;
subplot(2,4,1), imshow(im), title('Input image') ;
oim = uint8(oim*255) ;
subplot(2,4,2), imshow(oim), title('Output image') ;
tim = uint8(tim*255) ;
subplot(2,4,6), imshow(tim), title('True image') ;
%
% Here is the decorrelation procedure that computes eta, and then trans.
% (t) and the shading (l)
%

function [est_t est_l est_eta] = estimate(A, I)

IA = I * A / norm(A) ;
IR = sqrt(sum(I.^2,2) - IA.^2) ;
H = (norm(A) - IA) ./ IR ;

C = cov(IA, H) ./ cov(IR, H) ;
est_eta = C(1,2) ;

est_t = 1 - (IA - est_eta * IR) / norm(A) ;
est_l = IR ./ est_t ;