function t = transmission(RedImg, GreImg, BluImg, A)
% estimate optimal transmission map

[height, width] = size(RedImg);
win_height = floor(height * 0.03);
win_width = floor(width * 0.03);
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_R = RedImg(i * win_height + 1:(i + 1) * win_height, j * win_width + 1:(j + 1) * win_width);
        win_G = GreImg(i * win_height + 1:(i + 1) * win_height, j * win_width + 1:(j + 1) * win_width);
        win_B = BluImg(i * win_height + 1:(i + 1) * win_height, j * win_width + 1:(j + 1) * win_width);
        t_op = optimal(win_R, win_G, win_B, 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_R = RedImg(i * win_height + 1:(i + 1) * win_height, width_num * win_width + 1:width);
    win_G = GreImg(i * win_height + 1:(i + 1) * win_height, width_num * win_width + 1:width);
    win_B = BluImg(i * win_height + 1:(i + 1) * win_height, width_num * win_width + 1:width);
    t_op = optimal(win_R, win_G, win_B, 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_R = RedImg(height_num * win_height + 1:height, j * win_width + 1:(j + 1) * win_width);
    win_G = GreImg(height_num * win_height + 1:height, j * win_width + 1:(j + 1) * win_width); 
    win_B = BluImg(height_num * win_height + 1:height, j * win_width + 1:(j + 1) * win_width);
    t_op = optimal(win_R, win_G, win_B, 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_R = RedImg(height_num * win_height + 1:height, width_num * win_width + 1:width);
win_G = GreImg(height_num * win_height + 1:height, width_num * win_width + 1:width);
win_B = BluImg(height_num * win_height + 1:height, width_num * win_width + 1:width);
t_op = optimal(win_R, win_G, win_B, 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_R, win_G, win_B, win_num, A, a, b, c)

t_op = 0;
f_min = Inf;

for t = 6/32 : 1/32 : 31/32
	win_J_R = (win_R - A(1)) / t  + A(1);
    win_J_G = (win_G - A(2)) / t  + A(2);
    win_J_B = (win_B - A(3)) / t  + A(3);
    
	% f_std
	win_mean_R = mean(mean(win_J_R));
	f_std_R = sqrt(sum(sum((win_J_R - win_mean_R).^2)) / win_num);
    win_mean_G = mean(mean(win_J_G));
	f_std_G = sqrt(sum(sum((win_J_G - win_mean_G).^2)) / win_num);
    win_mean_B = mean(mean(win_J_B));
	f_std_B = sqrt(sum(sum((win_J_B - win_mean_B).^2)) / win_num);
	
    % f_unit
	[counts_J_R, x] = imhist(uint8(win_J_R));
	f_unit_R = sqrt(sum((1 / 256 - (counts_J_R / win_num)).^2));
    [counts_J_G, x] = imhist(uint8(win_J_G));
	f_unit_G = sqrt(sum((1 / 256 - (counts_J_G / win_num)).^2));
    [counts_J_B, x] = imhist(uint8(win_J_B));
	f_unit_B = sqrt(sum((1 / 256 - (counts_J_B / win_num)).^2));
	
    % f_loss
	alpha_R = ceil(A(1) * (1 - t));
	alpha_x_R = ones * (0:alpha_R - 1);
	beta_R = floor(255 * t + A(1) * (1 - t));
	beta_x_R = ones * (beta_R + 1:255);
	[counts_R, x] = imhist(uint8(win_R));
	f_loss_R = sum((counts_R(1:alpha_R))' .* (((alpha_x_R - A(1)) ./ t + A(1)).^2))+ sum((counts_R(beta_R + 2:256))' .* (((beta_x_R - A(1)) ./ t + A(1) - 255).^2));
    f_loss_R = max(1, f_loss_R);

    alpha_G = ceil(A(2) * (1 - t));
	alpha_x_G = ones * (0:alpha_G - 1);
	beta_G = floor(255 * t + A(2) * (1 - t));
	beta_x_G = ones * (beta_G + 1:255);
	[counts_G, x] = imhist(uint8(win_G));
	f_loss_G = sum((counts_G(1:alpha_G))' .* (((alpha_x_G - A(2)) ./ t + A(2)).^2))+ sum((counts_G(beta_G + 2:256))' .* (((beta_x_G - A(2)) ./ t + A(2) - 255).^2));
    f_loss_G = max(1, f_loss_G);
    
    alpha_B = ceil(A(3) * (1 - t));
	alpha_x_B = ones * (0:alpha_B - 1);
	beta_B = floor(255 * t + A(3) * (1 - t));
	beta_x_B = ones * (beta_B + 1:255);
	[counts_B, x] = imhist(uint8(win_B));
	f_loss_B = sum((counts_B(1:alpha_B))' .* (((alpha_x_B - A(3)) ./ t + A(3)).^2))+ sum((counts_B(beta_B + 2:256))' .* (((beta_x_B - A(3)) ./ t + A(3) - 255).^2));
    f_loss_B = max(1, f_loss_B);
	f = a * (f_unit_R + f_unit_G + f_unit_B) - b * (f_std_R + f_std_G + f_std_B) + c * log(f_loss_R * f_loss_G * f_loss_B);
    
	if f < f_min
		f_min = f;
		t_op = t;
	end
end

end

