% PSYCH 221
% Helmy Eltoukhy and Ali Ercan

% This is the main code for our project. Loads the previously measured and saved 
% variables and constructs the final images.

clear all;

load filters; % loads the power spectral densities of the filters
load qe; % loads the qe of the imager

% The wavelengths of qe is defined from 400 nm to 780 nm. Cut the psd's 
% of filters out of those wavelengths

specs(:,1) = speclight(6:101);

specs(:,2) = specyellow(6:101);

specs(:,3) = specred(6:101);

specs(:,4) = specgreen(6:101);

specs(:,5) = specblue(6:101);

specs(:,6) = specdaylight(6:101);

specs(:,7) = specbluegreen(6:101);

specs(:,8) = specorange(6:101);

% Now adjust the wavelengths such that they fit to wavelengths of the xyz 
% conversion matrix

wave = wave(1:83);
specs = specs(1:83,:);
qe = qe(1:83);
% As the values we get from spectrometer change with the strength of the light
% that shines on it, normalize the responses with the maximum of the pure light
% (without no filter)

maximum = max(specs(:,1));
specs = specs/maximum;

% Because spectum of different filters are measured with light, divide them 
% by the response of the light to get pure filter responses. 
% Below two lines are commented out in method 2 so that you have the light's 
% response in the filter response, you will be solving for surface reflectance 
% function only. 
light = specs(:,1);
for s= 1:8, specs(:,s) = specs(:,s)./light; end;

% Construct the matrix that multiplies with surface reflectance function 
% to give pixel value at the camera.
A = [];
filters = [2:5 7 8]; % The filters to be used. The first one is light it self,
        % not a filter, and 6th one is dayight filter, almost same response as 
        % light (with light), does not give any info.
for i = filters
    sp = specs(:,i);
    A = [A;qe.*sp'];
end

% Do the least norm for each point's surface reflectance function to find their 
% surface reflectance (SR)

load hdrImsNew; % Load the high DR images acquired in the camera. The last sample
    % before saturation algorithm is applied before saving hdrImsNew. 
load xyz; % Load the SR to XYZ conversion matrix
wavelength = wavelength(31:4:359); % Adjust the wavelengths to be compatible
xyz = xyz(31:4:359,:);             % with previous variables' wavelength range

[m n k] = size(im);
xyz_im = [];
psd=[];
Ad = pinv(A);
for i = 1:m
    for j = 1:n
        c = squeeze(im(i,j,filters)); % pixel i,j's HDR values for different filters
        psd(i,j,:) = (Ad*c); % least norm
        pix = xyz'*squeeze(psd(i,j,:)); % convert to XYZ
        xyz_im(i,j,:) = pix; % store the XYZ values for i,j
    end
end

xyz2rgb = [3.24 -1.54 -.5    % The XYZ to RGB conversion matrix we used. 
    -.97 1.88 .042 % This is the SRGB standard
    .056 .20 1.06];

for k = 1:m
    for l = 1:n
        pix = squeeze(xyz_im(k,l,:));
        pixrgb = xyz2rgb*pix;
        rgb_im(k,l,:) = pixrgb;
    end
end

rgb_im(find(rgb_im<0)) = 0; % Replace the negative RGB values with zero, 
                            % they do not make sense

rgb_mod = rgb_im/max(rgb_im(:)); %Scale the values so that everything in [0 1]

figure, imagesc(rgb_mod2); % Not color corrected RGB image (no gamma correction either)

% Do the color correction
% Find the mean R,G,B values for the white pathc in the color checker
r = mean2(squeeze(rgb_mod2(270:280,5:15,1)))
g = mean2(squeeze(rgb_mod2(270:280,5:15,2))) 
b = mean2(squeeze(rgb_mod2(270:280,5:15,3)))

% do the color correction
rgbcc(:,:,1) = rgb_mod2(:,:,1)/r;
rgbcc(:,:,2) = rgb_mod2(:,:,2)/g;
rgbcc(:,:,3) = rgb_mod2(:,:,3)/b;
rgbcc = rgbcc/max(rgbcc(:)); % everything has to be in [0 1] range for proper 
                            % display

% Do the gamma correction for both collor corrected and not color corrected image
% A gamma of 2.5 gave the best results for our screen
rgbgam =abs(1.055*rgb_mod.^(1/2.5)-.055);
testim(:,:,1) = rgbgam(:,:,1);
testim(:,:,2) = rgbgam(:,:,2);
testim(:,:,3) = rgbgam(:,:,3);
figure, imagesc(1/1*testim);

rgbgamcc =abs(1.055*rgbcc.^(1/2.5)-.055);
figure, imagesc(rgbgamcc);

