clc

clear

tic

global wave
global illuminant_name
global snr


global num_filters

global lb
global ub
global mean_range


%Desired sampling of wavelength in nm to test over:
wave=400:1:700;
%With the pig data you can probably go up to 800nm but not more because
%then it gets noisy.




%Specify sensor
%size(sensor)=number of channels x length(wave)
%The absolute scale of the sensor matrix is irrelevant.
%Units should be quanta so all values should be between 0 and 1.
tmp=load('NikonD100');
%sensor=interp1(tmp.wavelength,tmp.data,wave)';

figure
plot(tmp.wavelength,tmp.data)
title('NikonD100 SPD')
xlabel('Wavelength (nm)')
ylabel('Transmission')


values = load('pig_reflectances');

figure

original_spectra = interp1(values.wavelength,values.data,wave);
%plot(values.wavelength,values.data)
plot(wave,original_spectra)
title('Original Spectra')
xlabel('Wavelength (nm)')
ylabel('Reflectance')

% Specify illuminant
%You can choose other illuminant if desired.  Tungsten might be good since
%it has a lot of energy in the red region of the spectrum, just like tissue.
illuminant_name='Tungsten';



snrs=20:5:40;     %Reasonable values are 20-40.
% Specify SNR
snr=snrs(1)


options = optimset('MaxFunEvals',10^3, 'PlotFcns',@optimplotfval,'TolFun',1e-9,'TolCon',1e-09);


% possibilities
% blue_mean_range = 400:1:700;
% green_mean_range = 400:1:700;
% red_mean_range = 400:1:700;

mean_range = 400:1:700;


num_filters = 4;




