% function [xyz] = rgb2xyz(rgb)
% 
% // C-Code
% if ( var_R > 0.04045 ) var_R = ( ( var_R + 0.055 ) / 1.055 ) ^ 2.4
% else                   var_R = var_R / 12.92
% if ( var_G > 0.04045 ) var_G = ( ( var_G + 0.055 ) / 1.055 ) ^ 2.4
% else                   var_G = var_G / 12.92
% if ( var_B > 0.04045 ) var_B = ( ( var_B + 0.055 ) / 1.055 ) ^ 2.4
% else                   var_B = var_B / 12.92
%
% var_R = var_R * 100
% var_G = var_G * 100
% var_B = var_B * 100
%
% // Observer. = 2°, Illuminant = D65
% X = var_R * 0.4124 + var_G * 0.3576 + var_B * 0.1805
% Y = var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722
% Z = var_R * 0.0193 + var_G * 0.1192 + var_B * 0.9505
%
% Source http://www.easyrgb.com

function [xyz] = rgb2xyz(rgb)
[m,n] = size(rgb);

if ((m ~= 3) & (n == 3))
    rgb = rgb';
end;

% Observer = 2°, Illuminant = D65
t_mat = [[ 0.412453  0.357580  0.180423 ]
         [ 0.212671  0.715160  0.072169 ]
         [ 0.019334  0.119193  0.950227 ]];

rgb_k = f1(rgb);
xyz = t_mat * rgb_k;

function rgb_k = f1(rgb)
thresh_n = 0.04045 ;

[idx_m, idx_n] = find(rgb > thresh_n);
rgb_k(idx_m,idx_n) = ((rgb(idx_m,idx_n) + 0.055)/1.055).^(2.4);

[idx_m, idx_n] = find(rgb <= thresh_n);
rgb_k(idx_m,idx_n) = rgb(idx_m,idx_n)/12.92;