% function [lab] = xyz2cielab(xyz, ref_xyz)
%
% // C-Code
% var_X = X /  95.047          //Observer = 2°, Illuminant = D65
% var_Y = Y / 100.000
% var_Z = Z / 108.883
%
% if ( var_X > 0.008856 ) var_X = var_X ^ ( 1/3 )
% else                    var_X = ( 7.787 * var_X ) + ( 16 / 116 )
% if ( var_Y > 0.008856 ) var_Y = var_Y ^ ( 1/3 )
% else                    var_Y = ( 7.787 * var_Y ) + ( 16 / 116 )
% if ( var_Z > 0.008856 ) var_Z = var_Z ^ ( 1/3 )
% else                    var_Z = ( 7.787 * var_Z ) + ( 16 / 116 )
%
% CIE-L* = ( 116 * var_Y ) - 16
% CIE-a* = 500 * ( var_X - var_Y )
% CIE-b* = 200 * ( var_Y - var_Z )
%
% Source http://www.easyrgb.com

function [lab] = xyz2cielab(xyz, ref_xyz)

[m,n] = size(xyz);
if ((m ~= 3) & (n == 3))
    xyz = xyz'
end;

xyz(1,:) = xyz(1,:)./ref_xyz(1,:);
xyz(2,:) = xyz(2,:)./ref_xyz(2,:);
xyz(3,:) = xyz(3,:)./ref_xyz(3,:);
xyz_k = f1(xyz);

% L
lab(1,:) = 116*xyz_k(2,:) - 16;
% a
lab(2,:) = 500*(xyz_k(1,:) - xyz_k(2,:));
% b
lab(3,:) = 200*(xyz_k(2,:) - xyz_k(3,:));

function xyz_k = f1(xyz)
thresh_n = 0.008856;

[idx_m, idx_n] = find(xyz > thresh_n);
xyz_k(idx_m,idx_n) = xyz(idx_m,idx_n).^(1/3);

[idx_m, idx_n] = find(xyz <= thresh_n);
xyz_k(idx_m,idx_n) = 7.787 * xyz(idx_m,idx_n) + 16/116;