summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPaweł Redman <pawel.redman@gmail.com>2020-07-22 16:49:18 +0200
committerPaweł Redman <pawel.redman@gmail.com>2020-07-22 16:49:18 +0200
commit162aeb511f5bafc67b16fb920f9f8d037d4ab7ca (patch)
tree5fe7b45b087b9483e37d761399c43ccd52d773bb
parentdb0040dbe91729b9bd78c40b4be75068b1e6968a (diff)
Add a simple script for computing delta E's.
-rw-r--r--delta_E.py103
1 files changed, 103 insertions, 0 deletions
diff --git a/delta_E.py b/delta_E.py
new file mode 100644
index 0000000..cb47fd5
--- /dev/null
+++ b/delta_E.py
@@ -0,0 +1,103 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import tqdm
+
+from colour import (
+ STANDARD_OBSERVER_CMFS, SpectralShape, SpectralDistribution,
+ COLOURCHECKER_SDS, ILLUMINANT_SDS, ILLUMINANTS, sd_to_XYZ, XYZ_to_xy,
+ XYZ_to_Lab)
+from colour.difference import delta_E_CIE1976
+from colour.utilities import as_float_array
+from colour.plotting import plot_chromaticity_diagram_CIE1931
+
+from otsu2018 import load_Otsu2018_spectra
+from datasets.otsu2018 import *
+
+
+# Copied from https://github.com/enneract/colour/tree/feature/otsu2018
+def XYZ_to_sd_Otsu2018(
+ XYZ,
+ cmfs=STANDARD_OBSERVER_CMFS['CIE 1931 2 Degree Standard Observer']
+ .copy().align(OTSU_2018_SPECTRAL_SHAPE),
+ illuminant=ILLUMINANT_SDS['D65'].copy().align(
+ OTSU_2018_SPECTRAL_SHAPE),
+ clip=True):
+ XYZ = as_float_array(XYZ)
+ xy = XYZ_to_xy(XYZ)
+ cluster = select_cluster_Otsu2018(xy)
+
+ basis_functions = OTSU_2018_BASIS_FUNCTIONS[cluster]
+ mean = OTSU_2018_MEANS[cluster]
+
+ M = np.empty((3, 3))
+ for i in range(3):
+ sd = SpectralDistribution(
+ basis_functions[i, :],
+ OTSU_2018_SPECTRAL_SHAPE.range(),
+ )
+ M[:, i] = sd_to_XYZ(sd, illuminant=illuminant) / 100
+ M_inverse = np.linalg.inv(M)
+
+ sd = SpectralDistribution(mean, OTSU_2018_SPECTRAL_SHAPE.range())
+ XYZ_mu = sd_to_XYZ(sd, illuminant=illuminant) / 100
+
+ weights = np.dot(M_inverse, XYZ - XYZ_mu)
+ recovered_sd = np.dot(weights, basis_functions) + mean
+
+ if clip:
+ recovered_sd = np.clip(recovered_sd, 0, 1)
+
+ return SpectralDistribution(recovered_sd, OTSU_2018_SPECTRAL_SHAPE.range())
+
+
+
+if __name__ == '__main__':
+ print('Loading spectral data...')
+ data = load_Otsu2018_spectra('CommonData/spectrum_m.csv')
+ shape = SpectralShape(380, 730, 10)
+ sds = [SpectralDistribution(data[i, :], shape.range())
+ for i in range(data.shape[0])]
+
+
+ for name, colourchecker in COLOURCHECKER_SDS.items():
+ print('Adding %s...' % name)
+ sds += colourchecker.values()
+
+ D65 = ILLUMINANT_SDS['D65']
+ xy_w = ILLUMINANTS['CIE 1931 2 Degree Standard Observer']['D65']
+
+ x = []
+ y = []
+ errors = []
+ above_JND = 0
+ for i, sd in tqdm.tqdm(enumerate(sds), total=len(sds)):
+ XYZ = sd_to_XYZ(sd, illuminant=D65) / 100
+ xy = XYZ_to_xy(XYZ)
+ x.append(xy[0])
+ y.append(xy[1])
+ Lab = XYZ_to_Lab(XYZ, xy_w)
+
+ recovered_sd = XYZ_to_sd_Otsu2018(XYZ)
+ recovered_XYZ = sd_to_XYZ(recovered_sd, illuminant=D65) / 100
+ recovered_Lab = XYZ_to_Lab(recovered_XYZ, xy_w)
+
+ error = delta_E_CIE1976(Lab, recovered_Lab)
+ errors.append(error)
+ if error > 2.4:
+ above_JND += 1
+
+ print('Min. error: %g' % min(errors))
+ print('Max. error: %g' % max(errors))
+ print('Avg. error: %g' % np.mean(errors))
+ print('Errors above JND: %d (%.1f%%)'
+ % (above_JND, 100 * above_JND / len(sds)))
+
+ bins = [int((max(y) - min(y)) / 0.01), int((max(x) - min(x)) / 0.01)]
+ histogram, _, _ = np.histogram2d(x, y, bins=bins, weights=errors)
+
+ plot_chromaticity_diagram_CIE1931(standalone=False)
+ plt.imshow(histogram, extent=(min(x), max(x), min(y), max(y)),
+ interpolation='bicubic')
+ plt.colorbar()
+ plt.show()