1 files changed, 117 insertions, 0 deletions
diff --git a/crl/cri.py b/crl/cri.py
new file mode 100644
index 0000000..af62b1f
--- /dev/null
+++ b/crl/cri.py
@@ -0,0 +1,117 @@
+import numpy as np
+import scipy.spatial
+import crl.color as color
+import crl.tables as tables
+
+class CRIError(Exception):
+	pass
+
+class CRIResult:
+	def __init__(self, method_name):
+		self.method_name = method_name
+		self.SCRIs = []
+		self.tcs = []
+		self.tcs_ref = []
+
+	def __str__(self):
+		ret = f"CRI result for {self.method_name}\n"
+		ret += f"\tXYZ = {self.XYZ}\n"
+		ret += f"\tCCT = {self.CCT} K (delta uv = {self.duv})\n"
+		SCRIs = ", ".join("%.1f" % SCRI for SCRI in self.SCRIs)
+		ret += f"\tSCRIs = [{SCRIs}]\n"
+		ret += f"\tCRI = {self.CRI}"
+		return ret
+
+def CAT(sample, ref, test):
+	def c(u, v, _):
+		return (4 - u - 10 * v) / v
+	def d(u, v, _):
+		return (1.708 * v - 1.481 * u + 0.404) / v
+
+	cs, ds = c(*sample), d(*sample)
+	cr, dr = c(*ref), d(*ref)
+	ct, dt = c(*test), d(*test)
+
+	u = (10.872 + 0.404 * cr / ct * cs - 4 * dr / dt * ds) \
+	    / (16.518 + 1.481 * cr / ct * cs - dr / dt * ds)
+	v = 5.52 / (16.518 + 1.481 * cr / ct * cs - dr / dt * ds)
+
+	return color.uvY(u, v, sample[2])
+
+def cri_Ra(S, observer=color.Observer_2deg, ignore_errors=False):
+	res = CRIResult("CRI Ra")
+	res.S = S.normed()
+	res.XYZ = color.XYZ(res.S, observer=observer)
+	res.CCT, res.duv = color.cct(res.XYZ)
+
+	if res.duv > 5e-2 and not ignore_errors:
+		raise CRIError(f"spectrum not white enough with delta uv = {res.duv}")
+
+	if res.CCT < 5000:
+		white = color.BlackBodySpectrum(res.CCT).normed()
+	else:
+		if res.CCT > 25000 and ignore_errors:
+			res.CCT = 25000
+		white = color.CIEDaylightSpectrum(res.CCT).normed()
+
+	res.XYZ_ref = color.XYZ(white, observer=observer)
+
+	class UVWstarRef(color.UVWstar):
+		white_uvY = color.uvY(white, observer=observer)
+
+	class UVWstar(color.UVWstar):
+		white_uvY = color.uvY(S, observer=observer)
+
+	for i in range(14):
+		tcs = color.Spectrum(tables.tcs_ra[:, 0], tables.tcs_ra[:, 1 + i])
+
+		lit = color.uvY(color.CombinedSpectrum(tcs, res.S))
+		lit = CAT(lit, UVWstarRef.white_uvY, UVWstar.white_uvY)
+		lit = UVWstarRef(lit)
+		lit_ref = UVWstarRef(color.CombinedSpectrum(tcs, white), observer=observer)
+
+		res.tcs.append(lit)
+		res.tcs_ref.append(lit_ref)
+
+		dE = np.linalg.norm(lit.array - lit_ref.array)
+		res.SCRIs.append(100 - 4.6 * dE)
+
+	res.CRI = np.mean(res.SCRIs[:8])
+	res.ECRI = np.mean(res.SCRIs)
+	return res
+
+def cri_GAI(spec, observer=color.Observer_2deg):
+	points = []
+	for i in range(8):
+		tcs = color.Spectrum(color.tables.tcs_ra[:, 0], color.tables.tcs_ra[:, 1 + i])
+		u, v, _ = color.uvY76(color.CombinedSpectrum(tcs, spec), observer=observer)
+		points.append([u, v])
+	points = np.array(points)
+	hull = scipy.spatial.ConvexHull(points)
+
+	# The magic number is the hull.volume of the E illuminant
+	GAI = hull.volume * 100 / 0.007351716795637625
+
+	return GAI, hull.points[hull.vertices]
+
+def cri_FSI(spec):
+	wvl = np.arange(380, 730, 1)
+	Y = spec.Yi(wvl)
+	dwvl = wvl[1] - wvl[0]
+	Y = spec.Yi(wvl)
+	Y /= np.sum(Y * dwvl)
+
+	CEE = (wvl - 380) / (730 - 380)
+
+	Is = []
+	for i, _ in enumerate(wvl):
+		C = np.cumsum(Y * dwvl)
+		D = (C - CEE) ** 2
+		Is.append(sum(D * dwvl))
+		Y = np.roll(Y, 1)
+
+	FSI = np.mean(Is)
+	return FSI
+
+def cri_FSCI(spec):
+	return 100 - 5.1 * cri_FSI(spec)