function MCDMs = MCDMsFromMultipleSamples( ...
							ListOfColorMunkiCSVfiles, 	...
							OutputDirectory,			...
							NameKey,					...
							IllumObs);
% Purpose		Calculate a set of mean colour differences from the mean (MCDMs), with   
%				respect to CIEDE 2000.  The MCDM is a simple one-dimensional measure of the
%				variability in a spectrophotometer s object-colour measurements.  
%
% Description   Repeated spectrophotometer measurements of an object-colour sample (i.e.
%				measurements made with the identical spectrophotometer, of the identical
%				sample, but at different times) will show some variability.  For practical
%				purposes, it is good to have an estimate of this variability.  Standards
%				such as ASTM E2214 give some complicated, multi-dimensional measures for
%				this "repeatability."  The MCDM is a simpler, one-dimensional measure which
%				is easy to interpret.
%
%				The MCDM is used when the reflectance spectrum for one sample is measured 
%				multiple times.  The reflectance spectra measurements will vary.  As long
%				as the instrument is unbiased in each wavelength, the best estimate of the 
%				true reflectance spectrum is the average, or mean, of the reported spectra.
%				Each reported spectrum can then be compared to the mean spectrum, by using
%				a colour difference expression such as CIEDE 2000.  The result
%				will be a set of colour differences between the reported spectra and the
%				mean spectrum.  The mean of these colour differences is the MCDM.
%
%				Strictly speaking, an MCDM applies to one sample, that is measured multiple
%				times.  To evaluate a spectrophotometer s repeatability, the MCDM should
%				be calculated for a set of samples; preferably the set should cover a wide
%				gamut of colours, fairly uniformly.  The result will be a set of MCDMs.
%				The set can be plotted in a histogram, and summary statistics can be 
%				calculated, to describe the spectrophotometer s repeatability.
%
%				This routine produces just such output.  The input is a set of files, all
%				in the .csv format used by X-Rite s ColorMunki spectrophotometer.  The header
%				line of each file gives the wavelengths at which reflectances were measured.
%				Each subsequent line gives reflectances for one sample, at those wavelengths.
%				Each file was generated by measuring the same set of samples, in the same
%				order, so the ith row of the file will refer to the ith sample, across all
%				files.
%
%				The routine goes through the samples one by one.  Each sample will have been
%				measured n times (if there are n files), and there will be a different 
%				reflectance spectrum in each file.  Another routine is called to calculate
%				the MCDM for that sample.  That MCDM is added to a running list of MCDMs.
%				After producing a histogram and some summary statistics, the routine returns
%				the complete list of MCDMs. 
%
% Inputs		ListOfColorMunkiCSVfiles	A list (using {}) of strings, which are the filenames
%											of the ColorMunki .csv files to be averaged.  Each
%											file should have the same number of samples, in the
%											same order, and should be in ColorMunki .csv format
%
%				OutputDirectory		The directory where figures and other output will be saved
%
%				NameKey		A string that is appended to files and figures before they are saved
%
%				IllumObs	An illuminant/observer string, such as 'D50/2' or 'F12_64.' The
%							first part is a standard illuminant designation.  The symbol in the
%							middle can be either a forward slash or an underscore.  The number at
%							the end is either 2, 10, 31, or 64.  2 and 31 both correspond to the
%							standard 1931 2 degree observer, while 10 and 64 both correspond to 
%							the standard 1964 10 degree observer.  If no observer is indicated, then
%							the 1931 observer is assumed.  The input IllumObs is optional.
%
% Author		Paul Centore (January 5, 2015)
% Revision		Paul Centore (May 28, 2015)
%				----Added option to print out MCDM sequence, and spectra for high MCDMs
%				----Added NameKey and OutputDirectory as new inputs
%
% Copyright 2015 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/>.

% Default to Illuminant C with a 2 degree standard observer if no illuminant or observer 
% information is input
if ~exist('IllumObs')
    IllumObs = 'C/2'	;
end

DrawHistogram = false	;	% Flag for displaying histogram of MCDMs
DrawSequence  = true	;	% Flag for displaying MCDMs as a sequence over an ordinal x-scale
PlotOutliers  = true	;	% Flag for plotting reflectance spectra of individual samples
							% 	with MCDMs above an assigned threshold
MCDMThreshold = 2.0		;	% Plot reflectance spectra for samples whose MCDMs exceed this
							% 	threshold							

% Set formatting for some figures
AxisFontSizeInPoints  = 14	;
LabelFontSizeInPoints = 16	;

NumberOfFiles = length(ListOfColorMunkiCSVfiles)	;
% The files should all have a common set of wavelengths, so extract those wavelengths,
% using the first file
[Wavelengths, ~] = ColorMunkiCSVfileToOctaveFormat(ListOfColorMunkiCSVfiles{1});
% To avoid repeated file openings, extract all the reflectance information from each file
% once, and store it locally in a list of matrices.  
AllReflectances = {}	;
for ctr = 1:NumberOfFiles
	[~, Reflectances] = ColorMunkiCSVfileToOctaveFormat(ListOfColorMunkiCSVfiles{ctr});
	AllReflectances{ctr} = Reflectances		;
