%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  File: ColorTutorial.m
%%%  Original Lisp Code written by E.J. Chichilnisky, summer 1992 
%%%  Converted to Matlab by Geoff Boynton, summer 1996
%%%  Changed and Updates by hagit Hel-Or, winter 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% In this tutorial, we will introduce some basic spectral examples and calculations,
%% Color Matching and the CIEXYZ color system.

%% You will notice that much of the work below involves matrix algebra.  I
%% suggest that you pull out a pad of paper and draw matrix tableaus as we
%% proceed to help you follow the calculations.

%% We suggest you set your 'monitors' control panel  to show 
%% 'thousands' of colors.  Otherwise you can't view more than
%% one image at a time.

%% This tutorial uses m and mat files i directory color
path('color',path)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Examples of Spectral Distribution Functions (SPD) or spectra

%% The Macbeth Color Checker is a set of 24 surfaces commonly used to
%% evaluate color balancing systems.  The materials in the Color
%% Checker were selected to have the same surface reflectance
%% functions as various important surfaces that are commonly used in
%% television and film images.  We will read in a matrix of data.
%% Each row of the matrix is the reflectance function of one of the
%% Macbeth surfaces, measured at each of 31 samples in the visible
%% spectrum.  (Most of the visible spectrum is in the 400-700
%% nanometers wavelength region.)  Thus we get a 24x31 matrix of
%% surface spectra. Each row is one surface.

%% we'll define the spectrum as follows:
spectrum=linspace(400,700,31);

%% this loads in the macbeth color checker and the munsell 
%% standard surfaces

load surfaces  %macbeth, munsell

% Display an image of the Macbeth chart.
% WARNING : this is an image (from the internet) and is similar but not exactly
%  what the Macbeth spectra would look like on the display.
%  This can be calculated and will be done in a later script.

im = imread('macbeth.tif');
figure(1); 
image(im);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOOK AT THE SURFACE REFLECTANCE SPECTRA

%% The 8th row of this matrix is a surface that typically looks
%% greenish.  We plot the fractional reflectance (as a function of
%% wavelength) for this surface.

figure(1);clf
plot(spectrum,macbeth(8,:),'g');
xlabel('Wavelength (nm)');
ylabel('Surface Reflectance');
title('Reflectance of Macbeth Surfaces');

%% For example, about 28% of the light at wavelength 500 nm is
%% reflected by the "Green" surface.

green_500=macbeth(8,find(spectrum==500))

%% Here is a red surface:

hold on
plot(spectrum,macbeth(9,:),'r');

%% And here is a gray surface:

plot(spectrum,macbeth(20,:),'k');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOOK AT TYPICAL DAYLIGHT SPECTRA.

%% Load in the illuminants.
%% Each illuminants' entries represent the amount of light present
%% at each wavelength.  Again, the light is sampled at 31 points in the
%% visible spectrum, so we get a 31-vector.

load illuminants %cie_a,daylight_65,flourescent,illuminant_C,xenon_flash

%% Make a plot of daylight_65. Note: the energy units are arbitrary.

figure(1);clf
plot(spectrum,daylight_65,'y');
xlabel('Wavelength (nm)');
ylabel('Energy');
title('Energy spectrum of daylight 65 Illimunant');

%% daylight_65 is the standard illuminant which represents a mix of blue 
%% sky and clouds. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DIFFERENT LIGHT SOURCES

%% plot the energy spectrum of the xenon_flash illuminant in
%% green, and daylight_65 in yellow.

figure(1);clf
plot(spectrum,xenon_flash,'g');
hold on
plot(spectrum,daylight_65,'y');
hold off
xlabel('Wavelength (nm)');
ylabel('Energy');
title('Energy Spectrum of Illuminants');
legend('xenon_flash','daylight 65');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% READ IN THE HUMAN CONE SPECTRAL SENSITIVITIES.

%% Now we want to measure the photoresponses of the human cones to
%% these stimuli.  The human cone spectral sensitivities have been
%% measured and correspond closely with behavioral color matching data
%% known for many years.  We now read in a matrix of the human cone
%% sensitivities.  There are three classes of cone, the so-called L,M,
%% and S cones ((L)ong-, (M)iddle-, and (S)hort-wavelength peak
%% sensitivities). Again, we represent each cone by its sensitivity at
%% each of 31 sample wavelengths in the spectrum.  Thus, we have a
%% 3x31 matrix representing the human cone spectral sensitivities.

load cones  %cones

%% Look at the three cone classes superimposed.  Notice how close the
%% L and M cones are in terms of their peak sensitivities.

