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 .
% 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