function blend_luminance (fname)

% function un-gamma
function image_linear = ilinearize (i)
    image_linear = power(((i + 0.055) / 1.055), 2.4);
    image_linear (i < 0.04045) = i (i < 0.04045) / 12.92;
end

% function gamma
function image_gamma = igamma (i)
    image_gamma = 1.055 * power(i, 1/2.4) - 0.055;
    image_gamma (i <= 0.0031308) = i (i <= 0.0031308) * 12.92;
end

clc;

%fname           = 'IMG_2909';
fname_in        = strcat (fname, '.JPG');
fname_in_Y      = strcat (fname, '_Y.JPG');
fname_exp_out   = strcat (fname, '_exp_out.png');


I_in = im2double (imread(fname_in));
I_in_c2g = im2double (imread(fname_in_Y));
% convert if gray scale image
if (size (I_in_c2g, 3) == 1)
    I_in_c2g (:,:,2) = I_in_c2g (:,:,1);
    I_in_c2g (:,:,3) = I_in_c2g (:,:,1);
end

% prealloc
I_in_linear = I_in;
I_in_c2g_linear = I_in_c2g;
Y_in = I_in (:,:,2);
Y_in_c2g = I_in_c2g (:,:,2);

tic;
I_in_linear = ilinearize (I_in);
fprintf('degamma in %6.3f sec\n', toc); fflush(stdout);
tic;
I_in_c2g_linear = ilinearize (I_in_c2g);
fprintf('degamma c2g %6.3f sec\n', toc); fflush(stdout);


% calc Y ratio
tic
Y_in = 0.2126*I_in_linear (:,:,1) + 0.7152*I_in_linear (:,:,2) + 0.0722*I_in_linear (:,:,3);
Y_in_c2g = 0.2126*I_in_c2g_linear (:,:,1) + 0.7152*I_in_c2g_linear (:,:,2) + 0.0722*I_in_c2g_linear (:,:,3);

Y_ratio = Y_in_c2g ./ max (Y_in, 1e-9);

fprintf('ratio in %6.3f sec\n', toc); fflush(stdout);


tic;
I_out_linear (:,:,1) = I_in_linear (:,:,1) .* Y_ratio;
I_out_linear (:,:,2) = I_in_linear (:,:,2) .* Y_ratio;
I_out_linear (:,:,3) = I_in_linear (:,:,3) .* Y_ratio;
fprintf('mul %6.3f sec\n', toc); fflush(stdout);

tic
%scaling_max = max (max (max (I_out_linear)))
%I_out_linear = I_out_linear / scaling_max;
I_out_exposure = igamma (I_out_linear);
fprintf('gamma %6.3f sec\n', toc); fflush(stdout);

tic;
%figure;
imshow (I_out_exposure);
fprintf('plotting %6.3f sec\n', toc); fflush(stdout);

tic;
imwrite (im2uint8 (I_out_exposure), fname_exp_out);
fprintf('saving %6.3f sec\n', toc);

end
