close all;
clear all;
cd /local/class/psych221/tutorials
load macbethChart; 
greenChip = 11; redChip = 15; whiteChip = 19; grayChip = 20;
macbethImage = flipud(rot90(reshape(1:24,6,4)));

load D65; 
lgt = D65;

load cMatch/monitor; 			% This contains the monitor phosphors
load cMatch/monitorGam;  % This contains monitorInvGamma
load XYZ;
cd ~/proj;

% Factor for different colored casts on images
factor=1.5;

% Monitor gamma value
myGamma=2.7;

macbethXYZ = XYZ'* diag(lgt)* macbethChart;
xyz2rgb = inv(XYZ'*phosphors);
macbethLinearRGB = xyz2rgb*macbethXYZ;

% Darken surface 20
macbethLinearRGB(:,20)=0.4*macbethLinearRGB(:,20);

wht = macbethLinearRGB(:,whiteChip);
macbethLinearRGB = diag(1./wht)*macbethLinearRGB;

% Black out surfaces 19,16
macbethLinearRGB(:,16)=0;
macbethLinearRGB(:,19)=0;

macbethLinearRGB = scale(macbethLinearRGB,0,1);

% Perform Gray World assumption color balancing for display purposes.
grayLinearRGB = GWbalance(macbethLinearRGB');
grayLinearRGB = scale(grayLinearRGB,0,1);

grayRGB = (grayLinearRGB).^(1/myGamma);
figure; colormap(grayRGB); 
image(macbethImage);
title('Forced Gray World Assumption Image')



% Perform Gray World assumption color balancing on generated image.
newLinearRGB = GWbalance(macbethLinearRGB');

% Add in color filters here to grayworld image.
newLinearRGB(:,1) = factor*newLinearRGB(:,1);
newLinearRGB(:,2) = factor*newLinearRGB(:,2);
%newLinearRGB(:,3) = factor*newLinearRGB(:,3);

newLinearRGB = scale(newLinearRGB,0,1);
newRGB = (newLinearRGB).^(1/myGamma);
newColorMap = newRGB;

figure; colormap(newColorMap)
image(macbethImage);
title('Filtered Forced Gray World Assumption Image');


% Perform Gray World assumption color balancing on filtered Image.
testLinearRGB = GWbalance(newLinearRGB);
testLinearRGB = scale(testLinearRGB,0,1);
testRGB = (testLinearRGB).^(1/myGamma);
testColorMap = testRGB;

figure; colormap(testColorMap)
image(macbethImage);
title('Gray World on Filtered Image');

