function cielsh2 (fname)
clc;
fname           = 'IMG_7912s';
fname_in        = strcat (fname, '.jpg');
fname_in_Y      = strcat (fname, '_Y.jpg');
fname_LCH_out   = strcat (fname, '_blend_LCH.png');
fname_LSH_out   = strcat (fname, '_blend_LSH.png');
fname_exp_out   = strcat (fname, '_blend_exp.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

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

% 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

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


% calc Y ratio
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);

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);

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

tic;
% figure;
I_in_LCH = applycform(applycform(I_in, makecform('srgb2lab')), makecform('lab2lch'));

LCH_L = I_in_LCH (:,:,1);
LCH_c = I_in_LCH (:,:,2);
LCH_h = I_in_LCH (:,:,3);

% LCH2LSH
CIE_sat = max (0, min (1, LCH_c ./ sqrt (LCH_L.*LCH_L + LCH_c.*LCH_c)));

I_in_LCH = applycform(applycform(I_in_c2g, makecform('srgb2lab')), makecform('lab2lch'));

I_in_LCH (:,:,2) = LCH_c;
I_in_LCH (:,:,3) = LCH_h;

I_out_LCH = applycform(applycform(I_in_LCH, makecform('lch2lab')), makecform('lab2srgb'));

I_in_LCH (:,:,2) = min (I_in_LCH (:,:,1) .* CIE_sat ./ sqrt (1 - (CIE_sat .* CIE_sat)), 100);

I_out_LSH = applycform(applycform(I_in_LCH, makecform('lch2lab')), makecform('lab2srgb'));

fprintf('LCH %6.3f sec\n', toc);
tic;
figure;
subplot (2,2,1);
imshow (I_out_LCH);
subplot (2,2,2);
imshow (I_out_LSH);
subplot (2,2,3);
imshow (I_in);
subplot (2,2,4);
imshow (I_out_exposure);
fprintf('plotting %6.3f sec\n', toc); tic;

imwrite (I_out_LCH, fname_LCH_out);
imwrite (I_out_LSH, fname_LSH_out);
imwrite (im2uint8 (I_out_exposure), fname_exp_out);
fprintf('saving %6.3f sec\n', toc);

end
