% This script analyzes the data from the perfect-reflector
% experiment.

load graySpecData

indices = [1:9 12 13 16:24];

controlColors = graySpecData(indices, :, 4);
yellowColors = graySpecData(indices, :, 5);
blueColors = graySpecData(indices, :, 6);

figure
image(reshape(1:20,4,5));
colormap(controlColors);

% We now correct the yellow and blue images for exposure
% differences. A simple way of doing that is to make sure
% that (R+G+B) in all three images are the same.
% Another way of doing this is to equate the average Y
% values -- I avoided this approach because it scales
% each color channel differently.

average = sum(sum(controlColors));

yellowRGB = yellowColors.*average/sum(sum(yellowColors));
blueRGB = blueColors.*average/sum(sum(blueColors));

figure
image(reshape(1:20,4,5))
colormap(yellowRGB);

figure
image(reshape(1:20,4,5))
colormap(blueRGB);

yellowscalefactors = yellowRGB./controlColors;
bluescalefactors = blueRGB ./controlColors;

figure
plot(yellowscalefactors(:,1),'r')
hold on
axis([1 20 0 2])
grid on
plot(yellowscalefactors(:,2),'g')
plot(yellowscalefactors(:,3),'b')

plot(bluescalefactors(:,1),'r--')
plot(bluescalefactors(:,2),'g--')
plot(bluescalefactors(:,3),'b--')
hold off

title('Scale factors for R, G and B')
xlabel('Sample index')
ylabel('scale')

legend('Yellow highlight R', ...
    'Yellow highlight G',...
    'Yellow highlight B',...
    'Blue highlight R',...
    'Blue highlight G',...
    'Blue highlight B')

overallyellowscale = sum(yellowRGB)./sum(controlColors)
overallbluescale = sum(blueRGB)./sum(controlColors)