switch num_filters
    
    
    
    
    case 2
        
        num_trials = 50;
        
        for counter = 1:num_trials
            
            
            % randomly select "means"
            filter1_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            filter2_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            
            
            % selecting standard deviation for 1st filter
            if filter1_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter1_std = (filter1_mean - 398)/2;
                filter1_std = round(1 + (max_filter1_std-1)*rand(1));
                
            else
                max_filter1_std = 100;
                
                filter1_std = round(1 + (max_filter1_std-1)*rand(1));
                
            end
            
            
            % selecting standard deviation for 2nd filter
            if filter2_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter2_std = (filter2_mean - 398)/2;
                
                filter2_std = round(1 + (max_filter2_std-1)*rand(1));
                
                
            else
                max_filter2_std = 100;
                
                filter2_std = round(1 + (max_filter2_std-1)*rand(1));
                
            end
            
            
            initial_cond = [filter1_mean;filter2_mean;filter1_std;filter2_std];
            
            initial_sensor_values{counter} = initial_cond;
            
            lb = [400;400;1;1];
            
            ub = [700;700;max_filter1_std;max_filter2_std];
            
            
            
            [optimum_sensor_values, error] = fmincon(@myfun,initial_cond,[],[],[],[],lb,ub,[],options);
            
            
            
            
            cost(counter) = error
            
            
            sensor_values{counter} = optimum_sensor_values;
            
            
        end
        
        % finding index for minimum error
        min_cost_index = find(cost == min(cost),1);
        
        
        % initial condition that gives the smallest error
        filter1_initial_mean = initial_sensor_values{min_cost_index}(1);
        filter2_initial_mean = initial_sensor_values{min_cost_index}(2);
        filter1_initial_std = initial_sensor_values{min_cost_index}(3);
        filter2_initial_std = initial_sensor_values{min_cost_index}(4);
        
        % selecting range for standard deviation of 1st filter
        if filter1_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter1_initial_std = (filter1_initial_mean - 398)/2;
         
        else
         
            max_filter1_initial_std = 100;           
        end
        
        
        % selecting range for standard deviation of 2nd filter
        if filter2_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter2_initial_std = (filter2_initial_mean - 398)/2;
       
        else
            max_filter2_initial_std = 100;
         
        end
        
        
        lb = [400;400;1;1];
        
        ub = [700;700;max_filter1_initial_std;max_filter2_initial_std];
        
        initial_cond = [filter1_initial_mean;filter2_initial_mean;filter1_initial_std;filter2_initial_std];
        [optimum_sensor_values, error] = fmincon(@myfun,initial_cond,[],[],[],[],lb,ub,[],options);
        
        
        
        
        % plotting initial functions - 2 filters
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        filter1_initial = (exp(-(wave - initial_sensor_values{min_cost_index}(1)).^2/(2*initial_sensor_values{min_cost_index}(3)^2)));
        filter2_initial =(exp(-(wave - initial_sensor_values{min_cost_index}(2)).^2/(2*initial_sensor_values{min_cost_index}(4)^2)));
        
        figure
        hold on
        
        plot(wave, filter1_initial,'-');
        plot(wave, filter2_initial,'--','LineWidth',1.01);
        
        
        title('Initial Guess for Filter SPD')
        ylabel('Transmission')
        xlabel('Wavelength (nm)')
        
        hold off
        
        
        
        % plotting optimized functions - 2 filters
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        filter1_optimized = (exp(-(wave - sensor_values{min_cost_index}(1)).^2/(2*sensor_values{min_cost_index}(3)^2)));
        filter2_optimized =(exp(-(wave - sensor_values{min_cost_index}(2)).^2/(2*sensor_values{min_cost_index}(4)^2)));
        
        
        figure
        hold on
        
        plot(wave, filter1_optimized,'-');
        plot(wave, filter2_optimized,'--','LineWidth',1.01);
        
        
        title('Optimized Filter SPD')
        ylabel('Transmission')
        xlabel('Wavelength (nm)')
        
        hold off
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        
        % optimal sensor - 2 filters
        optimal_sensor = vertcat(filter1_optimized,filter2_optimized);
        [xhat,rmse,deltaE]=rectest(optimal_sensor,illuminant_name,snr,wave);
        
        
        
        
        figure
        plot(wave,xhat)
        title('Estimated Spectra')
        ylabel('Reflectance')
        xlabel('Wavelength (nm)')
        
        
        figure
        plot(deltaE)
        title('\DeltaE for each Surface')
        xlabel('Surface Number')
        ylabel('\DeltaE')
        
        
        min_deltaE_surface = find(deltaE == min(deltaE));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,min_deltaE_surface))
        plot(wave,xhat(:,min_deltaE_surface),'--')
        title('Surface with Minimum \DeltaE')
        xlabel('Wavelength (nm)')
        ylabel('\DeltaE')
        
        
        hold off
        
        
        max_deltaE_surface = find(deltaE == max(deltaE));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,max_deltaE_surface))
        plot(wave,xhat(:,max_deltaE_surface),'--')
        title('Surface with Maximum \DeltaE')
        xlabel('Wavelength (nm)')
        ylabel('\DeltaE')
        
        hold off
        
        
        
        
        
        
        
        
        
        figure
        plot(rmse)
        title('RMSE for each Surface')
        xlabel('Surface Number')
        ylabel('RMSE')
        
        
        
        min_rmse_surface = find(rmse == min(rmse));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,min_rmse_surface))
        plot(wave,xhat(:,min_rmse_surface),'--')
        title('Surface with Minimum RMSE')
        xlabel('Wavelength (nm)')
        ylabel('RMSE')
        
        hold off
        
        
        
        max_rmse_surface = find(rmse == max(rmse));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,max_rmse_surface))
        plot(wave,xhat(:,max_rmse_surface),'--')
        title('Surface with Maximum RMSE')
        xlabel('Wavelength (nm)')
        ylabel('RMSE')
        
        hold off
        
        
        
        
        
        % going through each wavelength for all surfaces
        for counter = 1:size(original_spectra,1)
            
            RMSE_lambda(counter) = sqrt(mean((original_spectra(counter,:) - xhat(counter,:)).^2));
            
        end
        
        figure
        
        hold on
        
        [AX H H1] = plotyy(wave,RMSE_lambda,wave,filter1_optimized)
        [AX H H2] = plotyy(wave,RMSE_lambda,wave,filter2_optimized)
       
        
        title('RMSE & Transmission vs. Wavelength')
        xlabel('Wavelength (nm)')
        
        
        set(get(AX(1),'Ylabel'),'String','RMSE','color','r')
        set(AX(1),'Ycolor','r')
        set(get(AX(2),'Ylabel'),'String','Transmission','color','b')
        set(AX(2),'Ycolor','b')
        
        
        
        
        set(H,'LineStyle','-','color','r')
        set(H1,'LineStyle','-','color','b')
        set(H2,'LineStyle','--','color','b','LineWidth',1.01)
        
        
        hold off
        
        
        
        
        
        
        
        
        
        
        
        
    case 3
        
        num_trials = 50;
        
        for counter = 1:num_trials
            
            
            
            
            
            % randomly select "means"
            filter1_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            filter2_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            filter3_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            
            
            if filter1_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter1_std = (filter1_mean - 398)/2;
                filter1_std = round(1 + (max_filter1_std-1)*rand(1));
                
            else
                max_filter1_std = 100;
                
                filter1_std = round(1 + (max_filter1_std-1)*rand(1));
                
            end
            
            
            
            if filter2_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter2_std = (filter2_mean - 398)/2;
                
                filter2_std = round(1 + (max_filter2_std-1)*rand(1));
                
                
            else
                max_filter2_std = 100;
                
                filter2_std = round(1 + (max_filter2_std-1)*rand(1));
                
            end
            
            
            if filter3_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter3_std = (filter3_mean - 398)/2;;
                filter3_std = round(1 + (max_filter3_std-1)*rand(1));
                
            else
                max_filter3_std = 100;
                
                filter3_std = round(1 + (max_filter3_std-1)*rand(1));
                
            end
            
            
            
            
            
            
            initial_cond = [filter1_mean;filter2_mean;filter3_mean;filter1_std;filter2_std;filter3_std];
            
            initial_sensor_values{counter} = initial_cond;
            
            lb = [400;400;400;1;1;1];
            
            ub = [700;700;700;max_filter1_std;max_filter2_std;max_filter3_std];
            
            
            
            [optimum_sensor_values, error] = fmincon(@myfun,initial_cond,[],[],[],[],lb,ub,[],options);
            
            
            
            
            cost(counter) = error
            
            
            sensor_values{counter} = optimum_sensor_values;
            
            
        end
        
        
        % finding index for minimum error
        min_cost_index = find(cost == min(cost),1);
        
        
        % initial condition that gives the smallest error
        filter1_initial_mean = initial_sensor_values{min_cost_index}(1);
        filter2_initial_mean = initial_sensor_values{min_cost_index}(2);
        filter3_initial_mean = initial_sensor_values{min_cost_index}(3);
        filter1_initial_std = initial_sensor_values{min_cost_index}(4);
        filter2_initial_std = initial_sensor_values{min_cost_index}(5);
        filter3_initial_std = initial_sensor_values{min_cost_index}(6);
        
        % selecting range for standard deviation of 1st filter
        if filter1_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter1_initial_std = (filter1_initial_mean - 398)/2;
         
        else
         
            max_filter1_initial_std = 100;           
        end
        
        
        % selecting range for standard deviation of 2nd filter
        if filter2_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter2_initial_std = (filter2_initial_mean - 398)/2;
            
        else
            max_filter2_initial_std = 100;
            
        end
        
        
        if filter3_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter3_initial_std = (filter3_initial_mean - 398)/2;
         
        else
            max_filter3_initial_std = 100;
        
        end
        
        
        lb = [400;400;400;1;1;1];
        
        ub = [700;700;700;max_filter1_initial_std;max_filter2_initial_std;max_filter3_initial_std];
        
        initial_cond = [filter1_initial_mean;filter2_initial_mean;filter3_initial_mean;filter1_initial_std;filter2_initial_std;filter3_initial_std];
        [optimum_sensor_values, error] = fmincon(@myfun,initial_cond,[],[],[],[],lb,ub,[],options);
        
        
        
        
        
        % plotting initial functions - 3 filters
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        filter1_initial = (exp(-(wave - initial_sensor_values{min_cost_index}(1)).^2/(2*initial_sensor_values{min_cost_index}(4)^2)));
        filter2_initial =(exp(-(wave - initial_sensor_values{min_cost_index}(2)).^2/(2*initial_sensor_values{min_cost_index}(5)^2)));
        filter3_initial = (exp(-(wave - initial_sensor_values{min_cost_index}(3)).^2/(2*initial_sensor_values{min_cost_index}(6)^2)));
        
        
        figure
        hold on
        
        plot(wave, filter1_initial,'-');
        plot(wave, filter2_initial,'--','LineWidth',1.01);
        plot(wave, filter3_initial,'-.');
        title('Initial Guess for Filter SPD')
        ylabel('Transmission')
        xlabel('Wavelength (nm)')
        
        hold off
        % plotting optimized functions - 3 filters
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        filter1_optimized = (exp(-(wave - sensor_values{min_cost_index}(1)).^2/(2*sensor_values{min_cost_index}(4)^2)));
        filter2_optimized =(exp(-(wave - sensor_values{min_cost_index}(2)).^2/(2*sensor_values{min_cost_index}(5)^2)));
        filter3_optimized = (exp(-(wave - sensor_values{min_cost_index}(3)).^2/(2*sensor_values{min_cost_index}(6)^2)));
        
        figure
        hold on
        
        plot(wave, filter1_optimized,'-');
        plot(wave, filter2_optimized,'--','LineWidth',1.01);
        plot(wave, filter3_optimized,'-.');
        title('Optimized Filter SPD')
        ylabel('Transmission')
        xlabel('Wavelength (nm)')
        
          
        
        hold off
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        
        % optimal sensor - 3 filters
        optimal_sensor = vertcat(filter1_optimized,filter2_optimized,filter3_optimized);
        [xhat,rmse,deltaE]=rectest(optimal_sensor,illuminant_name,snr,wave);
        
        
        
        
        figure
        plot(wave,xhat)
        title('Estimated Spectra')
        ylabel('Reflectance')
        xlabel('Wavelength (nm)')
        
        
        figure
        plot(deltaE)
        title('\DeltaE for each Surface')
        xlabel('Surface Number')
        ylabel('\DeltaE')
        
        
        min_deltaE_surface = find(deltaE == min(deltaE));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,min_deltaE_surface))
        plot(wave,xhat(:,min_deltaE_surface),'--')
        title('Surface with Minimum \DeltaE')
        xlabel('Wavelength (nm)')
        ylabel('\DeltaE')
        
        hold off
        
        
        max_deltaE_surface = find(deltaE == max(deltaE));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,max_deltaE_surface))
        plot(wave,xhat(:,max_deltaE_surface),'--')
        title('Surface with Maximum \DeltaE')
        xlabel('Wavelength (nm)')
        ylabel('\DeltaE')
        
        hold off
        
        
        
        
        
        
        
        
        figure
        plot(rmse)
        title('RMSE for each Surface')
        xlabel('Surface Number')
        ylabel('RMSE')
        
        min_rmse_surface = find(rmse == min(rmse));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,min_rmse_surface))
        plot(wave,xhat(:,min_rmse_surface),'--')
        title('Surface with Minimum RMSE')
        xlabel('Wavelength (nm)')
        ylabel('RMSE')
        
        hold off
        
        
        
        max_rmse_surface = find(rmse == max(rmse));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,max_rmse_surface))
        plot(wave,xhat(:,max_rmse_surface),'--')
        title('Surface with Maximum RMSE')
        xlabel('Wavelength (nm)')
        ylabel('RMSE')
        
        hold off
        
        
        
        
        
        % going through each wavelength for all surfaces
        for counter = 1:size(original_spectra,1)
            
            RMSE_lambda(counter) = sqrt(mean((original_spectra(counter,:) - xhat(counter,:)).^2));
            
        end
        
        figure
        
        hold on
        
        [AX H H1] = plotyy(wave,RMSE_lambda,wave,filter1_optimized)
        [AX H H2] = plotyy(wave,RMSE_lambda,wave,filter2_optimized)
        [AX H H3] = plotyy(wave,RMSE_lambda,wave,filter3_optimized)
        
        
        title('RMSE & Transmission vs. Wavelength')
        xlabel('Wavelength (nm)')
        
        
        set(get(AX(1),'Ylabel'),'String','RMSE','color','r')
        set(AX(1),'Ycolor','r')
        set(get(AX(2),'Ylabel'),'String','Transmission','color','b')
        set(AX(2),'Ycolor','b')
        
       
        
        set(H,'LineStyle','-','color','r')
        set(H1,'LineStyle','-','color','b')
        set(H2,'LineStyle','--','color','b','LineWidth',1.01)
        set(H3,'LineStyle','-.','color','b')
       
        
        
        
        hold off
        
        
        
        
        
        
        
        
        
    case 4
        
        
        
        
        
        num_trials = 50;
        
        for counter = 1:num_trials
            
            
            
            
            
            % randomly select "means"
            filter1_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            filter2_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            filter3_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            filter4_mean = round(mean_range(1) + (mean_range(length(mean_range))-mean_range(1))*rand(1));
            
            if filter1_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter1_std = (filter1_mean - 398)/2;
                filter1_std = round(1 + (max_filter1_std-1)*rand(1));
                
            else
                max_filter1_std = 100;
                
                filter1_std = round(1 + (max_filter1_std-1)*rand(1));
                
            end
            
            
            
            if filter2_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter2_std = (filter2_mean - 398)/2;
                
                filter2_std = round(1 + (max_filter2_std-1)*rand(1));
                
                
            else
                max_filter2_std = 100;
                
                filter2_std = round(1 + (max_filter2_std-1)*rand(1));
                
            end
            
            
            if filter3_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter3_std = (filter3_mean - 398)/2;
                filter3_std = round(1 + (max_filter3_std-1)*rand(1));
                
            else
                max_filter3_std = 100;
                
                filter3_std = round(1 + (max_filter3_std-1)*rand(1));
                
            end
            
            
            if filter4_mean < (mean_range(1) + mean_range(length(mean_range)))/2
                
                max_filter4_std = (filter4_mean - 398)/2;
                filter4_std = round(1 + (max_filter4_std-1)*rand(1));
                
            else
                max_filter4_std = 100;
                
                filter4_std = round(1 + (max_filter4_std-1)*rand(1));
                
            end
            
            
            
            
            
            initial_cond = [filter1_mean;filter2_mean;filter3_mean;filter4_mean;filter1_std;filter2_std;filter3_std;filter4_std];
            
            initial_sensor_values{counter} = initial_cond;
            
            lb = [400;400;400;400;1;1;1;1];
            
            ub = [700;700;700;700;max_filter1_std;max_filter2_std;max_filter3_std;max_filter4_std];
            
            
            
            [optimum_sensor_values, error] = fmincon(@myfun,initial_cond,[],[],[],[],lb,ub,[],options);
            
            
            
            
            cost(counter) = error
            
            
            sensor_values{counter} = optimum_sensor_values;
            
            
        end
        
        
        % finding index for minimum error
        min_cost_index = find(cost == min(cost),1);
        
        
        
        % initial condition that gives the smallest error
        filter1_initial_mean = initial_sensor_values{min_cost_index}(1);
        filter2_initial_mean = initial_sensor_values{min_cost_index}(2);
        filter3_initial_mean = initial_sensor_values{min_cost_index}(3);
        filter4_initial_mean = initial_sensor_values{min_cost_index}(4);
        filter1_initial_std = initial_sensor_values{min_cost_index}(5);
        filter2_initial_std = initial_sensor_values{min_cost_index}(6);
        filter3_initial_std = initial_sensor_values{min_cost_index}(7);
        filter4_initial_std = initial_sensor_values{min_cost_index}(8);
        
        % selecting range for standard deviation of 1st filter
        if filter1_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter1_initial_std = (filter1_mean - 398)/2;
         
        else
         
            max_filter1_initial_std = 100;           
        end
        
        
        % selecting range for standard deviation of 2nd filter
        if filter2_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter2_initial_std = (filter2_initial_mean - 398)/2;
            
        else
            max_filter2_initial_std = 100;
            
        end
        
        
        if filter3_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter3_initial_std = (filter3_initial_mean - 398)/2;
         
        else
            max_filter3_initial_std = 100;
            
        end
        
        if filter4_initial_mean < (mean_range(1) + mean_range(length(mean_range)))/2
            
            max_filter4_initial_std = (filter4_initial_mean - 398)/2;
          
        else
            max_filter4_initial_std = 100;
        
        end
        
        
        lb = [400;400;400;400;1;1;1;1];
        
        ub = [700;700;700;700;max_filter1_initial_std;max_filter2_initial_std;max_filter3_initial_std;max_filter4_initial_std];
        
        initial_cond = [filter1_initial_mean;filter2_initial_mean;filter3_initial_mean;filter4_initial_mean;filter1_initial_std;filter2_initial_std;filter3_initial_std;filter4_initial_std];
        [optimum_sensor_values, error] = fmincon(@myfun,initial_cond,[],[],[],[],lb,ub,[],options);
        
        
        
        
        
        
        % plotting initial functions - 4 filters
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        filter1_initial = (exp(-(wave - initial_sensor_values{min_cost_index}(1)).^2/(2*initial_sensor_values{min_cost_index}(5)^2)));
        filter2_initial =(exp(-(wave - initial_sensor_values{min_cost_index}(2)).^2/(2*initial_sensor_values{min_cost_index}(6)^2)));
        filter3_initial = (exp(-(wave - initial_sensor_values{min_cost_index}(3)).^2/(2*initial_sensor_values{min_cost_index}(7)^2)));
        filter4_initial = (exp(-(wave - initial_sensor_values{min_cost_index}(4)).^2/(2*initial_sensor_values{min_cost_index}(8)^2)));
        
        
        figure
        hold on
        
        plot(wave, filter1_initial,'-');
        plot(wave, filter2_initial,'--','LineWidth',1.01);
        plot(wave, filter3_initial,'-.');
        plot(wave, filter4_initial,':');
        title('Initial Guess for Filter SPD')
        ylabel('Transmission')
        xlabel('Wavelength (nm)')
        
        hold off
        
        % plotting optimized functions - 4 filters
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        filter1_optimized = (exp(-(wave - sensor_values{min_cost_index}(1)).^2/(2*sensor_values{min_cost_index}(5)^2)));
        filter2_optimized =(exp(-(wave - sensor_values{min_cost_index}(2)).^2/(2*sensor_values{min_cost_index}(6)^2)));
        filter3_optimized = (exp(-(wave - sensor_values{min_cost_index}(3)).^2/(2*sensor_values{min_cost_index}(7)^2)));
        filter4_optimized = (exp(-(wave - sensor_values{min_cost_index}(4)).^2/(2*sensor_values{min_cost_index}(8)^2)));
        
        
        figure
        hold on
        
        plot(wave, filter1_optimized,'-');
        plot(wave, filter2_optimized,'--','LineWidth',1.01);
        plot(wave, filter3_optimized,'-.');
        plot(wave, filter4_optimized,':');
        title('Optimized Filter SPD')
        ylabel('Transmission')
        xlabel('Wavelength (nm)')
        
        hold off
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        
        % optimal sensor - 4 filters
        optimal_sensor = vertcat(filter1_optimized,filter2_optimized,filter3_optimized,filter4_optimized);
        [xhat,rmse,deltaE]=rectest(optimal_sensor,illuminant_name,snr,wave);
        
        
        
        figure
        plot(wave,xhat)
        title('Estimated Spectra')
        ylabel('Reflectance')
        xlabel('Wavelength (nm)')
        
        
        figure
        plot(deltaE)
        title('\DeltaE for each Surface')
        xlabel('Surface Number')
        ylabel('\DeltaE')
        
        
        min_deltaE_surface = find(deltaE == min(deltaE));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,min_deltaE_surface))
        plot(wave,xhat(:,min_deltaE_surface),'--')
        title('Surface with Minimum \DeltaE')
        xlabel('Wavelength')
        ylabel('\DeltaE')
        
        hold off
        
        
        max_deltaE_surface = find(deltaE == max(deltaE));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,max_deltaE_surface))
        plot(wave,xhat(:,max_deltaE_surface),'--')
        title('Surface with Maximum \DeltaE')
        xlabel('Wavelength (nm)')
        ylabel('\DeltaE')
        
        hold off
        
        
        
        
        
        figure
        plot(rmse)
        title('RMSE for each Surface')
        xlabel('Surface Number')
        ylabel('RMSE')
        
        
        min_rmse_surface = find(rmse == min(rmse));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,min_rmse_surface))
        plot(wave,xhat(:,min_rmse_surface),'--')
        title('Surface with Minimum RMSE')
        xlabel('Wavelength (nm)')
        ylabel('RMSE')
        
        hold off
        
        
        
        max_rmse_surface = find(rmse == max(rmse));
        
        figure
        
        hold on
        plot(wave,original_spectra(:,max_rmse_surface))
        plot(wave,xhat(:,max_rmse_surface),'--')
        title('Surface with Maximum RMSE')
        xlabel('Wavelength (nm)')
        ylabel('RMSE')
        
        hold off
        
        
        
        
        
        % going through each wavelength for all surfaces
        for counter = 1:size(original_spectra,1)
            
            RMSE_lambda(counter) = sqrt(mean((original_spectra(counter,:) - xhat(counter,:)).^2));
            
        end
        
        figure
        
        hold on
        
        [AX H H1] = plotyy(wave,RMSE_lambda,wave,filter1_optimized)
        [AX H H2] = plotyy(wave,RMSE_lambda,wave,filter2_optimized)
        [AX H H3] = plotyy(wave,RMSE_lambda,wave,filter3_optimized)
        [AX H H4] = plotyy(wave,RMSE_lambda,wave,filter4_optimized)
        
        title('RMSE & Transmission vs. Wavelength')
        xlabel('Wavelength (nm)')
        
        
        set(get(AX(1),'Ylabel'),'String','RMSE','color','r')
        set(AX(1),'Ycolor','r')
        set(get(AX(2),'Ylabel'),'String','Transmission','color','b')
        set(AX(2),'Ycolor','b')
        
        set(H,'LineStyle','-','color','r')
        set(H1,'LineStyle','-','color','b')
        set(H2,'LineStyle','--','color','b','LineWidth',1.01)
        set(H3,'LineStyle','-.','color','b')
        set(H4,'LineStyle',':','color','b')
        
         
        hold off
        
       
        
end































%{


tmp=load(illuminant_name);
illuminant=interp1(tmp.wavelength,tmp.data,wave);  %illuminant

tmp=load('XYZ');
XYZ=interp1(tmp.wavelength,tmp.data,wave);  %XYZ standard observer


%spectra = xhat'.*repmat(illuminant,size(xhat,2),1);

spectra = xhat';



spectra_XYZ = spectra * XYZ;


xy = chromaticity(spectra_XYZ);


%Here xs and ys holds the chromaticity coordinates for each material.

xs = xy(:,1);
ys = xy(:,2);



background = 'white';
chromaticityDiagram(background)









for inputnum=1:size(xy,1)
    
    
    hold on
    
    
    if rmse(inputnum) < 0.01
        
        % white
        color_value=  [1 1 1];
        
        plot(xs(inputnum),ys(inputnum),'ko','MarkerSize',7,'color',color_value)
        
    else
        
        %black
        color_value = [0 0 0];
        
        plot(xs(inputnum),ys(inputnum),'ko','MarkerSize',7,'color',color_value)
        
    end
    
    
    
    
    
    
end

%}

toc