end
% Find the number of samples, using one of the reflectance matrices
TempReflectances = AllReflectances{1}		;
NumberOfSamples = size(TempReflectances,1)	;

% Calculate the MCDM for each sample, and save it in a master list
MCDMs = []	;
for ctr = 1:NumberOfSamples
	% Extract the reflectance spectra for one sample.  Each spectrum is in a different file
	SampleReflectanceSpectra = []	;
	for filectr = 1:NumberOfFiles
		TempReflectances         = AllReflectances{filectr}								;
		SampleReflectanceSpectra = [SampleReflectanceSpectra; TempReflectances(ctr,:)]	;
	end
	% Find the MCDM for just that sample
	[MCDM, ~] = MeanColourDifferenceFromTheMean( ...
											Wavelengths, ...
											SampleReflectanceSpectra, ...
											IllumObs);
	% Add the MCDM for that sample to the running list											
	MCDMs = [MCDMs, MCDM]	;
	
	if PlotOutliers		% Plot reflectance spectra, if desired, when the MCDM is high
		if MCDM > MCDMThreshold 
			disp(['Sample ',num2str(ctr),': MCDM is ',num2str(MCDM)])
		    fflush(stdout);

			FontName = 'Times New Roman'		;
			LabelFontWeight = 'normal'	;
			figure
			for filectr = 1:NumberOfFiles
				plot(Wavelengths, SampleReflectanceSpectra(filectr,:), 'k-')
				hold on
				text(Wavelengths(end),SampleReflectanceSpectra(filectr,end),num2str(filectr),...
					'horizontalalignment', 'left')
				hold on
			end

			set(gca, 'ylim', [0 1])
			set(gca, 'ytick', 0:0.2:1)
			set(gca, 'yticklabel', {'0', '20', '40', '60', '80', '100'}, ...
					 'fontname', FontName, ...
					 'fontweight', LabelFontWeight,...
					 'fontsize', AxisFontSizeInPoints)	
			xlabel('Wavelength (nm)', ...
					 'fontname', FontName, ...
					 'fontweight', LabelFontWeight,...
					 'fontsize', LabelFontSizeInPoints)
			ylabel('Reflectance (%)',...
					 'fontname', FontName, ...
					 'fontweight', LabelFontWeight,...
					 'fontsize', LabelFontSizeInPoints)
			figname = [NameKey, 'Sample',num2str(ctr)]	;
			set(gcf, 'Name', figname)
			figname = [OutputDirectory,'/',figname]	;
			print(gcf, [figname,'.eps'], '-depsc')		;
			print(gcf, [figname,'.png'], '-dpng')		;
			print(gcf, [figname,'.jpg'], '-djpg')		;
			print(gcf, [figname,'.pdf'], '-dpdf')		;
			fflush(stdout);
		end
	end
end

disp([num2str(NumberOfSamples),' colour samples, each measured ',num2str(NumberOfFiles),' times']);
disp(['Mean MCDM:   ',num2str(mean(MCDMs))])	;
disp(['Median MCDM: ',num2str(median(MCDMs))])	;
disp(['Min MCDM:    ',num2str(min(MCDMs))])		;
disp(['Max MCDM:    ',num2str(max(MCDMs))])		;

% Plot sequence of MCDMs, if desired
if DrawSequence
	figure
	plot([1:length(MCDMs)], MCDMs, 'k*')
	set(gca, 'ylim', [0, 1.2*max(MCDMs)])		
	figname = [NameKey, 'MCDMSequence']			;
	set(gcf, 'Name', figname)
	figname = [OutputDirectory,'/',figname]		;
	print(gcf, [figname,'.eps'], '-depsc')		;
	print(gcf, [figname,'.png'], '-dpng')		;
	print(gcf, [figname,'.jpg'], '-djpg')		;
	print(gcf, [figname,'.pdf'], '-dpdf')		;
end

% Print out histogram of MCDMs, if desired
if DrawHistogram
	EdgeVector = [0:0.25:2]					;
	MCDMCounts = histc(MCDMs, EdgeVector)	;
	figure
	stairs(EdgeVector, MCDMCounts, 'k')		;
	set(gca, 'xlim', [0 2])					;
	set(gca, 'xtick', 0:0.25:2)				;
	FontName = 'Times New Roman'			;
	LabelFontSizeInPoints = 14				;
	LabelFontWeight = 'normal'				;
	xlabel('MCDM', ...
			 'fontname', FontName, ...
			 'fontweight', LabelFontWeight,...
			 'fontsize', LabelFontSizeInPoints)
	ylabel('Number of Occurrences',...
			 'fontname', FontName, ...
			 'fontweight', LabelFontWeight,...
			 'fontsize', LabelFontSizeInPoints)
	figname = [NameKey,'MCDMHistogram']	;
	set(gcf, 'Name', figname)
	figname = [OutputDirectory,'/',figname]	;
	print(gcf, [figname,'.eps'], '-depsc')		;
	print(gcf, [figname,'.png'], '-dpng')		;
	print(gcf, [figname,'.jpg'], '-djpg')		;
	print(gcf, [figname,'.pdf'], '-dpdf')		;
end