Matlab code for HW1

function oc = colorconv(ic, itype, otype)
% OC = COLORCONV(IC, ITYPE, OTYPE)
% Performs colorspace conversion.  Transforms input color IC 
% into output color OC.
% The type of conversion is specified by the second and third
% parameters.  For now, only the default of RGB -> U*V*W* is 
% supported.

% Adam Berenzweig
% 2/03/01

% default types
itype = 'RGB';
otype = 'U*V*W*';

% force ic to colors down columns
[M,N] = size(ic);
if N == 3
   ic = ic';
elseif M == 3
   % no op
else
   error('No dimension is length 3. quitting');
end

if ~strncmp(upper(itype),'RGB',3)
   % convert to RGB
end

% go to CIE XYZ:
A = [0.4900    0.3100    0.2000; ...
     0.1770    0.8130    0.0110; ... 
          0    0.0100    0.9900];
% why do the rows sum to [1 1.001 1] and not [1 1 1]?  shouldn't Y (luminance) for
% reference white be 1??    
    
% I found a slightly different transform on the web: (what's RGB709)?
%[ X ] [ 0.412453  0.35758   0.180423 ] [ R709 ] 
%[ Y ]=[ 0.212671  0.71516   0.072169 ]*[ G709 ] 
%[ Z ] [ 0.019334  0.119193  0.950227 ] [ B709 ] 
   
XYZ = A * ic;
if strncmp(otype,'XYZ',3)
   oc = XYZ
   return
end

% go to CIE UCS (u,v,Y):
%X = oc(1); Y = oc(2); Z = oc(3);
% should we return tristimulus values or u,v,y??
%if strncmp(otype,'UVW',3)
%   oc = [(2*X/3) Y ((-X + 3*Y + Z)/2)];
%   return
%end
denom = [1 15 3]*XYZ + eps;	% (X + 15*Y + 3*Z)
%u = 4*X / (X + 15*Y + 3*Z);
%v = 6*Y / (X + 15*Y + 3*Z);  %I've seen this as 9*Y/(.) on the web!! [http://www.inforamp.net/~poynton/notes/colour_and_gamma/ColorFAQ.html]
uvY = [4*XYZ(1,:)./denom; 6*XYZ(2,:)./denom; XYZ(2,:)];
oc = uvY;

% calculate chromaticity of reference white:
white = sum(A');	% RGB -> XYZ
white = white ./ sum(white);	x = white(1); y = white(2); % XYZ -> xyz
u0 = 4*x / (-2*x + 12*y + 3);
v0 = 6*y / (-2*x + 12*y + 3);
%u0 = 0.2104; v0 = 0.3159;	% chromaticity of reference white

% go to U*V*W*
% this may be wrong for luminance < .01 (e.g. black)
Ws = 25 * (100 * uvY(3,:)) .^ (1/3) - 17;  % I've also seen "lightness" L* = 116*(Y/Yn)^(1/3) - 16 [same colorFAQ above]
UsVsWs = [13 * Ws .* (uvY(1,:) - u0); 13 * Ws .* (uvY(2,:) - v0); Ws];

oc = UsVsWs;


function D = distances(P)
% D = DISTANCES(P)
% computes distances between points in P.
% The columns of P are points in an M-dimensional euclidean space,
% where MxN is the size of P.
% DISTANCES returns the NxN upper triangular matrix D where Dij is
% the Euclidean distance between points Pi and Pj.

% Adam Berenzweig
% February 4, 2001

[M,N] = size(P);

for i = 1:N
   for j = i:N
      D(i,j) = euclid_dist( P(:,i), P(:,j) );
   end
end

function d = euclid_dist(n,m)
% helper function for distances

d = sqrt(sum((n-m).^2));