% calibrate.m
%
% Depository for random code, until I figure out how to partition
% it better.
%

% These values depend on the set of images
% TakeTwo/cam1
saturationValue = 180;
minValue = 10;
macBethSpectrum='TakeTwo/cam1/t2mbSpecIndiv.mat';
imageDir='TakeTwo/cam1';

% Load macBeth color checker photospectrometer data
% This loads mbSpec, the intensity per wavelength for each
% checker, and mbWave, the wavelengths at which the intensities
% were measured. The checkers are ordered by row, starting with 
% the grayscale row first. Here are a few for reference:
% Checker 1: white 6: black
% checker 7: blue  8: green 9: red

load(macBethSpectrum);
load XYZ;

% This XYZ has values for 370:730
% Spectrometer output is 380:780
% Grab just the portion of XYZ that we need (380-728)
% and portion of spectrometer output to match XYZ input range

neededIndices = [11:4:359];
cropXYZ = XYZ(neededIndices,:);
macSpec = mbSpec(1:88,:);
macWave = mbWave(1:88);

% Determine XYZ values of color checkers,
% then scale Y values so white square has Y of 1.0
macXYZ = cropXYZ' * macSpec;
macXYZ = macXYZ / macXYZ(2,1);

% plot(mbWave,mbSpec(:,1))  % plot the white spectra -- reddish light!
% plot(mbWave,mbSpec(:,7))  % plot the blue spectra

% Test this against the photospectrometer output (correct).
% macChroma = chromaticity(macXYZ);

% Load the images.
% Grabbing the color patches out of each image set is a bit
% different, so I've put a function called "grabEachPatch" in the
% directories of each set of images. It averages the dark images,
% subtracts that value from the macBeth images, 

cd(imageDir);
measRGB=grabEachPatch;
cd ../..

% Doh! This exercise shows two things--first, that the sensor is
% saturating, and second, it's response is linear until it saturates.
% Redo measurements with non-saturating values.
% This is a function of the lighting, too--the red sensor was
% pegging first on the grayscale color checkers.

