function CIEDE = CompareDuplicateColorMunkiMeasurements(ColorMunkiCSVfile1, ColorMunkiCSVfile2);
% Purpose		Compare two sets of ColorMunki measurements.
%
% Description   This routine is intended to test the consistency of two sets of measurements,
%				made by the ColorMunki spectrophotometer, of one set of colour samples.
%				It could also be used to measure the differences in two sets of 
%				corresponding samples, under the assumption that the measurement device and
%				process is consistent.  
%
%				The input consists of two files of ColorMunki measurements, that have been 
%				exported to a comma-separated (.csv) file format.  The files are presumed to
%				consist of duplicate measurements of the same set of colour samples, in the
%				same order.  Measurement errors can result from ColorMunki variability, 
%				differences in measuring procedures, or even changes in the samples (if the
%				set of samples is measured at different times).  If all these factors have
%				been accounted for, then the corresponding sample measurements should not
%				differ by much.  The routine returns a vector of differences, expressed with
%				respect to CIE DE 2000.  A histogram, summary statistics, and a list of 
%				largest differences can also be displayed, if desired.
%
%				Since this routine was written for analysis of the Munsell system, the CIE
%				coordinates for the samples are calculated with respect to Illuminant C.  In
%				many applications, Illuminant D50 or D65 might be more appropriate.  For most
%				samples, however, the DE values calculated should not vary much with the
%				illuminant.
%
%				ColorMunkiCSVfile1, ColorMunkiCSVfile2		File names (as text strings) 
%						for two ColorMunki .csv files containing duplicate measurements of
%						the same set of samples.
%
%				CIEDE	A vector of the DE2000 differences between duplicate measurements
%						of the same set of samples.
%
% Author		Paul Centore (September 29, 2012)
% Revised		Paul Centore (August 31, 2013)  
%				 ---Moved from MunsellToolbox program to MunsellAndKubelkaMunkToolbox.
% Revised		Paul Centore (January 1, 2014)
%				 ---Calculated white point for Illuminant C and 2 deg observer, and passed to revised
%				    routine for calculating CIE DE 2000
%				 ---Replaced called routine ReflectanceSpectrumToCIEIllumC with routine
%				    ReflectancesToCIEwithWhiteY100, to insure consistency in CIE calculations
% Revised		Paul Centore (December 2, 2014)
%				 ---Changed -deps option to -depsc to preserve colour in .eps files
%
% Copyright 2012, 2014 Paul Centore
%
%    This file is part of MunsellAndKubelkaMunkToolbox.
%
%    MunsellAndKubelkaMunkToolbox is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    MunsellAndKubelkaMunkToolbox is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with MunsellAndKubelkaMunkToolbox.  If not, see <http://www.gnu.org/licenses/>.

% Set a flag to produce graphical output if desired
DisplayHistogram = false						;

% Open the first CSV file generated by the ColorMunki Design
fid  = fopen(ColorMunkiCSVfile1, 'r')		;

% Read the file into a matrix; strings in the file will not appear as strings, but
% numbers will be preserved.
AllFileData = dlmread(fid, ',')				;
fclose(fid)									;

[row, col] = size(AllFileData)				;

% The first row of the file is a header row.  The first four entries of the header
% row are Name, L*, a*, and b*.  The remainder of the header row is a list of 
% wavelengths at which reflectances are measured.  Extract these wavelengths.
Wavelengths = AllFileData(1, 5:col)			;

% Every row of the file except the first is data for one colour sample.  The first
% four entries of each row are a name, an L* value, an a* value, and a b* value.  The
% remaining entries are the reflectances for the wavelengths listed in the first row.
Reflectances = AllFileData(2:row, 5:col)	;
% Find CIE coordinates for the first file (the following line was modified Jan. 1, 2014)
CIEcoords1 = ReflectancesToCIEwithWhiteY100(Wavelengths, Reflectances, 'C/2');

% Perform the same operations on the second file
fid          = fopen(ColorMunkiCSVfile2, 'r')	;
AllFileData  = dlmread(fid, ',')				;
fclose(fid)										;
[row, col]   = size(AllFileData)				;
Wavelengths  = AllFileData(1, 5:col)			;
Reflectances = AllFileData(2:row, 5:col)		;
% The following line was modified Jan. 1, 2014
CIEcoords2   = ReflectancesToCIEwithWhiteY100(Wavelengths, Reflectances, 'C/2')	;

% Now that the two files have been read in and converted to CIE coordinates,
% a vector of CIE DE2000 values can be computed for corresponding samples.
[rows1,~] = size(CIEcoords1)		;
[rows2,~] = size(CIEcoords2)		;
CIEDE     = []						;
WhitePointXYZ = WhitePointWithYEqualTo100('C/2')	;	% Line added Jan. 1, 2014, to calculate white point
for ind = 1:min(rows1,rows2)
	XYZ1   = CIEcoords1(ind,1:3)		;
	XYZ2   = CIEcoords2(ind,1:3)		;
	DE2000 = CIEDE2000ForXYZ(XYZ1, XYZ2, WhitePointXYZ); % Modified Dec. 14, 2013
	CIEDE  = [CIEDE, DE2000]			;
end

if DisplayHistogram
    EdgeVector = [0:0.1:5]					;
	figure
    Counts = histc(CIEDE, EdgeVector)		;
	HistogramData = [EdgeVector; Counts]	;
	stairs(EdgeVector, Counts, 'k')				;
    set(gca, 'xlim', [0,5])					;
%	set(gca, 'xtick', [0:0.5:3])
%	set(gca, 'xticklabel', ['0.0'; '0.5'; '1.0'; '1.5'; '2.0'; '2.5'; '3.0'])
	set(gca, 'xtick', [0:0.5:5])
	set(gca, 'xticklabel', ['0.0'; '0.5'; '1.0'; '1.5'; '2.0'; '2.5'; '3.0'; '3.5'; '4.5'; '4.5'; '5.0'])
    set(gcf, 'Name', ['Differences in Measurements']);
    figname = ['MeasurementDifferences']	;
    print(gcf, [figname,'.eps'], '-depsc')	;
    print(gcf, [figname,'.png'], '-dpng')	;
    print(gcf, [figname,'.jpg'], '-djpg')	;
    print(gcf, [figname,'.pdf'], '-dpdf')	;
	disp(['Number of sample pairs:                         ', num2str(length(CIEDE))])	;
	disp(['Mean, median, minimum, and maximum difference:  ', ...
	      num2str(mean(CIEDE)),', ',...
	      num2str(median(CIEDE)),', ',...
		  num2str(min(CIEDE)),', ',...
	      num2str(max(CIEDE))])			;
		  
    % List most extreme differences
    [SortedDEs, idx] = sort(transpose(CIEDE), 'descend')			;		
	WorstCases       = [transpose(1:length(idx)), idx, SortedDEs]	;
	disp([' ']);
	disp(['      Worst Cases             '])					;
	disp(['   File Position    DE'])					;
	for ctr = 1:min([rows1, rows2, 30])
	    disp([num2str(WorstCases(ctr,1)),'         ', ...
		      num2str(WorstCases(ctr,2)),'         ', ...
			  num2str(WorstCases(ctr,3))])	;
	end
	disp([' ']);
end