%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  File: T_IST05_Color2_Illumination.m
%%%  Source from Wandell Labs - Stanford
%%%  Changed and Updated by hagit Hel-Or, winter 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In this Tutorial we will study the Image Formation model, 
% Illuminant and surface linear models, surface and illuminant estimation,
% white balance and Illumination Correction.

%% 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)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

spectrum=linspace(400,700,31);

% load illuminants and surfaces
%D65 = daylight, color temperature 6500K. 
%CIE Illuminant_A = tungsten, color temperature 2856K
%CIE Illuminant_C = filtered CIE Illuminant A operating at 6700K
%fluorescent = cool white. 
% Add pure white illumination

load illuminants %cie_a,daylight_65,flourescent,illuminant_C,xenon_flash
illuminants = [cie_a./2 ; daylight_65 ; flourescent; illuminant_C ; 90*ones(1,length(cie_a))];
illuminantNames = {};
[illuminantNames{1:5,1}] = deal('cie a' , 'daylight 65' , 'flourescent' , 'cie C' , 'white')

load surfaces  %macbeth, munsell


% Image Formation Model
% Sensor responses are equal to Illuminant times surface times sensors
% sensors in this case will be the Human Cones

load cones  %cones

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DIFFERENT SURFACES UNDER DIFFERENT LIGHT SOURCES

% Calculate and display the cone responses for different surfaces under 
% different illuminants.


selectedSurfaces = [1 4 5 14 24]; %surfaces to plot from Munsell

figure;
subplot(5,5,1);
for j = 1:5
   coneResponses=cones*diag(illuminants(j,:))*macbeth';
   for i = 1:5
      subplot(5,5,(j-1)*5+i);
      bar(coneResponses(:,selectedSurfaces(i)));
   end;
   subplot(5,5,(j-1)*5+1);
   ylabel(illuminantNames(j));
end;
% write the titles
for i = 1:5
   subplot(5,5,20+i);
   xlabel(['Surface ' num2str(selectedSurfaces (i))]);
end;
subplot(5,5,3);
title('Cone Sensor Responses to Surfaces under Different Illuminants','Fontsize',18);


% plot color patch that approximates the visual perception of  these 
% cone responses. Use function cones2monitor. This function assumes
% standard phosphors, so color patch might not look exactly as it should.
% To get correct color patch, monitor calibration must be performed.

