function A = airlight(RedImg, GreImg, BluImg)
% estimating the global airlight by channal separately

thres = 50;
score = 0;
SquareSum_temp = Inf;
[height, width] = size(RedImg);

if max(height, width) <= thres
    for i = 1:height
        for j = 1:width
            SquareSum = (RedImg(i,j) - 255).^2 + (GreImg(i,j) - 255).^2 + (BluImg(i,j) - 255).^2;
            if SquareSum < SquareSum_temp
                SquareSum_temp = SquareSum;
                i_temp = i;
                j_temp = j;
            end
        end
    end
    A = [RedImg(i_temp,j_temp), GreImg(i_temp,j_temp), BluImg(i_temp,j_temp)];
else
    Image_R = RedImg(1:floor(height / 2),1:floor(width / 2));
    Image_G = GreImg(1:floor(height / 2),1:floor(width / 2));
    Image_B = BluImg(1:floor(height / 2),1:floor(width / 2));
    num = floor(height / 2) * floor(width / 2);
    average_R = mean(mean(Image_R));
    deviation_R = sqrt(sum(sum((Image_R - average_R).^2)) / num);
    average_G = mean(mean(Image_G));
    deviation_G = sqrt(sum(sum((Image_G - average_G).^2)) / num);
    average_B = mean(mean(Image_B));
    deviation_B = sqrt(sum(sum((Image_B - average_B).^2)) / num);
    if average_R + average_G + average_B - deviation_R - deviation_G - deviation_B > score
        score = average_R + average_G + average_B - deviation_R - deviation_G - deviation_B;
        Image_temp_R = Image_R;
        Image_temp_G = Image_G;
        Image_temp_B = Image_B;
    end

    Image_R = RedImg(1:floor(height / 2),floor(width / 2) + 1:width);
    Image_G = GreImg(1:floor(height / 2),floor(width / 2) + 1:width);
    Image_B = BluImg(1:floor(height / 2),floor(width / 2) + 1:width);
    num = floor(height / 2) * (width - floor(width / 2));
    average_R = mean(mean(Image_R));
    deviation_R = sqrt(sum(sum((Image_R - average_R).^2)) / num);
    average_G = mean(mean(Image_G));
    deviation_G = sqrt(sum(sum((Image_G - average_G).^2)) / num);
    average_B = mean(mean(Image_B));
    deviation_B = sqrt(sum(sum((Image_B - average_B).^2)) / num);
    if average_R + average_G + average_B - deviation_R - deviation_G - deviation_B > score
        score = average_R + average_G + average_B - deviation_R - deviation_G - deviation_B;
        Image_temp_R = Image_R;
        Image_temp_G = Image_G;
        Image_temp_B = Image_B;
    end

    Image_R = RedImg(floor(height / 2) + 1:height,1:floor(width / 2));
    Image_G = GreImg(floor(height / 2) + 1:height,1:floor(width / 2));
    Image_B = BluImg(floor(height / 2) + 1:height,1:floor(width / 2));
    num = (height - floor(height / 2)) * floor(width / 2);
    average_R = mean(mean(Image_R));
    deviation_R = sqrt(sum(sum((Image_R - average_R).^2)) / num);
    average_G = mean(mean(Image_G));
    deviation_G = sqrt(sum(sum((Image_G - average_G).^2)) / num);
    average_B = mean(mean(Image_B));
    deviation_B = sqrt(sum(sum((Image_B - average_B).^2)) / num);
    if average_R + average_G + average_B - deviation_R - deviation_G - deviation_B > score
        score = average_R + average_G + average_B - deviation_R - deviation_G - deviation_B;
        Image_temp_R = Image_R;
        Image_temp_G = Image_G;
        Image_temp_B = Image_B;
    end

    Image_R = RedImg(floor(height / 2) + 1:height,floor(width / 2) + 1:width);
    Image_G = GreImg(floor(height / 2) + 1:height,floor(width / 2) + 1:width);
    Image_B = BluImg(floor(height / 2) + 1:height,floor(width / 2) + 1:width);
    num = (height - floor(height / 2)) * (width - floor(width / 2));
    average_R = mean(mean(Image_R));
    deviation_R = sqrt(sum(sum((Image_R - average_R).^2)) / num);
    average_G = mean(mean(Image_G));
    deviation_G = sqrt(sum(sum((Image_G - average_G).^2)) / num);
    average_B = mean(mean(Image_B));
    deviation_B = sqrt(sum(sum((Image_B - average_B).^2)) / num);
    if average_R + average_G + average_B - deviation_R - deviation_G - deviation_B > score
        score = average_R + average_G + average_B - deviation_R - deviation_G - deviation_B;
        Image_temp_R = Image_R;
        Image_temp_G = Image_G;
        Image_temp_B = Image_B;
    end
    
    A = airlight(Image_temp_R, Image_temp_G, Image_temp_B);
end

end