figure(1)
clf
plot(spectrum,cones(1,:),'r');
hold on
plot(spectrum,cones(2,:),'g');
plot(spectrum,cones(3,:),'b');
xlabel('Wavelength (nm)');
ylabel('Relative Sensitivity');
title('Cone Spectral Sensitivity Functions');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CALCULATE THE HUMAN CONE RESPONSES TO THE MACBETH SURFACES
%% AND TO THE ILLUMINANTS.

%% We can describe the receptor encoding as a matrix multiplication
%% too. The matrix product of the receptor sensitivites and the
%% transposed spectral signals gives the photoisomerizations in
%% each receptor class due to each of the surfaces. Hence, we have a
%% 3x24 matrix of receptor responses.

cone_signals_surfaces=cones*macbeth';
cone_signals_illuminant=cones*daylight_65'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Color matching 
%% Load CIE-RGB Color Matching Functions
load ciergb10

%% Calculate the Tristimulus values for the Macbeth surfaces and 
%% for the illuminant
tristim_RGB_surfaces=ciergb10*macbeth';
tristim_RGB_illuminant=ciergb10*daylight_65'



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CIEXYZ
%% Load and display the CIEXYZ Color Matching Functions.
%%

load xyz

%% Calculate the Tristimulus values for the Macbeth surfaces and 
%% for the illuminant
tristim_XYZ_surfaces=ciexyz*macbeth';
tristim_XYZ_illuminant=ciexyz*daylight_65'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Metamers
%% Cones, CIERGB and CIEXYZ preserve metamers

%%
%% Load two Metameric Spectras
load metamers

figure(1)
clf
plot(spectrum,metamers(1,:),'r');
hold on
plot(spectrum,metamers(2,:),'g');
xlabel('Wavelength (nm)');
ylabel('Energy');
title('Spectral Distribution Functions of two metamers');

%% Compare Tristimulus values
cones * metamers'
ciergb10 * metamers'
ciexyz * metamers'


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Find the linear transformation (3x3 matrix) that maps:
%% the cone spectral sensitivity to the cieRGB and cieXYZ CMF
%% and the two CMF. 
%% That is find the 3x3 matrix cmf12cmf2 such that:
%%             triStimulusVals2 = cmf12cmf2 * triStimulusVals1
%% To do this we use the Pseudo Inverse function of Matlab.
%% This will give the best Least Squares Estimate of the linear transformation.



RGB2lms =  cones * pinv(ciergb10);
lms2XYZ = ciexyz * pinv(cones);
XYZ2RGB = ciergb10 * pinv(ciexyz);


%% These linear transformations should also map the tristimulus values of the corresponding
%% Color Space. We test this on the illuminant:

[cone_signals_illuminant   RGB2lms*tristim_RGB_illuminant]
[tristim_XYZ_illuminant    lms2XYZ*cone_signals_illuminant]
[tristim_RGB_illuminant    XYZ2RGB*tristim_XYZ_illuminant]


%% We can measure the error over many such measurements
sqrt(sum(sum(abs(lms2XYZ*cone_signals_surfaces - tristim_XYZ_surfaces))))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CIE Diagram 
%%
%% Draw the CIE Diagram boundary by calculating and plotting
%% the xy coordinates of the monochromatic spectra. 
%% xy are calculated from XYZ of the monochromatic lights.
%% The XYZ of the monochromatic lights are in the columns of the XYZ cmf.
%%

s = sum(ciexyz);
xyCoor = ciexyz([1,2],:) ./ (ones(2,1)*s);

figure(1); clf;
plot(xyCoor(1,:),xyCoor(2,:),'k');
%plot the nonspectral locus line
hold on
line([xyCoor(1,1),xyCoor(1,end)],[xyCoor(2,1),xyCoor(2,end)],'Color',[0 0 0]);
tick = [0:.2:.8];
axis equal
set(gca,'xlim',[0 .8],'ylim',[0 .85],'xtick',tick,'ytick',tick)
grid on


%% Given a spectra, plot its corresponding point in the CIE diagram.
%% Plot the points for the illuminants:
%%  cie_a,daylight_65,floursecent,illuminant_C,xenon_flash

illuminants = [cie_a ; daylight_65 ; flourescent; illuminant_C ; xenon_flash];
XYZCoor = ciexyz * illuminants';
xyCoor = XYZCoor([1,2],:) ./ (ones(2,1)*sum(XYZCoor));

plot(xyCoor(1,:),xyCoor(2,:),'ro');


% Notice the chromaticity coordinates of the illuminants fall on
%% the balck body radiator line.


%%% ADD Conversions between spaces - later or next tutorial