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