function t = transmission(GrayImg, A)
% estimate optimal transmission map

[height, width] = size(GrayImg);
win_height = floor(height*0.01);
win_width = floor(width*0.01);
height_num = floor(height / win_height);
width_num = floor(width / win_width);

a = 100;
b = 1;
c = 5;
t = ones(height, width);

win_num_1 = win_height * win_width;
for i = 0 : height_num - 1
    for j = 0 : width_num - 1
        win = GrayImg(i * win_height + 1:(i + 1) * win_height, j * win_width + 1:(j + 1) * win_width);
        t_op = optimal(win, win_num_1, A, a, b, c);
        t(i * win_height + 1:(i + 1) * win_height, j * win_width + 1:(j + 1) * win_width) = t_op * ones(win_height, win_width);
    end
end

win_num_2 = win_height * (width - width_num * win_width);
for i = 0 : height_num - 1
    win = GrayImg(i * win_height + 1:(i + 1) * win_height, width_num * win_width + 1:width);
    t_op = optimal(win, win_num_2, A, a, b, c);
    t(i * win_height + 1:(i + 1) * win_height, width_num * win_width + 1 : width) = t_op * ones(win_height, width - width_num * win_width);
end

win_num_3 = (height - height_num * win_height) * win_width;
for j = 0 : width_num - 1
    win = GrayImg(height_num * win_height + 1:height, j * win_width + 1:(j + 1) * win_width);
    t_op = optimal(win, win_num_3, A, a, b, c);
    t(height_num * win_height + 1:height, j * win_width + 1:(j + 1) * win_width) = t_op * ones(height - height_num * win_height, win_width);
end

win_num_4 = (height - height_num * win_height) * (width - width_num * win_width);
win = GrayImg(height_num * win_height + 1:height, width_num * win_width + 1:width);
t_op = optimal(win, win_num_4, A, a, b, c);
t(height_num * win_height + 1:height, width_num * win_width + 1:width) = t_op * ones(height - height_num * win_height, width - width_num * win_width);

end

function t_op = optimal(win, win_num, A, a, b, c)

t_op = 0;
f_min = Inf;

for t = 6/32 : 1/32 : 31/32
	win_R = (win - A) / t  + A;
    
	% f_std
	win_mean = mean(mean(win_R));
	f_std = sqrt(sum(sum((win_R - win_mean).^2)) / win_num);
	
    % f_unit
	[counts_R, x] = imhist(uint8(win_R));
	f_unit = sqrt(sum((1 / 256 - (counts_R / win_num)).^2));
	
    % f_loss
	alpha = ceil(A * (1 - t));
	alpha_x = ones * (0:alpha - 1);
	beta = floor(255 * t + A * (1 - t));
	beta_x = ones * (beta + 1:255);
	[counts, x] = imhist(uint8(win));
	f_loss = sum((counts(1:alpha))' .* (((alpha_x - A) ./ t + A).^2))+ sum((counts(beta + 2:256))' .* (((beta_x - A - 1) ./ t + A - 255).^2));

	f = a * f_unit - b * f_std + c * log(f_loss);
    % f = sigma*(a * f_unit - b * f_std) + c * log(f_loss);
    
	if f < f_min
		f_min = f;
		t_op = t;
	end
end

end

