ceshi

% void CLASS cam_xyz_coeff (float rgb_cam[3][4], double cam_xyz[4][3])
% {
% double cam_rgb[4][3], inverse[4][3], num;
% int i, j, k;
%
% for (i=0; i < colors; i++) /* Multiply out XYZ colorspace /
% for (j=0; j < 3; j++)
% for (cam_rgb[i][j] = k=0; k < 3; k++)
% cam_rgb[i][j] += cam_xyz[i][k] * xyz_rgb[k][j];
%
% for (i=0; i < colors; i++) { / Normalize cam_rgb so that /
% for (num=j=0; j < 3; j++) / cam_rgb * (1,1,1) is (1,1,1,1) */
% num += cam_rgb[i][j];
% for (j=0; j < 3; j++)
% cam_rgb[i][j] /= num;
% pre_mul[i] = 1 / num;
% }
% pseudoinverse (cam_rgb, inverse, colors);
% for (i=0; i < 3; i++)
% for (j=0; j < colors; j++)
% rgb_cam[i][j] = inverse[j][i];
% }

% void CLASS pseudoinverse (double (in)[3], double (out)[3], int size)
% {
% double work[3][6], num;
% int i, j, k;
%
% for (i=0; i < 3; i++) {
% for (j=0; j < 6; j++)
% work[i][j] = j == i+3;
% for (j=0; j < 3; j++)
% for (k=0; k < size; k++)
% work[i][j] += in[k][i] * in[k][j];
% }
% for (i=0; i < 3; i++) {
% num = work[i][i];
% for (j=0; j < 6; j++)
% work[i][j] /= num;
% for (k=0; k < 3; k++) {
% if (k==i) continue;
% num = work[k][i];
% for (j=0; j < 6; j++)
% work[k][j] -= work[i][j] * num;
% }
% }
% for (i=0; i < size; i++)
% for (j=0; j < 3; j++)
% for (out[i][j]=k=0; k < 3; k++)
% out[i][j] += work[j][k+3] * in[i][k];
% }

% const double xyz_rgb[3][3] = { /* XYZ from RGB */
% { 0.412453, 0.357580, 0.180423 },
% { 0.212671, 0.715160, 0.072169 },
% { 0.019334, 0.119193, 0.950227 } };
clc
clear all
cmsrc0 = [1.895316, -0.688256,-0.20706;
-0.180953,1.625399,-0.444446;
-0.0000147,-0.524012,1.524159;
0,0,0]

cam_xyz = [0.617188, -0.140625, -0.0390625;
-0.460938, 1.19531, 0.226563;
-0.0859375, 0.1953113, 0.545875;
0,0,0
];

cam_rgb = zeros(4,3);
inverse = zeros(4,3);
xyz_rgb = [0.412453, 0.357580, 0.180423,0;
0.212671, 0.715160, 0.072169,0;
0.019334, 0.119193, 0.950227,0];
cam_rgb = cam_xyz * xyz_rgb;
for i =1:3
cam_rgb(i,:) = cam_rgb(i,:) / sum(cam_rgb(i,:));
end
cam_rgb = cam_rgb(:,1:3);
colors = 3;

x11 = cam_rgb(1,1);
x12 = cam_rgb(1,2);
x13 = cam_rgb(1,3);
x21 = cam_rgb(2,1);
x22 = cam_rgb(2,2);
x23 = cam_rgb(2,3);
x31 = cam_rgb(3,1);
x32 = cam_rgb(3,2);
x33 = cam_rgb(3,3);
work = zeros(3,6);
work(1,1) = x11 *x11 + x21 *x21 + x31 * x31;
work(1,2) = x11 *x12 + x21 *x22 + x31 * x32;
work(1,3) = x11 *x13 + x21 *x23 + x31 * x33;
work(2,1) = x11 *x12 + x21 *x22 + x31 * x32;
work(2,2) = x12 *x12 + x22 *x22 + x32 * x32;
work(2,3) = x12 *x13 + x22 *x23 + x32 * x33;
work(3,1) = x11 *x13 + x21 *x23 + x31 * x33;
work(3,2) = x12 *x13 + x22 *x23 + x32 * x33;
work(3,3) = x13 *x13 + x23 *x23 + x33 * x33;
work(:,4:6) = [1,0,0;0,1,0;0,0,1];

% i =1
num1 = work(1,1);
work(1,1) = work(1,1) / num1;
work(1,2) = work(1,2) / num1;
work(1,3) = work(1,3) / num1;
work(1,4) = work(1,4) / num1;
work(1,5) = work(1,5) / num1;
work(1,6) = work(1,6) / num1;

% k = 2
num2 = work(2,1);
work(2,1) = work(2,1) - num2 * work(1,1);
work(2,2) = work(2,2) - num2 * work(1,2);
work(2,3) = work(2,3) - num2 * work(1,3);
work(2,4) = work(2,4) - num2 * work(1,4);
work(2,5) = work(2,5) - num2 * work(1,5);
work(2,6) = work(2,6) - num2 * work(1,6);