figure;
hold on
plot(macXYZ(2,1:6),measRGB(1,1:6),'ro-');
plot(macXYZ(2,1:6),measRGB(2,1:6),'go-');
plot(macXYZ(2,1:6),measRGB(3,1:6),'bo-');
title('Measured R, G and B vs. measured Y for gray scale checkers');
xlabel('Luminance of checker;);
ylabel('Pixel intensity');
hold off  

% Pick out the color checkers that aren't saturated or too low.
maxVals = max(measRGB);
minVals = min(measRGB);
validCheckers = find(maxVals < saturationValue & minVals > minValue)
validXYZ = macXYZ(:,validCheckers)
validRGB = measRGB(:,validCheckers)


%%% MODEL FITTING
%%% Computing errors to data used for fit.

% Ok, now it's finally time to do some model fitting. Start with
% a 3x3 transformation from RGB to XYZ. See how good the fit
% is. Does it matter which direction we go with these
% transformations? I'm more interested in going from RGB back to
% XYZ.

% First Order Model, 3x3 Transform
% validRGB'*colorTrans' = validXYZ'
colorTrans = validRGB'\validXYZ';
colorTrans = colorTrans'
sensorXYZ=colorTrans*validRGB;
mse1=3*mse(sensorXYZ,validXYZ)
%This is mean squared Euclidean distance in XYZ space.

% Now let's calculate the error is CIE-Lab space.
% Use the measured white point from the photospectrometer.

whitePoint = macXYZ(:,1);
validLab=xyz2lab(validXYZ',whitePoint');
measLab=xyz2lab(sensorXYZ',whitePoint');
mseLab1=3*mse(measLab,validLab)


% Second fitting method
% R G B RG GB BR
[colorTrans2] = Model2(validRGB,validXYZ);
[rgb2] = Model2Rgb(validRGB);
sensorXYZ2=colorTrans2*rgb2;
mse2=3*mse(sensorXYZ2,validXYZ) 
measLab2=xyz2lab(sensorXYZ2',whitePoint');
mseLab2=3*mse(measLab2,validLab)


% Third fitting method
% Use R G B RG GB BR RR GG BB
[colorTrans3] = Model3(validRGB,validXYZ);
[rgb3] = Model3Rgb(validRGB);
sensorXYZ3=colorTrans3*rgb3;
mse3 = 3*mse(sensorXYZ3,validXYZ)
measLab3=xyz2lab(sensorXYZ3',whitePoint');
mseLab3=3*mse(measLab3,validLab)


% So now the multimillion dollar question is,...
% How bad is the error for colors not in the calibration set?
% I fear to calibrate with few patches, so I'm going to leave one
% patch out each time, fit the model, then calculate the error in
% the computed XYZ value for the omitted patch.
% Sum up squared XYZ errors in firstOrderXYZError.
% measRGB is original full set of checkers.
% validCheckers are indices into measRGB.

numIters=size(validCheckers,2);
indices=[1:numIters];
checkerSubset=indices(2:numIters);
computedXYZ1=zeros(size(validXYZ));
computedXYZ2=computedXYZ1;
computedXYZ3=computedXYZ1;
% i = number of the checker that was left out
for i = 1:numIters
  meterXYZ=validXYZ(:,checkerSubset);
  cameraRGB=validRGB(:,checkerSubset);

  % Model 1
  colorTrans = cameraRGB'\meterXYZ';
  colorTrans = colorTrans';
  computedXYZ1(:,i)=colorTrans*validRGB(:,i);

  % Model 2
  cameraRGB2=Model2Rgb(cameraRGB);
  colorTrans2=Model2(cameraRGB,meterXYZ);
  computedXYZ2(:,i)=colorTrans2*rgb2(:,i);
    
  % Model 3
  cameraRGB3=Model3Rgb(cameraRGB);
  colorTrans3=Model3(cameraRGB,meterXYZ);
  computedXYZ3(:,i)=colorTrans3*rgb3(:,i);
  
  if(i<numIters)
    checkerSubset(i)=indices(i);
  end
end

fairMse=zeros(1,3);
fairLabMse=zeros(1,3);

% Compute errors
fairMse(1)=3*mse(computedXYZ1,validXYZ);
computedLab1=xyz2lab(computedXYZ1',whitePoint');
fairMseLab(1)=3*mse(computedLab1,validLab);

fairMse(2)=3*mse(computedXYZ2,validXYZ);
computedLab2=xyz2lab(computedXYZ2',whitePoint');
fairMseLab(2)=3*mse(computedLab2,validLab);

fairMse(3)=3*mse(computedXYZ3,validXYZ);
computedLab3=xyz2lab(computedXYZ3',whitePoint');
fairMseLab(3)=3*mse(computedLab3,validLab);


% Picture Time!
% Let's take a full color checker shot, correct it, and
% display it.
% xyz2lab take Nx3, lab2xyz takes 3xN. Arg.

img=double(imread('TakeTwo/cam1/mb22.bmp'));

correctMB1=TransImg1(img,colorTrans,whitePoint);
figure;imshow(correctMB1);
imwrite(correctMB1,'correctMB1.bmp','BMP');

correctMB2=TransImg2(img,colorTrans2,whitePoint);
figure;imshow(correctMB2);
imwrite(correctMB2,'correctMB2.bmp','BMP');

correctMB3=TransImg3(img,colorTrans3,whitePoint);
figure;imshow(correctMB3);
imwrite(correctMB3,'correctMB3.bmp','BMP');

img=double(imread('TakeTwo/cam1/clutter.bmp'));

correctMB1=TransImg1(img,colorTrans,whitePoint);
figure;imshow(correctMB1);
imwrite(correctMB1,'clutter1.bmp','BMP');

correctMB2=TransImg2(img,colorTrans2,whitePoint);
figure;imshow(correctMB2);
imwrite(correctMB2,'clutter2.bmp','BMP');

correctMB3=TransImg3(img,colorTrans3,whitePoint);
figure;imshow(correctMB3);
imwrite(correctMB3,'clutter3.bmp','BMP');