% create colormap with all colors of the patches.
% Use the function cone2monitor. This function will be discussed in a
% seperate tutorial.
cmap = [];
for j = 1:5
   coneResponses=cones*diag(illuminants(j,:))*macbeth';
   monitorRGB = cones2monitor(coneResponses(:,selectedSurfaces));
   cmap= [cmap; monitorRGB'];
end;
figure;
image(reshape(1:25,5,5)');
% colormap must be between 0 and 1
if min(cmap(:))<0 
   cmap = cmap-min(cmap(:))';
end;
cmap = cmap./max(cmap(:));
% normalize so that all have same intensity
%cmap = cmap./(sum(cmap')'*ones(1,3));
colormap(cmap);
set(gca,'YTick',[],'XTick',[]);
title('Surfaces under Different Illuminants','Fontsize',18);
xlabel('Surfaces','Fontsize',18);
ylabel('Illuminants','Fontsize',18);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SHOW COLOR CONSTANCY WORKING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create scene with large bg and several patches.
% 

% image under illuminant1
coneResponses=cones*diag(illuminants(5,:))*macbeth';

cmap = [];
cmap = [cmap; cones2monitor(coneResponses(:,23))'] ; % bg
cmap = [cmap; cones2monitor(coneResponses(:,1))'];   % patch1
cmap = [cmap; cones2monitor(coneResponses(:,5))'];   % patch2
cmap = [cmap; cones2monitor(coneResponses(:,11))'];  % patch3
cmap = [cmap; cones2monitor(coneResponses(:,14))'];  % patch3

im = ones(20,20);
im(8,8) = 2;
im(8,12) = 3;
im(12,8) = 4;
im(12,12) = 5;


% image under illuminant2
coneResponses=cones*diag(illuminants(2,:))*macbeth';
cmap = [cmap; cones2monitor(coneResponses(:,23))'] ;% bg
cmap = [cmap; cones2monitor(coneResponses(:,1))'] ;  % patch1
cmap = [cmap; cones2monitor(coneResponses(:,5))'];   % patch2
cmap = [cmap; cones2monitor(coneResponses(:,11))'] ; % patch3
cmap = [cmap; cones2monitor(coneResponses(:,14))'];  % patch3

im2 = ones(20,20).*6;
im2(8,8) = 7;
im2(8,12) = 8;
im2(12,8) = 9;
im2(12,12) = 10;

figure;
subplot(1,2,1)
image(im);
axis off;
subplot(1,2,2)
image(im2);
cmap = cmap./max(cmap(:));
colormap(cmap)
axis off;


% STRETCH FIGURE TO FILL SCREEN AND FOCUS ON THE 
% MID POINT BETWEEN THE TWO FIGURES.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LINEAR MODELS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LINEAR MODEL FOR SURFACES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% We constructed this display by calculating the full spectral signal
%% from the surface reflectances and light spectrum, then calculating
%% the receptor responses and recreating those responses on the
%% monitor.  Unfortunately, 31 spectral samples are a rather expensive
%% way to represent the spectral properties of surfaces.  How can we
%% compress the information in the reflectances and still achieve an
%% accurate rendition of the surfaces?  A common method is to use the
%% principal components of the surface set.  Below, we use the
%% singular value decomposition (SVD) to reconstruct a
%% lower-dimensional representation of the full surface set. The variable
%% 'num_bases' determines the dimensionality of the reconctructed
%% macbeth surfaces.


num_bases=4;

[u,s,v]=svd(macbeth'*macbeth);

% the Basis vectors are first num_components columns in u
surfaceBasis = u;

% plot the Basis SPDs

figure; 
plot(u(:,1:num_bases));
legend('s1','s2','s3','s4');


% calculate the coefficients of the basis vectors for each surface
% by projecting the surface spd onto the basis.

sigmaSurfaces = macbeth*u(:,1:num_bases);

% Calculate the approximated SPD by multiplying the coefficients with
% the basis vectors.

approxMacbeth = sigmaSurfaces * u(:,1:num_bases)';


% compare
figure;
for i = 1:size(macbeth,1)
   subplot(4,6,i);
   plot(spectrum,macbeth(i,:),'r');
   hold on;
   plot(spectrum,approxMacbeth(i,:),'g');
end;

subplot(4,6,3);
title(['    Macbeth Surfaces Approximated by ' int2str(num_bases) ' basis SPDs'],'Fontsize',18);

   
% NOW GO BACK AND RERUN THE CODE WITH A DIFFERENT num_bases.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LINEAR MODEL FOR Illuminants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the same as above for a collection of illuminants.
% Instead of repeating above, here is a trick:

num_bases = 3;

[u,s,v]=svd(illuminants');

% the Basis vectors are first num_components columns in u
illuminantBasis = u;

% plot basis SPDs
figure;
plot(illuminantBasis(:,1:num_bases));
legend('s1','s2','s3');

% zero out all but the first num_components eigenvalues in s
s(num_bases+1:end,:) = 0;

% reconstruct approximated surfaces by multiplying SVD components.
approxIlluminants = (u*s*v')';

% compare
figure;
for i = 1:size(illuminants,1)
   subplot(1,5,i);
   plot(spectrum,illuminants(i,:),'r');
   hold on;
   plot(spectrum,approxIlluminants(i,:),'g');
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ILLUMINATION CORRECTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Change image acquired under one illumination, to
% appear as if taken under a different illumination.

im1 = imread('color/tungsten.tif');
im2 = imread('color/sun.tif');

f1 = figure;
image(im1)
axis off

f2 = figure;
image(im2)
axis off

% change type of variable for computation
im1 = double(im1);
im2 = double(im2);


% Correct im1 to look as if illuminated as im2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I.
% Correct using Von Cries Approach. Thus a single surface viewed in both 
% images is required.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% corner of teacher's table was chosen
x1 = 303; y1 = 706;
x2 = 310; y2 = 726; 

% or choose your own matching surface points: 
% figure(f1); [x1,y1] = ginput(1);   
% figure(f2); [x2,y2] = ginput(1);  


% Look at the two RGB values. They differ due to illumination
RGB1 = squeeze(im1(round(y1),round(x1),:))'
RGB2 = squeeze(im2(round(y2),round(x2),:))'


im1new = im1;
im1new(:,:,1) = im1new(:,:,1) * RGB2(1)./RGB1(1);
im1new(:,:,2) = im1new(:,:,2) * RGB2(2)./RGB1(2);
im1new(:,:,3) = im1new(:,:,3) * RGB2(3)./RGB1(3);

% display new image
figure;
image(uint8(im1new));
axis off


% Try the above code using a surface which is NOT achromatic:
x1 = 796; y1 = 288;
x2 = 812; y2 = 301; 
% OR
x1 = 196; y1 = 481;
x2 = 185; y2 = 497; 

% Notice that correction based on chromatic surfaces is not so good.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% II.
% Can we do better with chromatic surfaces only? Find more than
%  1 chromatic surface in both images and find the best transformation matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% choose 5 chromatic surfaces
x1 = [206   446   793   838   727]';
y1 = [479   324   292   270   533]';
x2 = [185   446   798   831   725]';
y2 = [501   346   304   283   551]';

% or choose your own matching surface points: 
% figure(f1); [x1,y1] = ginput(5);   
% figure(f2); [x2,y2] = ginput(5);  

% show the points in the images
figure(f1);
hold on;
plot(x1,y1,'m*')

figure(f2);
hold on;
plot(x2,y2,'c*')


% Create matrices of RGB values of the two images
RGB1 = [diag(im1(round(y1),round(x1),1)) diag(im1(round(y1),round(x1),2)) ...
           diag(im1(round(y1),round(x1),3))]';
RGB2 = [diag(im2(round(y2),round(x2),1)) diag(im2(round(y2),round(x2),2)) ...
           diag(im2(round(y2),round(x2),3))]';

% Find matrix M that transforms im1 values to those illuminated as in im2:
% im2 = M * im1;
% Use psuedo-inverse:  M = im2 * pinv(im1);
M = RGB2 * pinv(RGB1);

% Apply M to all pixels in im1:

% create matrix with all RGB values of im1
im1RGB = reshape(im1(:),length(im1(:))/3,3);

% transform im1 colors and reshape 
im1RGBnew = (M * im1RGB')';
im1new = reshape(im1RGBnew(:),size(im1,1),size(im1,2),3);

% make sure values are not out of range
im1new(find(im1new>255)) = 255;
im1new(find(im1new<0)) = 0;


figure; 
image(uint8(round(im1new)))
axis off




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% III.
% Gray World Assumption - mean of all surfaces are gray.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Look at the two mean RGB values. They differ due to illumination
RGB1 = mean(reshape(im1(:),length(im1(:))/3,3))
RGB2 = mean(reshape(im2(:),length(im2(:))/3,3))

% Change im1 so that its mean equals that of im2.
im1new = im1;
im1new(:,:,1) = im1new(:,:,1) * RGB2(1)./RGB1(1);
im1new(:,:,2) = im1new(:,:,2) * RGB2(2)./RGB1(2);
im1new(:,:,3) = im1new(:,:,3) * RGB2(3)./RGB1(3);

% display new image
figure;
image(uint8(im1new));
axis off





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXTRA - delete when script done

% show macbeth patches under white illuminant
coneResponses=cones*diag(illuminants(5,:))*macbeth';
cmap = [];
for j = 1:24
   monitorRGB = cones2monitor(coneResponses(:,j));
   cmap= [cmap; monitorRGB'];
end;
figure;
image(reshape(1:24,6,4)');
cmap(find(cmap<0)) = 0;
%cmap = cmap - min(cmap(:));
cmap = cmap./max(cmap(:));
colormap(cmap);
axis off;
title('Macbeth Surfaces','Fontsize',18);