% k = 3
num3 = work(3,1);
work(3,1) = work(3,1) - num3 * work(1,1);
work(3,2) = work(3,2) - num3 * work(1,2);
work(3,3) = work(3,3) - num3 * work(1,3);
work(3,4) = work(3,4) - num3 * work(1,4);
work(3,5) = work(3,5) - num3 * work(1,5);
work(3,6) = work(3,6) - num3 * work(1,6);

% i = 2
num1 = work(2,2);
work(2,1) = work(2,1) / num1;
work(2,2) = work(2,2) / num1;
work(2,3) = work(2,3) / num1;
work(2,4) = work(2,4) / num1;
work(2,5) = work(2,5) / num1;
work(2,6) = work(2,6) / num1;

num2 = work(1,2);
work(1,1) = work(1,1) - num2 * work(2,1);
work(1,2) = work(1,2) - num2 * work(2,2);
work(1,3) = work(1,3) - num2 * work(2,3);
work(1,4) = work(1,4) - num2 * work(2,4);
work(1,5) = work(1,5) - num2 * work(2,5);
work(1,6) = work(1,6) - num2 * work(2,6);

num3 = work(3,2);
work(3,1) = work(3,1) - num3 * work(2,1);
work(3,2) = work(3,2) - num3 * work(2,2);
work(3,3) = work(3,3) - num3 * work(2,3);
work(3,4) = work(3,4) - num3 * work(2,4);
work(3,5) = work(3,5) - num3 * work(2,5);
work(3,6) = work(3,6) - num3 * work(2,6);

% i = 3
num1 = work(3,3);
work(3,1) = work(3,1) / num1;
work(3,2) = work(3,2) / num1;
work(3,3) = work(3,3) / num1;
work(3,4) = work(3,4) / num1;
work(3,5) = work(3,5) / num1;
work(3,6) = work(3,6) / num1;

num2 = work(1,3);
work(1,1) = work(1,1) - num2 * work(3,1);
work(1,2) = work(1,2) - num2 * work(3,2);
work(1,3) = work(1,3) - num2 * work(3,3);
work(1,4) = work(1,4) - num2 * work(3,4);
work(1,5) = work(1,5) - num2 * work(3,5);
work(1,6) = work(1,6) - num2 * work(3,6);

num3 = work(2,3);
work(2,1) = work(2,1) - num3 * work(3,1);
work(2,2) = work(2,2) - num3 * work(3,2);
work(2,3) = work(2,3) - num3 * work(3,3);
work(2,4) = work(2,4) - num3 * work(3,4);
work(2,5) = work(2,5) - num3 * work(3,5);
work(2,6) = work(2,6) - num3 * work(3,6);

out = zeros(3,3);
out(1,1) = work(1,4) * cam_rgb(1,1) + work(1,5) * cam_rgb(1,2) + work(1,6) * cam_rgb(1,3);
out(1,2) = work(2,4) * cam_rgb(1,1) + work(2,5) * cam_rgb(1,2) + work(2,6) * cam_rgb(1,3);
out(1,3) = work(3,4) * cam_rgb(1,1) + work(3,5) * cam_rgb(1,2) + work(3,6) * cam_rgb(1,3);
out(2,1) = work(1,4) * cam_rgb(2,1) + work(1,5) * cam_rgb(2,2) + work(1,6) * cam_rgb(2,3);
out(2,2) = work(2,4) * cam_rgb(2,1) + work(2,5) * cam_rgb(2,2) + work(2,6) * cam_rgb(2,3);
out(2,3) = work(3,4) * cam_rgb(2,1) + work(3,5) * cam_rgb(2,2) + work(3,6) * cam_rgb(2,3);
out(3,1) = work(1,4) * cam_rgb(3,1) + work(1,5) * cam_rgb(3,2) + work(1,6) * cam_rgb(3,3);
out(3,2) = work(2,4) * cam_rgb(3,1) + work(2,5) * cam_rgb(3,2) + work(2,6) * cam_rgb(3,3);
out(3,3) = work(3,4) * cam_rgb(3,1) + work(3,5) * cam_rgb(3,2) + work(3,6) * cam_rgb(3,3);
out'

% work = zeros(3,6);
% for i =1:3
% for j = 1:6
% if j == i+3
% work(i,j) = 1;
% else
% work(i,j) = 0;
% end
% end
%
% for j = 1:3
% for k = 1:colors
% work(i,j) = work(i,j) + cam_rgb(k,i) * cam_rgb(k,j);
% end
% end
% end
%
% for i =1:3
% num = work(i,i);
% for j = 1:6
% work(i,j) = work(i,j) / num;
% end
% for k = 1:3
% if(k == i)
% continue;
% end
% num = work(k,i);
% for j = 1:6
% work(k,j) = work(k,j) - work(i,j) * num;
% end
% end
% end

% out = zeros(3,3);
% for i = 1:colors
% for j =1:3
% out(i,j) =0;
% for k = 1:3
% out(i,j) = out(i,j) + work(j,k+3) * cam_rgb(i,k);
% end
% end
% end
%
% out = out'