% This script analyzes the data from the gray-world assumption
% experiment.

load /local/class/psych221/tutorials/cMatch/monitor;
load XYZ
rgb2xyz = XYZ'*phosphors;
xyz2rgb = inv(rgb2xyz);

load graySpecData

% Examine the four corner gray control patches to detect color
% shifts.

lightGrayIndex = [1 24];
darkGrayIndex = [4 21];
otherColorsIndex = [2 3 5:20 22 23];

controlLightGrays = graySpecData(lightGrayIndex, :, 1);
controlDarkGrays = graySpecData(darkGrayIndex, :, 1);
controlColors = graySpecData(otherColorsIndex,:,1);

controlLightGrays = sum(controlLightGrays)/2;
controlDarkGrays = sum(controlDarkGrays)/2;

% Analyze only relative differences. The automatic exposure
% mechanism of the camera makes it impossible to take images
% under the same conditions each time.

% Yellow-shift
% ------------
% The test-pattern's R and G components were scaled by a
% factor of 1.3. The gray patches remain unchanged.

yellowLightGrays = graySpecData(lightGrayIndex, :, 2);
yellowDarkGrays = graySpecData(darkGrayIndex, :, 2);
yellowControlColors = graySpecData(otherColorsIndex,:,2);

yellowLightGrays = sum(yellowLightGrays)/2;
yellowDarkGrays = sum(yellowDarkGrays)/2;

% Now compute scale factors.

lightscale = yellowLightGrays./controlLightGrays;
darkscale = yellowDarkGrays./controlDarkGrays;
colorsscale = yellowControlColors./controlColors;

% Dividing throughout by the blue scale factor

yellowlightscale = lightscale./lightscale(3)
yellowdarkscale = darkscale./darkscale(3)
yellowlightgeometricmean = sqrt(yellowlightscale(2).*yellowlightscale(1))
yellowdarkgeometricmean = sqrt(yellowdarkscale(2).*yellowdarkscale(1))
darkreciprocal = 1./yellowdarkgeometricmean
lightreciprocal = 1./yellowlightgeometricmean

colorsscale = colorsscale./(colorsscale(:,3)*[1 1 1]);
figure
plot(colorsscale(:,1),'r')
hold on
axis([1 20 0 2])
grid on
plot(colorsscale(:,2),'g')
h = line([1 20],[1.3 1.3])
set(h,'Color',[.4 .4 .4])
text(2.1,1.35,'1.3')
xlabel('Sample index')
ylabel('Scale factor')
title('Scale factor variation (Yellow)')
hold off

% Blue-shift
% ------------
% Repeat the analysis for the test pattern where only the
% blue component is increased by a factor of 1.3.

blueLightGrays = graySpecData(lightGrayIndex, :, 3);
blueDarkGrays = graySpecData(darkGrayIndex, :, 3);
blueControlColors = graySpecData(otherColorsIndex,:,3);

blueLightGrays = sum(blueLightGrays)/2;
blueDarkGrays = sum(blueDarkGrays)/2;

% Now compute scale factors.

lightscale = blueLightGrays./controlLightGrays;
darkscale = blueDarkGrays./controlDarkGrays;
colorsscale = blueControlColors./controlColors;

% Dividing throughout by the red scale factor

bluelightscale1 = lightscale./lightscale(1)
bluedarkscale1 = darkscale./darkscale(1)
darkreciprocal = 1./bluedarkscale1(3)
lightreciprocal = 1./bluelightscale1(3)

colorsscale = colorsscale./(colorsscale(:,1)*[1 1 1]);
figure
plot(colorsscale(:,3),'b')
hold on
axis([1 20 0 2])
grid on
plot(colorsscale(:,2),'g')
h = line([1 20],[1.3 1.3])
set(h,'Color',[.4 .4 .4])
text(2.1,1.35,'1.3')
xlabel('Sample index')
ylabel('Scale factor')
title('Scale factor variation (Blue, holding Red constant)')
hold off


% Dividing throughout by the green scale factor

bluelightscale2 = lightscale./lightscale(2)
bluedarkscale2 = darkscale./darkscale(2)
darkreciprocal2 = 1./bluedarkscale2(3)
lightreciprocal2 = 1./bluelightscale2(3)

colorsscale = colorsscale./(colorsscale(:,2)*[1 1 1]);
figure
plot(colorsscale(:,1),'r')
hold on
axis([1 20 0 2])
grid on
plot(colorsscale(:,3),'b')
h = line([1 20],[1.3 1.3])
set(h,'Color',[.4 .4 .4])
text(2.1,1.35,'1.3')
xlabel('Sample index')
ylabel('Scale factor')
title('Scale factor variation (Blue, holding Green constant)')
hold off







