summaryrefslogtreecommitdiff
path: root/crl/cri.py
blob: af62b1f4ec3bf8c8fd6e8ebe2c7ef0dd82d71907 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
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)