function [rmse,deltaE]=evaluateEstimates(reflectances,xhat,wave)

% EVALUATEESTIMATES returns RMSE and Delta E for reflectance estimates
%
%[rmse,deltaE]=evaluateEstimates(reflectances,xhat,wave)
%
%INPUTS:
%   reflectances: Matrix giving the data (columns are for a single surface)
%                size(reflectances)=length(wave) x number of surfaces
%   xhat:        Matrix giving the estimates of the reflectances
%                size(xhat)=length(wave) x number of surfaces
%   wave:        Vector giving the wavelength samples in nm
%
%OUTPUTS:
%   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


tmp=load('D65');
D65=interp1(tmp.wavelength,tmp.data,wave);  %D65 light

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

spectra1=reflectances'.*repmat(D65,size(reflectances,2),1); %estimated spectra under D65
spectra2=xhat'.*repmat(D65,size(xhat,2),1);     %predicted spectra under D65

deltaE=deltaEab(spectra1*XYZ,spectra2*XYZ,D65*XYZ); %deltaE color difference metric


rmse=sqrt(mean((reflectances-xhat).^2));