% PSYCH221/EE362
% Applied Vision and Image Systems
% Final Project
% March 2006
% Kevin Clinton, Jessica Faruque, Cecylia Wati

% The code is divided into 4 main parts as follows:
% Part 1
% In part 1, the objective is to derive a linear model from MCC (Macbeth
% Color Checker) radiance. The basis functions derived is then used to
% reconstruct the original radiance, and we will see that the basis
% functions do a pretty good  at reconstructing the original data.
% In this part we also check the optimal number of bases to derive the
% linear model (nBases).
%
% Part 2
% In this part we would like to see how well the camera response predicts
% coefficients for MCC radiance. The camera has 9 color channels which are
% attained from condition with no filter, with a yellow filter, and with a
% blue filter.
% 
% Part 3
% In this part, using the basis functions derived from part 1, we would
% like to see how well these basis functions predict skin radiance data.
% The LightAndSkin.mat consists of 8 sets of data, the first two sets are
% the lighting conditions under which we took the measurements (we took two
% measurements for the same lighting condition).
%  
% Part 4
% Lastly we would like to know how well the camera response predicts
% coefficients for skin radiance data. The method is similar to part 2, but
% we're just using different set of radiance data (skin, instead of MCC).

%add ISET dir
cd '~/psych221/ISET-1.0'
isetPath(pwd);

cd '~/psych221/project/MultispectralImaging'
addpath(genpath(pwd));

clear all; close all;
dataDir = ['CreateSceneRadianceImages/Data'];   % Directory where data is stored
wavelength = 400:5:700;  % 61 wavelength samples between 400 and 700 nm

% The following routine will get the reflectance of MCC.
allSurfaces = {'macbethChart.mat'};
surfaces = [];
for ii=1:length(allSurfaces)
    fullName = fullfile(dataDir,'Surfaces',allSurfaces{ii}); % location of surface reflectance data
    theseSurfaces = vcReadSpectra(fullName,wavelength);
    surfaces = [surfaces, theseSurfaces];
end
% plot(wavelength,surfaces);


% Create low-dimensional linear model for surface radiance data
%
%The camera captures the radiance of a scene, not the reflectance. 
% Since we would like to be able to compare our data with that from the
% camera, we need to find the radiance instead of reflectance. 
% We can multiply the reflectance with the light to get the radiance
% data (radiance=reflectance*light).

fullName = fullfile(dataDir,'Illuminants','LightNoFilter'); % location of measured light data
SceneIllum = vcReadSpectra(fullName,wavelength);
radiance=diag(SceneIllum)*surfaces;
%figure;plot(radiance);

nBases = 6;
if ~exist('nBases','var'),  nBases = []; end
[sBasis,sValues] = lmComputeBases(radiance,0,nBases);
% figure; plot(sBasis);
coef = sBasis'*radiance; % [6xN] where N is the number of surfaces
% Computing the mse between the linear model and the surface reflectances
% plot(sBasis);
for kk=1:nBases
    dd = radiance - sBasis(:,1:kk)*coef(1:kk,:);
    dd = sum(dd.^2);
    serr(kk) = mean(dd);
end
%figure; semilogy(serr/sum(serr),'*-');

% Plot for part 1
% How well do N basis functions predict MCC
predictedMCC=sBasis*coef;
figure;
for i=1:size(radiance,2)
    plot(wavelength,radiance(:,i),'b');
    hold on; plot(wavelength, predictedMCC(:,i),'r');
    legend('MCC Radiance','Predicted MCC (From Linear Model)');
end
grid on;
title('Part 1: How well does linear model predict coefficients for MCC');
%%%%%%%%%%%%%%%%%%%%%End of part 1%%%%%%%%%%%%%%%%%%%%%%%

% Read in the D70 sensor spectral response 
fullName = fullfile(dataDir,'Sensors','NikonD70.mat'); % location of Nikon sensor data
D70 = vcReadSpectra(fullName,wavelength);
%figure; plot(D70);

