function [xhat,rmse,deltaE]=rectest(sensor,illuminant_name,snr,wave,savefilename)

% RECTEST simulates measurements of surfaces and estimates reflectances
%
%[xhat,rmse,deltaE]=rectest(sensor,illuminant_name,snr,wave,savefilename)
%
%INPUTS:
%   sensor:     Matrix giving the sensor's spectral sensitivity (quantum
%               efficiency).  Units are in photons
%               size(sensor)= num channels x length(wave)
%   illuminant_name:   String giving the name of the file where the
%                      illuminant is stored
%   snr:        Scalar giving the signal to noise ratio in dB of the
%               measurements  (20-40 is reasonable)
%   wave:        Vector giving the wavelength samples in nm
%   savefilename:  (Optional) String giving the name of the file 
%                  where results will be saved
%
%OUTPUTS:
%   xhat:        Matrix giving the estimates of the reflectances
%                size(xhat)=length(wave) x number of surfaces
%   rmse:        Matrix giving root mean squared error of estimate over all
%                wavelengths for each surface
%                size(rmse)=1 x number of surfaces
%   deltaE:      Vector giving CIELab delta E for each surface 
%                size(deltaE)=1 x number of surfaces

randn('state',0)    %introduced so random numbers will always be the same         

%% Load Data

%Load training data:
tmp=load('pig_reflectances');
reflectances=interp1(tmp.wavelength,tmp.data,wave);  %Adjust reflectances to desired wavelength sampling

%Load illuminant:
tmp=load(illuminant_name);
illuminant=interp1(tmp.wavelength,tmp.data,wave);   %Adjust illuminant to desired wavelength sampling

%Illuminants are stored in terms of energy per wavelength but 
%we need it in quanta (# of photons) per wavelength.
illuminant=Energy2Quanta(wave,illuminant')';


%% Simulate measurements

A=sensor*diag(illuminant);     %linear measurement system (combined sensor-illuminant)
ys_nonoise=A*reflectances;     %observed measurements for reconstruction
c=mean(ys_nonoise(:))/mean(ys_nonoise(:).^2)*10^(snr/10);   %scaling constant to achieve desired SNR
A=c*A;  %scale measurement matrix to get desired SNR using noise formula below, can think of as scaling the light level

ys_nonoise=A*reflectances;     %observed measurements for reconstruction

%The following noise formula is an approximation to photon shot noise.
%Actual shot noise is Poisson distributed but for large values can be
%modeled as an additive Gaussian noise with standard deviation given by the
%square root of the noise free value.
ys=ys_nonoise+sqrt(ys_nonoise).*randn(size(ys_nonoise));

%Formula for SNR to check it is similar to desired value:
% actual_snr=10*log10(mean(ys_nonoise.^2)/mean((ys_nonoise-ys).^2));   


%% Estimate reflectance from Measurements and Evaluate

%Following is estimated variance of noise for each measurement.  
%We could find the actual value but that is not possible in practice.
sigma2=abs(ys);

xhat=wienerfilter(reflectances,ys,A,sigma2);    %Estimated reflectances

[rmse,deltaE]=evaluateEstimates(reflectances,xhat,wave);    %Quality evaluation

if nargin==5 & ~isempty(savefilename)
    save(savefilename,'wave','sensor','illuminant_name','A','xhat','rmse','deltaE','c','snr')
end