function [C, M, Y] = xyz2cmy (point, space)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% xyz2cmy.m                               %
% Written 3/1/00 by Ben Olding            %
% Given a data set relating xyz values    %
% to cmy space, calculates cmy values     %
% for a point in xyz space using local    %
% linear interpolation                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Requires: findtetra.m, tetravolume.m,
%           coplanar.m

% Space is [X Y Z C M Y]
% Point is [x y z]

[a, b, c, d, V]=findtetra(point,space);

% weight vertices of surrounding tetrahedron
if V >0,
  w1=tetravolume (b(1:3),c(1:3),d(1:3),point);
  w2=tetravolume (a(1:3),c(1:3),d(1:3),point);
  w3=tetravolume (a(1:3),b(1:3),d(1:3),point);
  w4=tetravolume (a(1:3),b(1:3),c(1:3),point);
  W = [w1 w2 w3 w4]/V;
  C = W*[a(4);b(4);c(4);d(4)];
  M = W*[a(5);b(5);c(5);d(5)];
  Y = W*[a(6);b(6);c(6);d(6)];

% if point is equal to a point in data space, equate
elseif V==0,
  C = a(4);
  M = a(5);
  Y = a(6);

% if point is outside gamut, clip to nearest plane/line/point
else
  close3=[a; b; c];
  
  % create plane out of closest three points in gamut
  v1=close3(2,1:3)-close3(1,1:3);
  v2=close3(3,1:3)-close3(1,1:3);
  nrm= cross(v1,v2);
  d = nrm*close3(1,1:3)';
  n = (d-point(1)-point(2)-point(3))/(nrm(1)+nrm(2)+nrm(3));
  
  % project point onto plane using normal
  newpoint = point+n*nrm;
  ptlist = [newpoint; close3(1,1:3); close3(2,1:3); close3(3,1:3)];
  ptlist = sort(ptlist);
  
  % if new point is outside triangle of 3 points, project to line
  if ((newpoint(1)==ptlist(1,1) | newpoint(1)==ptlist(4,1)) | ...
      (newpoint(2)==ptlist(1,2) | newpoint(2)==ptlist(4,2)) | ...
      (newpoint(3)==ptlist(1,3) | newpoint(3)==ptlist(4,3))),
    
    % Project point onto nearest line
    v1=point-close3(1,1:3);
    v2=close3(2,1:3)-close3(1,1:3);
    n = dot(v1,v2);
    newpoint = close3(1,1:3)+n*v2;
    
    % if new point is outside line segment, use nearest point
    if n<0 | n>1,
      C=close3(1,4);
      M=close3(1,5);
      Y=close3(1,6);
    
    % if new point is inside line segment, weight with endpoints
    else
      p1=close3(1,:);
      p2=close3(2,:);
      np=newpoint;
      l  = sqrt((p1(1)-p2(1)).^2+(p1(2)-p2(2)).^2+(p1(3)-p2(3)).^2);
      w1 = sqrt((np(1)-p2(1)).^2+(np(2)-p2(2)).^2+(np(3)-p2(3)).^2)/l;
      w2 = sqrt((np(1)-p1(1)).^2+(np(2)-p1(2)).^2+(np(3)-p1(3)).^2)/l;
      C=w1*p1(4)+w2*p2(4);
      M=w1*p1(5)+w2*p2(5);
      Y=w1*p1(6)+w2*p2(6);    
    end;
  
  % if new point is within triangle, weight with vertices
  else
    p1=close3(1,1:3);
    p2=close3(2,1:3);
    p3=close3(3,1:3);
    np=newpoint;
    a =trianglearea(p1,p2,p3);
    w1=trianglearea(p2,p3,np)/a;
    w2=trianglearea(p1,p3,np)/a;
    w3=trianglearea(p1,p2,np)/a;
    C = w1*close3(1,4)+w2*close3(2,4)+w3*close3(3,4);
    M = w1*close3(1,5)+w2*close3(2,5)+w3*close3(3,5);
    Y = w1*close3(1,6)+w2*close3(2,6)+w3*close3(3,6);
  end;
end;      