% Read in the measured colored lights.  
allLights = {'LightNoFilter','LightFilter809','LightFilter849'}; % place holder for real data
lights = [];
for ii=1:length(allLights)
    fullName = fullfile(dataDir,'Illuminants',allLights{ii}); % location of measured light data
    theseLights = vcReadSpectra(fullName,wavelength);
    lights = [lights, theseLights];
end

% Now, we estimate the filter transmittance.  The first light has no
% filter, so we divide the measurements by the first light to get the
% transmissivities.
% lights is 61x12 matrix, we took 12 kinds of lights
% diag(1./lights(:,1)) will create diagonal matrix of size 61x61
% the first command of filters will create 61x12 normalized matrix
filters = diag(1./lights(:,1))*lights;
filters = filters(:,2:end);  % These are the filter transmissivities
nFilters = size(filters,2);
% figure; plot(wavelength,filters)

noiseSD = 0.00001;
ListOfFilters = filters(:,[1,2]);
% We combine the filters to get 9 color channels
sensor = mcCombineSensorsFilters(D70,ListOfFilters); 

% Camera response to radiance
cameraR   = sensor'*radiance;  
cameraR   = cameraR + noiseSD*randn(size(cameraR)); % 9xN - 9 color channels and N surfaces

%cameraR is approximately equal to sensor'*sBasis*coef;
estCoef=pinv(sensor'*sBasis)*cameraR;  % [6xN] where N is the number of color signals

%figure; plot(coef,estCoef);
cameraPredictedMCC = sBasis *estCoef ;

% Plot for part 2
% How well does camera response predict coefficients for MCC
figure;
for i=1:size(radiance,2)
    plot(wavelength,radiance(:,i),'b');
    hold on; plot(wavelength, cameraPredictedMCC(:,i),'r');
    legend('MCC Radiance','Predicted MCC (From Camera)');
end
grid on;
title('Part 2:How well does camera response predict coefficients for MCC');
%%%%%%%%%%%%%%%%%%%%%End of part 2%%%%%%%%%%%%%%%%%%%%%%%

 SkinFile = 'LightAndSkin.mat'
 Skinfullname = fullfile(dataDir,'Surfaces',SkinFile);
 SkinRadiance = vcReadSpectra(Skinfullname,wavelength);
 SkinRadiance = SkinRadiance(:,3:8);
 
 SkinCoef=sBasis'*SkinRadiance;
 predictedSkinRadiance=sBasis*SkinCoef;
 
% Plot for part 3
% How well does linear model for Macbeth radiance predict skin radiance
% data
figure;
for i=1:size(SkinRadiance,2)
    plot(wavelength,SkinRadiance(:,i),'b');
    hold on; plot(wavelength, predictedSkinRadiance(:,i),'r');
    legend('Skin Radiance','Predicted Skin Radiance (From Linear Model)');
end
grid on;
title('Part 3:How well does linear model predict coefficients for skin radiance data');
%%%%%%%%%%%%%%%%%%%%%End of part 3%%%%%%%%%%%%%%%%%%%%%%%

% Camera response to radiance
cameraR   = sensor'*SkinRadiance;  
cameraR   = cameraR + noiseSD*randn(size(cameraR)); % 9xN - 9 color channels and N surfaces

%cameraR approx= sensor'*sBasis*coef;
estSkinCoef=pinv(sensor'*sBasis)*cameraR;  % [6xN] where N is the number of color signals
   
cameraPredictedSkinRadiance=sBasis*estSkinCoef;

% Plot for part 4
% How well does camera response predict coefficients for skin radiance data

%figure; plot(coef,estCoef,'+');
figure;
for i=1:size(SkinRadiance,2)
    plot(wavelength,SkinRadiance(:,i),'b');
    hold on; plot(wavelength, cameraPredictedSkinRadiance(:,i),'r');
    legend('Skin Radiance','Predicted Skin Radiance (From Camera)');
end
grid on;
title('Part 4:How well does camera response predict coefficients for skin radiance data');
%%%%%%%%%%%%%%%%%%%%%End of part 4%%%%%%%%%%%%%%%%%%%%%%%

