Crosscheck - skinoptics.colors and colour-science

Victor Lima
victorporto@ifsc.usp.br
victor.lima@ufscar.br
victorportog.github.io

Release date:
01 April 2025
Last modification:
01 April 2025

This notebook compares the skinoptics.colors module with the colour-science package (https://www.colour-science.org/) for numerical calculations of CIE XYZ, sRGB and CIE L*a*b* color space coordinates from human skin reflectance spectra.

The skin reflectance spectra used here are from the publicly available database by Xiao et al. [X*16].

All calculations were performed assuming the CIE standard illuminant D65 and the CIE 10 degree standard observer (except for the sRGB coordinates, for which the CIE standard illuminant D65 and the CIE 2 degree standard observer were assumed, as default).

CIE XYZ coordinates were computed by numerical integration of the reflectance spectrum in the region between 360 and 830 nm. A cubic spline interpolation was carried out from 360 to 740 nm and a constant extrapolation was carried out from 740 to 830 nm. Then, sRGB and CIE L*a*b* coordinates were obtained from CIE XYZ coordinates.

This file is part of SkinOptics documentation.

References:

[X*16] Xiao, Zhu, Li, Connah, Yates & Wuerger 2016.
Improved method for skin reflectance reconstruction from camera images.
https://doi.org/10.1364/OE.24.014934

[F23] Feranández-Barral 2023.
Best practices for colour blind friendly publications & descriptions.
https://pos.sissa.it/guidelines.pdf

import time
import numpy as np
import matplotlib
from matplotlib import rcParams, style
import matplotlib.pyplot as plt
import pandas as pd

print('numpy version:', np.__version__)
print('matplotlib version:', matplotlib.__version__)
print('pandas version:', pd.__version__)

numpy version: 1.26.4
matplotlib version: 3.9.0
pandas version: 2.2.2

rcParams.update({'font.size': 12})
rcParams.update({'axes.labelsize': 12})
rcParams.update({'axes.titlesize': 12})
rcParams.update({'xtick.labelsize': 12})
rcParams.update({'ytick.labelsize': 12})
rcParams.update({'legend.fontsize': 10})
rcParams.update({'axes.grid': True})
plt.rcParams["figure.figsize"] = (6, 4)

color_cycle =  ['#CC78BC', '#029E73', '#56B4E9']

A bigger font size (12.) than plt deafult (10.) was set and some colors from 'seaborn-colorbind' palette were chosen.
This color palette may be more inclusive for people with color vision deficiencies [F23].

Importing tools from skinoptics.colors and colour-science

import skinoptics
import colour
from skinoptics.colors import *
from colour import *

print('skinoptics version:', skinoptics.__version__)
print('colour-science version:', colour.__version__)

skinoptics version: 0.0.1
colour-science version: 0.4.6

Getting data from the database by Xiao et al. 2016

all_lambda = np.arange(360, 740 + 10, 10)
path = os.getcwd()
parent_path = os.path.dirname(path)
reflectance = np.array(pd.read_excel(parent_path + '\skinoptics\datasets\spectra\Xiao2016\skindatabaseSpectra\skin spectra data.xlsx'))[:,1:]
n = len(reflectance)
print(n)

4392

Visualizing the cubic spline interpolation and the constant extrapolation

plt.scatter(all_lambda, reflectance[0,:]*100, color = 'k', marker = '.', label = 'data')
plt.plot(all_lambda, interp1d(all_lambda, reflectance[0,:]*100, kind = 'cubic',
                              bounds_error = False,
                              fill_value = (reflectance[0,0]*100, reflectance[0,-1]*100))(all_lambda),
         color = 'b', lw = 2., label = 'interpolation')
x = np.arange(740, 830 + 1, 1)
plt.plot(x, interp1d(all_lambda, reflectance[0,:]*100, kind = 'cubic',
                     bounds_error = False, fill_value = (reflectance[0,0]*100, reflectance[0,-1]*100))(x),
         color = 'b', ls = '--', lw = 2., label = 'extrapolation')
plt.xlabel('wavelength [nm]')
plt.ylabel('reflectance [%]')
plt.title('spectrum #1')
plt.legend(loc = 'upper left')
plt.xlim(360, 830)
plt.ylim(0, 80)
plt.show()

Calculation using skinoptics.colors

X_sc, Y_sc, Z_sc, R_sc, G_sc, B_sc, L_sc, a_sc, b_sc = np.zeros([9, n])
t0 = time.time()
for i in range(n):
    X_sc[i], Y_sc[i], Z_sc[i] = XYZ_from_spectrum(all_lambda, reflectance[i,:]*100)
    R_sc[i], G_sc[i], B_sc[i] = sRGB_from_spectrum(all_lambda, reflectance[i,:]*100)
    L_sc[i], a_sc[i], b_sc[i] = Lab_from_spectrum(all_lambda, reflectance[i,:]*100)
print('execution time: {} seconds (skinoptics.colors)'.format(round(time.time() - t0, 2)))

execution time: 24.65 seconds (skinoptics.colors)

Calculation using colour-science

X_cs, Y_cs, Z_cs, R_cs, G_cs, B_cs, L_cs, a_cs, b_cs = np.zeros([9, n])
t0 = time.time()
for i in range(n):
    sd = SpectralDistribution(dict_from_arrays(array_keys = all_lambda, array_values = reflectance[i, :]), name = 'Sample')
    cmfs2 = MSDS_CMFS['CIE 1931 2 Degree Standard Observer']
    cmfs10 = MSDS_CMFS['CIE 1964 10 Degree Standard Observer']
    illuminant = SDS_ILLUMINANTS['D65']
    X_cs[i], Y_cs[i], Z_cs[i] = sd_to_XYZ(sd, cmfs10, illuminant)/100
    R_cs[i], G_cs[i], B_cs[i] = XYZ_to_sRGB(sd_to_XYZ(sd, cmfs2, illuminant)/100)
    L_cs[i], a_cs[i], b_cs[i] = XYZ_to_Lab(sd_to_XYZ(sd, cmfs10, illuminant)/100, 
                                           illuminant = CCS_ILLUMINANTS["CIE 1964 10 Degree Standard Observer"]["D65"])
print('execution time: {} seconds (colour-science)'.format(round(time.time() - t0, 2)))

execution time: 11.9 seconds (colour-science)

Comparing CIE XYZ coordinates

rel_diff_X_sc_cs = np.abs(X_sc/X_cs - 1)
rel_diff_Y_sc_cs = np.abs(Y_sc/Y_cs - 1)
rel_diff_Z_sc_cs = np.abs(Z_sc/Z_cs - 1)

g, axs = plt.subplots(3, 3,
                        gridspec_kw = {'wspace': 0.35, 'hspace': 0.4,
                                       'width_ratios': [2.5, 2.5, 1]}, figsize = (14, 12))

axs[0,0].scatter(X_sc, Y_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.', label = 'skinoptics.colors')
axs[0,0].scatter(X_cs, Y_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x',label = 'colour-science')
axs[0,0].set_xlabel('X')
axs[0,0].set_ylabel('Y')
axs[0,0].set_title('(a)', loc = 'left')
axs[0,0].legend(loc = 'upper left')

axs[1,0].scatter(X_sc, Z_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.')
axs[1,0].scatter(X_cs, Z_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x')
axs[1,0].set_xlabel('X')
axs[1,0].set_ylabel('Z')
axs[1,0].set_title('(b)', loc = 'left')

axs[2,0].scatter(Y_sc, Z_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.')
axs[2,0].scatter(Y_cs, Z_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x')
axs[2,0].set_xlabel('Y')
axs[2,0].set_ylabel('Z')
axs[2,0].set_title('(c)', loc = 'left')

axs[0,1].scatter(np.arange(n), rel_diff_X_sc_cs*1E5, color = color_cycle[0], marker = '.')
axs[0,1].set_xlabel('spectrum index [-]')
axs[0,1].set_ylabel('relative diff. in X [10$^{-5}$]')
axs[0,1].set_title('(d)', loc = 'left')
axs[0,1].set_xlim(0, n)
axs[0,1].set_ylim(0, rel_diff_X_sc_cs.max()*1E5,)

axs[1,1].scatter(np.arange(n), rel_diff_Y_sc_cs*1E5, color = color_cycle[1], marker = '.')
axs[1,1].set_xlabel('spectrum index [-]')
axs[1,1].set_ylabel('relative diff. in Y [10$^{-5}$]')
axs[1,1].set_title('(e)', loc = 'left')
axs[1,1].set_xlim(0, n)
axs[1,1].set_ylim(0, rel_diff_Y_sc_cs.max()*1E5,)

axs[2,1].scatter(np.arange(n), rel_diff_Z_sc_cs*1E5, color = color_cycle[2], marker = '.')
axs[2,1].set_xlabel('spectrum index [-]')
axs[2,1].set_ylabel('relative diff. in Z [10$^{-5}$]')
axs[2,1].set_title('(f)', loc = 'left')
axs[2,1].set_xlim(0, n)
axs[2,1].set_ylim(0, rel_diff_Z_sc_cs.max()*1E5,)

axs[0,2].hist(rel_diff_X_sc_cs*1E5, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[0], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[0,2].hlines(rel_diff_X_sc_cs.mean()*1E5, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[0,2].set_xlabel('count [-]')
axs[0,2].set_ylabel('relative diff. in X [10$^{-5}$]')
axs[0,2].set_title('(g)', loc = 'left')
axs[0,2].set_xlim(0, 400)
axs[0,2].set_ylim(0, rel_diff_X_sc_cs.max()*1E5,)

axs[1,2].hist(rel_diff_Y_sc_cs*1E5, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[1], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[1,2].hlines(rel_diff_Y_sc_cs.mean()*1E5, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[1,2].set_xlabel('count [-]')
axs[1,2].set_ylabel('relative diff. in Y [10$^{-5}$]')
axs[1,2].set_title('(h)', loc = 'left')
axs[1,2].set_xlim(0, 400)
axs[1,2].set_ylim(0, rel_diff_Y_sc_cs.max()*1E5,)

axs[2,2].hist(rel_diff_Z_sc_cs*1E5, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[2], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[2,2].hlines(rel_diff_Z_sc_cs.mean()*1E5, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[2,2].set_xlabel('count [-]')
axs[2,2].set_ylabel('relative diff. in Z [10$^{-5}$]')
axs[2,2].set_title('(i)', loc = 'left')
axs[2,2].set_xlim(0, 400)
axs[2,2].set_ylim(0, rel_diff_Z_sc_cs.max()*1E5,)

plt.show()

pd.DataFrame(np.array([[rel_diff_X_sc_cs.mean()*1E5],
                       [rel_diff_Y_sc_cs.mean()*1E5],
                       [rel_diff_Z_sc_cs.mean()*1E5]]),
             index = ['mean relative diff. in X [10$^{-5}$]',
                      'mean relative diff. in Y [10$^{-5}$]',
                      'mean relative diff. in Z [10$^{-5}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
mean relative diff. in X [10$^{-5}$]	0.239637
mean relative diff. in Y [10$^{-5}$]	0.604465
mean relative diff. in Z [10$^{-5}$]	1.636001

pd.DataFrame(np.array([[rel_diff_X_sc_cs.max()*1E5],
                       [rel_diff_Y_sc_cs.max()*1E5],
                       [rel_diff_Z_sc_cs.max()*1E5]]),
             index = ['max relative diff. in X [10$^{-5}$]',
                      'max relative diff. in Y [10$^{-5}$]',
                      'max relative diff. in Z [10$^{-5}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
max relative diff. in X [10$^{-5}$]	2.137114
max relative diff. in Y [10$^{-5}$]	1.713784
max relative diff. in Z [10$^{-5}$]	7.118888

pd.DataFrame(np.array([[rel_diff_X_sc_cs.min()*1E5],
                       [rel_diff_Y_sc_cs.min()*1E5],
                       [rel_diff_Z_sc_cs.min()*1E5]]),
             index = ['min relative diff. in X [10$^{-5}$]',
                      'min relative diff. in Y [10$^{-5}$]',
                      'min relative diff. in Z [10$^{-5}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
min relative diff. in X [10$^{-5}$]	0.000164
min relative diff. in Y [10$^{-5}$]	0.002226
min relative diff. in Z [10$^{-5}$]	0.000510

Comparing sRGB coordinates

rel_diff_R_sc_cs = np.abs(R_sc/R_cs - 1)
rel_diff_G_sc_cs = np.abs(G_sc/G_cs - 1)
rel_diff_B_sc_cs = np.abs(B_sc/B_cs - 1)

fig, axs = plt.subplots(3, 3,
                        gridspec_kw = {'wspace': 0.35, 'hspace': 0.4,
                                       'width_ratios': [2.5, 2.5, 1]}, figsize = (14, 12))

axs[0,0].scatter(R_sc, G_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.', label = 'skinoptics.colors')
axs[0,0].scatter(R_cs, G_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x', label = 'colour-science')
axs[0,0].set_xlabel('R')
axs[0,0].set_ylabel('G')
axs[0,0].set_title('(a)', loc = 'left')
axs[0,0].legend(loc = 'upper left')

axs[1,0].scatter(R_sc, B_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.')
axs[1,0].scatter(R_cs, B_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x')
axs[1,0].set_xlabel('R')
axs[1,0].set_ylabel('B')
axs[1,0].set_title('(b)', loc = 'left')

axs[2,0].scatter(G_sc, B_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.')
axs[2,0].scatter(G_cs, B_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x')
axs[2,0].set_xlabel('G')
axs[2,0].set_ylabel('B')
axs[2,0].set_title('(c)', loc = 'left')

axs[0,1].scatter(np.arange(n), rel_diff_R_sc_cs*1E5, color = color_cycle[0], marker = '.')
axs[0,1].set_xlabel('spectrum index [-]')
axs[0,1].set_ylabel('relative diff. in R [10$^{-5}$]')
axs[0,1].set_title('(d)', loc = 'left')
axs[0,1].set_xlim(0, n)
axs[0,1].set_ylim(0, rel_diff_R_sc_cs.max()*1E5)

axs[1,1].scatter(np.arange(n), rel_diff_G_sc_cs*1E5, color = color_cycle[1], marker = '.')
axs[1,1].set_xlabel('spectrum index [-]')
axs[1,1].set_ylabel('relative diff. in G [10$^{-5}$]')
axs[1,1].set_title('(e)', loc = 'left')
axs[1,1].set_xlim(0, n)
axs[1,1].set_ylim(0, rel_diff_G_sc_cs.max()*1E5)

axs[2,1].scatter(np.arange(n), rel_diff_B_sc_cs*1E5, color = color_cycle[2], marker = '.')
axs[2,1].set_xlabel('spectrum index [-]')
axs[2,1].set_ylabel('relative diff. in B [10$^{-5}$]')
axs[2,1].set_title('(f)', loc = 'left')
axs[2,1].set_xlim(0, n)
axs[2,1].set_ylim(0, rel_diff_B_sc_cs.max()*1E5)

axs[0,2].hist(rel_diff_R_sc_cs*1E5, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[0], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[0,2].hlines(rel_diff_R_sc_cs.mean()*1E5, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[0,2].set_xlabel('count [-]')
axs[0,2].set_ylabel('relative diff. in R [10$^{-5}$]')
axs[0,2].set_title('(g)', loc = 'left')
axs[0,2].set_xlim(0, 400)
axs[0,2].set_ylim(0, rel_diff_R_sc_cs.max()*1E5)

axs[1,2].hist(rel_diff_G_sc_cs*1E5, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[1], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[1,2].hlines(rel_diff_G_sc_cs.mean()*1E5, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[1,2].set_xlabel('count [-]')
axs[1,2].set_ylabel('relative diff. in G [10$^{-5}$]')
axs[1,2].set_title('(h)', loc = 'left')
axs[1,2].set_xlim(0, 400)
axs[1,2].set_ylim(0, rel_diff_G_sc_cs.max()*1E5)

axs[2,2].hist(rel_diff_B_sc_cs*1E5, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[2], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[2,2].hlines(rel_diff_B_sc_cs.mean()*1E5, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[2,2].set_xlabel('count [-]')
axs[2,2].set_ylabel('relative diff. in B [10$^{-5}$]')
axs[2,2].set_title('(i)', loc = 'left')
axs[2,2].set_xlim(0, 400)
axs[2,2].set_ylim(0, rel_diff_B_sc_cs.max()*1E5)

plt.show()

pd.DataFrame(np.array([[rel_diff_R_sc_cs.mean()*1E5],
                       [rel_diff_G_sc_cs.mean()*1E5],
                       [rel_diff_B_sc_cs.mean()*1E5]]),
             index = ['mean relative diff. in R [10$^{-5}$]',
                      'mean relative diff. in G [10$^{-5}$]',
                      'mean relative diff. in B [10$^{-5}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
mean relative diff. in R [10$^{-5}$]	1.382633
mean relative diff. in G [10$^{-5}$]	1.305374
mean relative diff. in B [10$^{-5}$]	1.132943

pd.DataFrame(np.array([[rel_diff_R_sc_cs.max()*1E5],
                       [rel_diff_G_sc_cs.max()*1E5],
                       [rel_diff_B_sc_cs.max()*1E5]]),
             index = ['max relative diff. in R [10$^{-5}$]',
                      'max relative diff. in G [10$^{-5}$]',
                      'max relative diff. in B [10$^{-5}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
max relative diff. in R [10$^{-5}$]	2.978674
max relative diff. in G [10$^{-5}$]	4.741567
max relative diff. in B [10$^{-5}$]	4.778786

pd.DataFrame(np.array([[rel_diff_R_sc_cs.min()*1E5],
                       [rel_diff_G_sc_cs.min()*1E5],
                       [rel_diff_B_sc_cs.min()*1E5]]),
             index = ['min relative diff. in R [10$^{-5}$]',
                      'min relative diff. in G [10$^{-5}$]',
                      'min relative diff. in B [10$^{-5}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
min relative diff. in R [10$^{-5}$]	0.055746
min relative diff. in G [10$^{-5}$]	0.027578
min relative diff. in B [10$^{-5}$]	0.000062

Comparing CIE L*a*b* coordinates

rel_diff_L_sc_cs = np.abs(L_sc/L_cs - 1)
rel_diff_a_sc_cs = np.abs(a_sc/a_cs - 1)
rel_diff_b_sc_cs = np.abs(b_sc/b_cs - 1)

fig, axs = plt.subplots(3, 3,
                        gridspec_kw = {'wspace': 0.35, 'hspace': 0.4, 'width_ratios': [2.5, 2.5, 1]}, figsize = (14, 12))

axs[0,0].scatter(L_sc, a_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.', label = 'skinoptics.colors')
axs[0,0].scatter(L_cs, a_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x', label = 'colour-science')
axs[0,0].set_xlabel('L*')
axs[0,0].set_ylabel('a*')
axs[0,0].set_title('(a)', loc = 'left')
axs[0,0].legend(loc = 'upper left')

axs[1,0].scatter(L_sc, b_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.')
axs[1,0].scatter(L_cs, b_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x')
axs[1,0].set_xlabel('L*')
axs[1,0].set_ylabel('b*')
axs[1,0].set_title('(b)', loc = 'left')

axs[2,0].scatter(a_sc, b_sc, color = np.dstack((R_sc, G_sc, B_sc))[0,:,:], marker = '.')
axs[2,0].scatter(a_cs, b_cs, color = np.dstack((R_cs, G_cs, B_cs))[0,:,:], marker = 'x')
axs[2,0].set_xlabel('a*')
axs[2,0].set_ylabel('b*')
axs[2,0].set_title('(c)', loc = 'left')

axs[0,1].scatter(np.arange(n), rel_diff_L_sc_cs*1E4, color = color_cycle[0], marker = '.')
axs[0,1].set_xlabel('spectrum index [-]')
axs[0,1].set_ylabel('relative diff. in L* [10$^{-4}$]')
axs[0,1].set_title('(d)', loc = 'left')
axs[0,1].set_xlim(0, n)
axs[0,1].set_ylim(0, rel_diff_L_sc_cs.max()*1E4)

axs[1,1].scatter(np.arange(n), rel_diff_a_sc_cs*1E4, color = color_cycle[1], marker = '.')
axs[1,1].set_xlabel('spectrum index [-]')
axs[1,1].set_ylabel('relative diff. in a* [10$^{-4}$]')
axs[1,1].set_title('(e)', loc = 'left')
axs[1,1].set_xlim(0, n)
axs[1,1].set_ylim(0, rel_diff_a_sc_cs.max()*1E4)

axs[2,1].scatter(np.arange(n), rel_diff_b_sc_cs*1E4, color = color_cycle[2], marker = '.')
axs[2,1].set_xlabel('spectrum index [-]')
axs[2,1].set_ylabel('relative diff. in b* [10$^{-4}$]')
axs[2,1].set_title('(f)', loc = 'left')
axs[2,1].set_xlim(0, n)
axs[2,1].set_ylim(0, rel_diff_b_sc_cs.max()*1E4)

axs[0,2].hist(rel_diff_L_sc_cs*1E4, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[0], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[0,2].hlines(rel_diff_L_sc_cs.mean()*1E4, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[0,2].set_xlabel('count [-]')
axs[0,2].set_ylabel('relative diff. in L* [10$^{-4}$]')
axs[0,2].set_title('(g)', loc = 'left')
axs[0,2].set_xlim(0, 400)
axs[0,2].set_ylim(0, rel_diff_L_sc_cs.max()*1E4)

axs[1,2].hist(rel_diff_a_sc_cs*1E4, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[1], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[1,2].hlines(rel_diff_a_sc_cs.mean()*1E4, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[1,2].set_xlabel('count [-]')
axs[1,2].set_ylabel('relative diff. in a* [10$^{-4}$]')
axs[1,2].set_title('(h)', loc = 'left')
axs[1,2].set_xlim(0, 400)
axs[1,2].set_ylim(0, rel_diff_a_sc_cs.max()*1E4)

axs[2,2].hist(rel_diff_b_sc_cs*1E4, bins = 'fd', orientation = 'horizontal', lw = 1.5,
              color = color_cycle[2], alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[2,2].hlines(rel_diff_b_sc_cs.mean()*1E4, 0, 400, color = 'k', ls = '--', lw = 2.)
axs[2,2].set_xlabel('count [-]')
axs[2,2].set_ylabel('relative diff. in b* [10$^{-4}$]')
axs[2,2].set_title('(i)', loc = 'left')
axs[2,2].set_xlim(0, 400)
axs[2,2].set_ylim(0, rel_diff_b_sc_cs.max()*1E4)

plt.show()

pd.DataFrame(np.array([[rel_diff_L_sc_cs.mean()*1E4],
                       [rel_diff_a_sc_cs.mean()*1E4],
                       [rel_diff_b_sc_cs.mean()*1E4]]),
             index = ['mean relative diff. in L* [10$^{-4}$]',
                      'mean relative diff. in a* [10$^{-4}$]',
                      'mean relative diff. in b* [10$^{-4}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
mean relative diff. in L* [10$^{-4}$]	0.025814
mean relative diff. in a* [10$^{-4}$]	0.955426
mean relative diff. in b* [10$^{-4}$]	0.660821

pd.DataFrame(np.array([[rel_diff_L_sc_cs.max()*1E4],
                       [rel_diff_a_sc_cs.max()*1E4],
                       [rel_diff_b_sc_cs.max()*1E4]]),
             index = ['max relative diff. in L* [10$^{-4}$]',
                      'max relative diff. in a* [10$^{-4}$]',
                      'max relative diff. in b* [10$^{-4}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
max relative diff. in L* [10$^{-4}$]	0.093207
max relative diff. in a* [10$^{-4}$]	6.104785
max relative diff. in b* [10$^{-4}$]	3.197917

pd.DataFrame(np.array([[rel_diff_L_sc_cs.min()*1E4],
                       [rel_diff_a_sc_cs.min()*1E4],
                       [rel_diff_b_sc_cs.min()*1E4]]),
             index = ['min relative diff. in L* [10$^{-4}$]',
                      'min relative diff. in a* [10$^{-4}$]',
                      'min relative diff. in b* [10$^{-4}$]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
min relative diff. in L* [10$^{-4}$]	0.000092
min relative diff. in a* [10$^{-4}$]	0.001069
min relative diff. in b* [10$^{-4}$]	0.000639

Color difference

fig, axs = plt.subplots(1, 2,
                        gridspec_kw = {'wspace': 0.4, 'hspace': 0.3,
                                       'width_ratios': [2, 1]}, figsize = (9, 4))

Delta_E_sc_cs = Delta_E(L_sc, a_sc, b_sc, L_cs, a_cs, b_cs)
axs[0].scatter(np.arange(n), Delta_E_sc_cs, color = 'gray', marker = '.', label = 'skinoptics.colors & colour-science')
axs[0].set_xlabel('spectrum index [-]')
axs[0].set_ylabel('$\Delta$E*')
axs[0].legend(loc = 'upper left')
axs[0].set_xlim(0, n)
axs[0].set_ylim(0, Delta_E_sc_cs.max())
axs[0].set_title('(a)', loc = 'left')

axs[1].hist(Delta_E_sc_cs, bins = 'fd', orientation = 'horizontal', lw = 1.5,
            color = 'gray', alpha = 1., edgecolor = 'k', histtype = 'stepfilled')
axs[1].hlines(Delta_E_sc_cs.mean(), 0, 400, color = 'k', ls = '--', lw = 2.)
axs[1].set_xlabel('count [-]')
axs[1].set_ylabel('$\Delta$E*')
axs[1].set_xlim(0, 400)
axs[1].set_ylim(0, Delta_E_sc_cs.max())
axs[1].set_title('(b)', loc = 'left')

plt.show()

pd.DataFrame(np.array([[Delta_E_sc_cs.mean()],
                       [Delta_E_sc_cs.max()],
                       [Delta_E_sc_cs.min()]]),
             index = ['mean $\\Delta$E* [-]', 'max $\\Delta$E* [-]', 'min $\\Delta$E* [-]'],
             columns = ['skinoptics.colors & colour-science'])

	skinoptics.colors & colour-science
mean $\Delta$E* [-]	0.001345
max $\Delta$E* [-]	0.003466
min $\Delta$E* [-]	0.000527

Extra: ideal grayscale diffuse reflectance standards (DRSs)

\begin{equation} R_d(\lambda) = \alpha \end{equation}

$\alpha \in [0, 1]$, $\lambda \in [360, 830]$ nm

all_alpha = np.arange(100 + 1)
all_lambda = np.arange(360, 830 + 1, 1)
Rd = np.zeros(len(all_lambda))

X_DRS_sc, Y_DRS_sc, Z_DRS_sc, R_DRS_sc, G_DRS_sc, B_DRS_sc, L_DRS_sc, a_DRS_sc, b_DRS_sc = np.ones((9, 100 + 1))
for i in range(100 + 1):
    Rd = i*np.ones(len(all_lambda))
    X_DRS_sc[i], Y_DRS_sc[i], Z_DRS_sc[i] = XYZ_from_spectrum(all_lambda, Rd)
    R_DRS_sc[i], G_DRS_sc[i], B_DRS_sc[i] = sRGB_from_spectrum(all_lambda, Rd)
    L_DRS_sc[i], a_DRS_sc[i], b_DRS_sc[i] = Lab_from_spectrum(all_lambda, Rd)

X_DRS_cs, Y_DRS_cs, Z_DRS_cs, R_DRS_cs, G_DRS_cs, B_DRS_cs, L_DRS_cs, a_DRS_cs, b_DRS_cs = np.zeros([9, 100 + 1])
for i in range(100 + 1):
    Rd = i*0.01*np.ones(len(all_lambda))
    sd = SpectralDistribution(dict_from_arrays(array_keys = all_lambda, array_values = Rd), name = 'Sample')
    cmfs2 = MSDS_CMFS['CIE 1931 2 Degree Standard Observer']
    cmfs10 = MSDS_CMFS['CIE 1964 10 Degree Standard Observer']
    illuminant = SDS_ILLUMINANTS['D65']
    X_DRS_cs[i], Y_DRS_cs[i], Z_DRS_cs[i] = sd_to_XYZ(sd, cmfs10, illuminant)/100
    R_DRS_cs[i], G_DRS_cs[i], B_DRS_cs[i] = XYZ_to_sRGB(sd_to_XYZ(sd, cmfs2, illuminant)/100)
    L_DRS_cs[i], a_DRS_cs[i], b_DRS_cs[i] = XYZ_to_Lab(sd_to_XYZ(sd, cmfs10, illuminant)/100, 
                                                       illuminant = CCS_ILLUMINANTS["CIE 1964 10 Degree Standard Observer"]["D65"])

fig, axs = plt.subplots(3, 3,
                        gridspec_kw = {'wspace': 0.4, 'hspace': 0.45}, figsize = (12, 12))

axs[0,0].plot(all_alpha, X_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[0,0].plot(all_alpha, X_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[0,0].set_xlabel('$\\alpha$ [%]')
axs[0,0].set_ylabel('X (D65, 10 degree)')
axs[0,0].set_xlim(0, 100)
axs[0,0].set_ylim(0, 1.2)
axs[0,0].legend(loc = 'upper left')
axs[0,0].set_title('(a)', loc = 'left')

axs[1,0].plot(all_alpha, Y_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[1,0].plot(all_alpha, Y_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[1,0].set_xlabel('$\\alpha$ [%]')
axs[1,0].set_ylabel('Y (D65, 10 degree)')
axs[1,0].set_xlim(0, 100)
axs[1,0].set_ylim(0, 1.2)
axs[1,0].legend(loc = 'upper left')
axs[1,0].set_title('(b)', loc = 'left')

axs[2,0].plot(all_alpha, Z_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[2,0].plot(all_alpha, Z_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[2,0].set_xlabel('$\\alpha$ [%]')
axs[2,0].set_ylabel('Z (D65, 10 degree)')
axs[2,0].set_xlim(0, 100)
axs[2,0].set_ylim(0, 1.2)
axs[2,0].legend(loc = 'upper left')
axs[2,0].set_title('(c)', loc = 'left')

axs[0,1].plot(all_alpha, nonlinear_corr_sRGB(all_alpha/100), 'k', lw = 2., label = 'analytical R = $\\gamma$($\\alpha$)')
axs[0,1].plot(all_alpha, R_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[0,1].plot(all_alpha, R_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[0,1].set_xlabel('$\\alpha$ [%]')
axs[0,1].set_ylabel('R (D65, 2 degree)')
axs[0,1].set_xlim(0, 100)
axs[0,1].set_ylim(0, 1.)
axs[0,1].legend(loc = 'lower right')
axs[0,1].set_title('(d)', loc = 'left')

axs[1,1].plot(all_alpha, nonlinear_corr_sRGB(all_alpha/100), 'k', lw = 2., label = 'analytical G = $\\gamma$($\\alpha$)')
axs[1,1].plot(all_alpha, G_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[1,1].plot(all_alpha, G_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[1,1].set_xlabel('$\\alpha$ [%]')
axs[1,1].set_ylabel('G (D65, 2 degree)')
axs[1,1].set_xlim(0, 100)
axs[1,1].set_ylim(0, 1.)
axs[1,1].legend(loc = 'lower right')
axs[1,1].set_title('(e)', loc = 'left')

axs[2,1].plot(all_alpha, nonlinear_corr_sRGB(all_alpha/100), 'k', lw = 2., label = 'analytical B = $\\gamma$($\\alpha$)')
axs[2,1].plot(all_alpha, B_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[2,1].plot(all_alpha, B_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[2,1].set_xlabel('$\\alpha$ [%]')
axs[2,1].set_ylabel('B (D65, 2 degree)')
axs[2,1].set_xlim(0, 100)
axs[2,1].set_ylim(0, 1.)
axs[2,1].legend(loc = 'lower right')
axs[2,1].set_title('(f)', loc = 'left')

axs[0,2].plot(all_alpha, 116*f_Lab_from_XYZ(all_alpha/100) - 16, 'k', lw = 2., label = 'analytical\nL$^∗$ = 116 f($\\alpha$) − 16')
axs[0,2].plot(all_alpha, L_DRS_sc, '--r', lw = 2., label = 'skinoptics.colors')
axs[0,2].plot(all_alpha, L_DRS_cs, ':b', lw = 2., label = 'colour-science')
axs[0,2].set_xlabel('$\\alpha$ [%]')
axs[0,2].set_ylabel('L* (D65, 10 degree)')
axs[0,2].set_xlim(0, 100)
axs[0,2].set_ylim(0, 100)
axs[0,2].legend(loc = 'lower right')
axs[0,2].set_title('(g)', loc = 'left')

axs[1,2].plot(all_alpha, np.zeros(len(all_alpha)), 'k', lw = 2., label = 'analytical a$^*$ = 0')
axs[1,2].plot(all_alpha, a_DRS_sc*1E3, '--r', lw = 2., label = 'skinoptics.colors')
axs[1,2].plot(all_alpha, a_DRS_cs*1E3, ':b', lw = 2., label = 'colour-science')
axs[1,2].set_xlabel('$\\alpha$ [%]')
axs[1,2].set_ylabel('a* (D65, 10 degree) [10$^{-3}$]')  
axs[1,2].set_xlim(0, 100)
axs[1,2].set_ylim(-3, 3)
axs[1,2].legend(loc = 'lower right')
axs[1,2].set_title('(h)', loc = 'left')

axs[2,2].plot(all_alpha, np.zeros(len(all_alpha)), 'k', lw = 2., label = 'analytical b$^*$ = 0')
axs[2,2].plot(all_alpha, b_DRS_sc*1E3, '--r', lw = 2., label = 'skinoptics.colors')
axs[2,2].plot(all_alpha, b_DRS_cs*1E3, ':b', lw = 2., label = 'colour-science')
axs[2,2].set_xlabel('$\\alpha$ [%]')
axs[2,2].set_ylabel('b* (D65, 10 degree) [10$^{-3}$]') 
axs[2,2].set_xlim(0, 100)
axs[2,2].set_ylim(-3, 3)
axs[2,2].legend(loc = 'lower right')
axs[2,2].set_title('(i)', loc = 'left')

plt.show()

plt.plot(all_alpha, np.abs(L_DRS_sc - (116*f_Lab_from_XYZ(all_alpha/100) - 16))*1E14,
            color = color_cycle[0], marker = '.', lw = 2., label = 'L*')
plt.plot(all_alpha, np.abs(a_DRS_sc - np.zeros(len(all_alpha)))*1E14,
            color = color_cycle[1], marker = '.', lw = 2., label = 'a*')
plt.plot(all_alpha, np.abs(b_DRS_sc - np.zeros(len(all_alpha)))*1E14,
            color = color_cycle[2], marker = '.', lw = 2., label = 'b*')
plt.xlabel('$\\alpha$ [%]')
plt.ylabel('absolute difference [10$^{-14}$]') 
plt.title('skinoptics.colors & analytical')
plt.xlim(0, 100)
plt.ylim(0, 16)
plt.legend(loc = 'upper right')
plt.show()