% Calculate XYZ values of MacBeth Chart

load macbethChart


% Li-Yi Calculated the XYZ values for each patch under each
% illumination and stored them in a 3x24x18 matrix.

load macbethXYZ_filters

% LAB values are calculated from the matrix using xyz2lab.  The
% 19th value on the MacBeth chart is white.

for i=1:18
     macbethLAB_filters(:,:,i) = xyz2lab(macbethXYZ_filters(:,:,i), macbethXYZ_filters(:,19,i));
end

% Li-Yi extracted the RGB values from each photograph from each
% camera and stored them in a 18x3x24 matrix
% (illumationxRGBxcolor). 

load sony
load olympus

% Make these 18x3x24 matrices into 3x24x18 matrices

for i=1:18
  for j=1:3
    for k=1:24
      sony(j,k,i) = sonyInverseGammaColors(i,j,k);
      olympus(j,k,i) = olympusInverseGammaColors(i,j,k);
    end
  end
end

% scale the RGB values from 0-1 instead of 0-255

for i=1:3
  for j=1:24
    for k=1:18
      sony(i,j,k) = sony(i,j,k) / 255;
      olumpus(i,j,k) = olympus(i,j,k) / 255;
    end
  end
end

% Convert the RGB values to XYZ

load /local/class/psych221/tutorials/cMatch/monitor;
load xyz;

rgb2xyz = xyz'*phosphors;

for i=1:18
   for j=1:24
      sonyXYZ(:,j,i) = rgb2xyz * sony(:,j,i);
      olympusXYZ(:,j,i) = rgb2xyz * olympus(:,j,i);
   end
end

% The Y value for white should be 100, and all other values
% scaled accordingly

for i=1:18
  sonyXYZ(:,:,i) = sonyXYZ(:,:,i) * 100 / sonyXYZ(2,19,i);
  olympusXYZ(:,:,i) = olympusXYZ(:,:,i) * 100 / olympusXYZ(2,19,i);
end
     
% Convert the XYZ values to LAB

for i=1:18
     sonyLAB(:,:,i) = xyz2lab(sonyXYZ(:,:,i), sonyXYZ(:,19,i));
     olympusLAB(:,:,i) = xyz2lab(sonyXYZ(:,:,i), sonyXYZ(:,19,i));
end

% Calculate the deltaE values

sonydiff = macbethLAB_filters - sonyLAB;
olympusdiff = macbethLAB_filters - olympusLAB;

for i=1:24
  for j=1:18
    sonyDeltaEab(i,j) = norm(sonydiff(:,i,j));
    olympusDeltaEab(i,j) = norm(olympusdiff(:,i,j));
  end
end
