summaryrefslogtreecommitdiff
path: root/crl/color.py
diff options
context:
space:
mode:
Diffstat (limited to 'crl/color.py')
-rw-r--r--crl/color.py566
1 files changed, 566 insertions, 0 deletions
diff --git a/crl/color.py b/crl/color.py
new file mode 100644
index 0000000..d6240ff
--- /dev/null
+++ b/crl/color.py
@@ -0,0 +1,566 @@
+import numpy as np, scipy.interpolate, scipy.optimize
+import crl.tables as tables
+import crl.tables2 as tables2
+
+class Spectrum:
+ def __init__(self, X, Y):
+ self.X = X
+ self.Y = Y
+ self.Yi = scipy.interpolate.interp1d(self.X, self.Y, \
+ "linear", fill_value=0, bounds_error=False)
+
+ def dot(A, B, Xmin=380, Xmax=780, Xstep=0.01):
+ if not isinstance(B, Spectrum):
+ raise TypeError("B is an instance of %s, not Spectrum" \
+ % type(B))
+
+ if isinstance(A, PointSpectrum):
+ if isinstance(B, PointSpectrum):
+ raise TypeError("dotting two instances of PointSpectum isn't supported")
+ # ^ could be, but who cares?
+ return B.dot(A)
+
+ if isinstance(B, PointSpectrum):
+ return A.Yi(B.x_0) * B.S
+
+ X = np.arange(Xmin, Xmax, Xstep)
+ dX = X[1] - X[0]
+ return dX * np.sum(A.Yi(X) * B.Yi(X))
+
+ def normed(self, norm=100):
+ return NormedSpectrum(self, norm)
+
+ def __str__(self):
+ return "a tabulated Spectrum"
+
+class NormedSpectrum(Spectrum):
+ def __init__(self, spec, norm=100):
+ self.spec = spec
+
+ _, Y, _ = XYZ(self.spec)
+ self.factor = norm / Y
+
+ def Yi(self, wvl):
+ return self.factor * self.spec.Yi(wvl)
+
+ def __repr__(self):
+ return "NormedSpectrum(" + repr(self.spec) + ")"
+
+ def __str__(self):
+ return "normed " + str(self.spec)
+
+class PointSpectrum(Spectrum):
+ # the "S" stands for "sum" (being the integral of the Dirac delta)
+ def __init__(self, x_0, S=1):
+ self.x_0 = x_0
+ self.S = S
+
+ # a crappy hack for easy plotting
+ self.X = np.array([380, x_0 - 0.001, x_0, x_0 + 0.001, 780])
+ self.Y = np.array([0, 0, S, 0, 0])
+
+ def __repr__(self):
+ return "PointSpectrum(x_0=%g, S=%g)" % (self.x_0, self.S)
+
+ def __str__(self):
+ return self.__repr__()
+
+class BlackBodySpectrum(Spectrum):
+ def __init__(self, T):
+ self.T = T
+
+ def Yi(self, wvl):
+ h = 6.62607015e-34 # kg m^2 / s
+ c = 299792458 # m / s
+ k_B = 1.3806485e-23 # kg m^2 / K s^2
+ return 2 * h * c ** 2 / ((wvl / 1e9) ** 5) \
+ / (np.exp(h * c / (wvl / 1e9 * k_B * self.T)) - 1) \
+ * 1e-9 # kg / (nm s^3)
+
+ def __repr__(self):
+ return "BlackBodySpectrum(T=%g)" % self.T
+
+ def __str__(self):
+ return repr(self)
+
+# Schanda, page 41
+class CIEDaylightSpectrum(Spectrum):
+ S0 = Spectrum(tables.D_eigenvectors[:, 0], tables.D_eigenvectors[:, 1])
+ S1 = Spectrum(tables.D_eigenvectors[:, 0], tables.D_eigenvectors[:, 2])
+ S2 = Spectrum(tables.D_eigenvectors[:, 0], tables.D_eigenvectors[:, 3])
+
+ def __init__(self, T):
+ if T < 4000 or T > 25000:
+ raise ValueError("D illuminants are not defined for CCTs outside the 4000 K-25000 K range")
+ elif T <= 7000:
+ x = -4.6070e9 / T**3 \
+ + 2.9678e6 / T**2 \
+ + 0.09911e3 / T \
+ + 0.244063
+ else:
+ x = -2.0064e9 / T**3 \
+ + 1.9018e6 / T**2 \
+ + 0.0247483 / T \
+ + 0.237040
+
+ y = -3 * x**2 + 2.870 * x - 0.275
+
+ self.M1 = (-1.3515 - 1.7703 * x + 5.9114 * y) \
+ / (0.0241 + 0.2562 * x - 0.7341 * y)
+ self.M2 = (0.0300 - 31.4424 * x + 30.0717 * y) \
+ / (0.0241 + 0.2562 * x - 0.7341 * y)
+
+ self.M1 = round(self.M1, 3)
+ self.M2 = round(self.M2, 3)
+
+ def Yi(self, wvl):
+ return self.S0.Yi(wvl) \
+ + self.M1 * self.S1.Yi(wvl) \
+ + self.M2 * self.S2.Yi(wvl)
+
+class CombinedSpectrum(Spectrum):
+ def __init__(self, A, B):
+ self.A = A
+ self.B = B
+
+ def Yi(self, X):
+ return self.A.Yi(X) * self.B.Yi(X)
+
+ def __repr__(self):
+ return "CombinedSpectrum(%s, %s)" % (repr(self.A), repr(self.B))
+
+def combine_spectra(A, B):
+ if isinstance(B, PointSpectrum):
+ A, B = B, A
+
+ if isinstance(A, PointSpectrum):
+ if isinstance(B, PointSpectrum):
+ # technically this would make sense but would be pretty useless
+ raise TypeError("can't combine two instances of PointSpectrum")
+
+ return PointSpectrum(A.x_0, A.S * B.Ef(A.x_0))
+
+ return CombinedSpectrum(A, B)
+
+#
+# Tables, standard observers, basic colors
+#
+
+class Illuminants:
+ # 0 l,nm
+ # 1 Standard Illuminant A
+ # 2 Standard Illuminant D65
+ # 3 Illuminant C
+ # 4 IlluminantD50
+ # 5 IlluminantD55
+ # 6 IlluminantD75
+ A = Spectrum(tables.illuminants[:, 0], tables.illuminants[:, 1])
+ D50 = Spectrum(tables.illuminants[:, 0], tables.illuminants[:, 4])
+ D55 = Spectrum(tables.illuminants[:, 0], tables.illuminants[:, 5])
+ D65 = Spectrum(tables.illuminants[:, 0], tables.illuminants[:, 2])
+ D75 = Spectrum(tables.illuminants[:, 0], tables.illuminants[:, 6])
+
+ _wvl = np.linspace(380, 780, 801)
+ _E = np.ones(_wvl.size)
+ E = Spectrum(_wvl, _E)
+
+ _E = np.zeros(_wvl.size)
+ _E[120] = 3 # 440 nm
+ _E[376] = 5 # 568 nm
+ RedmanTest1 = Spectrum(_wvl, _E) # 6500K, CRI -12
+
+class Observer:
+ def __init__(self, name, wvl, cmf_x, cmf_y, cmf_z):
+ self.name = name
+ self.cmf_x = Spectrum(wvl, cmf_x)
+ self.cmf_y = Spectrum(wvl, cmf_y)
+ self.cmf_z = Spectrum(wvl, cmf_z)
+
+ def __str__(self):
+ return self.name
+
+Observer_2deg = Observer("2-deg Standard Observer",
+ *(tables.cmf[:, i] for i in range(4)))
+Observer_10deg = Observer("10-deg Standard Observer",
+ *(tables.cmf_1964[:, i] for i in range(4)))
+
+class Color:
+ full_name = "Color"
+ observer = Observer_2deg
+
+ # Call either:
+ # __init__(self, x, y, z)
+ # __init__(self, color)
+ # __init__(self, spectrum)
+ def __init__(self, *args, **kwargs):
+ observer = kwargs.get("observer", Observer_2deg)
+
+ if len(args) == 3:
+ self.array = np.array(args)
+ elif len(args) == 1:
+ if isinstance(args[0], Color):
+ self.array = args[0].to(self).array
+ elif isinstance(args[0], Spectrum):
+ self.array = self.from_spectrum(args[0], observer).array
+ else:
+ raise TypeError("argument is not a Color or a Spectrum")
+ else:
+ raise TypeError("bad number of arguments")
+
+ def __repr__(self):
+ return "%s(%g, %g, %g)" % (type(self).__name__, *self)
+
+ def __str__(self):
+ return "%s (%g, %g, %g), %s" % (self.full_name, *self,
+ self.observer)
+
+ @classmethod
+ def from_spectrum(cls, S, observer):
+ return cls.from_XYZ(XYZ.from_spectrum(S, observer))
+
+ def to(self, CS):
+ return CS.from_XYZ(self.to_XYZ())
+
+ @staticmethod
+ def from_XYZ(C):
+ raise NotImplementedError
+
+ def to_XYZ(C):
+ raise NotImplementedError
+
+ def __getitem__(self, index):
+ return self.array[index]
+
+ def __setitem__(self, index, value):
+ self.array[index] = value
+
+ def __iter__(self):
+ return self.array.__iter__()
+
+ def __add__(A, B):
+ return type(A)(*(A.array + B.array))
+
+class sRGB(Color):
+ full_name = "sRGB"
+ matrix = np.array([
+ [ 3.240479, -1.537150, -0.498353],
+ [-0.969256, 1.875992, 0.041556],
+ [ 0.055648, -0.204043, 1.057331]])
+ matrix_inverse = np.linalg.inv(matrix)
+
+ @classmethod
+ def from_XYZ(cls, C):
+ def gamma(x):
+ if x <= 0.0031308:
+ return 12.92 * x
+ return 1.055 * x ** (1 / 2.4) - 0.055
+
+ RGB = np.dot(cls.matrix, C.array)
+ return sRGB(*[gamma(x) for x in RGB])
+
+ def to_XYZ(self):
+ def gamma_inverse(x):
+ if x <= 0.040449936:
+ return x / 12.92;
+ return ((x + 0.055) / 1.055) ** 2.4
+
+ RGB = [gamma_inverse(x) for x in self]
+ return XYZ(*np.dot(self.matrix_inverse, RGB))
+
+class sRGB100(sRGB):
+ full_name = "sRGB (with Y=100 as white)"
+ matrix = 0.01 * sRGB.matrix
+ matrix_inverse = np.linalg.inv(matrix)
+
+#
+# CIE 1931
+#
+
+class XYZ(Color):
+ full_name = "CIE 1931 XYZ"
+
+ @staticmethod
+ def from_XYZ(C):
+ return C
+
+ def to_XYZ(self):
+ return self
+
+ @staticmethod
+ def from_spectrum(S, observer):
+ X = S.dot(observer.cmf_x)
+ Y = S.dot(observer.cmf_y)
+ Z = S.dot(observer.cmf_z)
+ return XYZ(X, Y, Z)
+
+ def normed(self, norm=100):
+ X, Y, Z = self
+ return XYZ(X * norm / Y, norm, Z * norm / Y)
+
+class xyY(Color):
+ full_name = "CIE 1931 xyY"
+
+ @staticmethod
+ def from_XYZ(C):
+ x = C[0] / sum(C)
+ y = C[1] / sum(C)
+ return xyY(x, y, C[1])
+
+ def to_XYZ(self):
+ x, y, Y = self
+ try:
+ X = Y / y * x
+ Z = Y / y * (1 - x - y)
+ except ZeroDivisionError:
+ X = 0
+ Z = 0
+ return XYZ(X, Y, Z)
+
+#
+# CIE 1960
+#
+
+class UVW(Color):
+ full_name = "CIE 1960 UCS (McAdam) UVW"
+
+ @staticmethod
+ def from_XYZ(C):
+ X, Y, Z = C
+ U = 2 * X / 3
+ V = Y
+ W = 1/2 * (-X + 3 * Y + Z)
+ return UVW(U, V, W)
+
+ def to_XYZ(self):
+ U, V, W = self
+ X = 3/2 * U
+ Y = V
+ Z = 3/2 * U - 3 * V + 2 * W
+ return XYZ(X, Y, Z)
+
+ def to_uvY(self):
+ U, V, W = self
+ u = U / sum(UVW)
+ v = V / sum(UVW)
+ return uvY(u, v, V)
+
+class uvY(Color):
+ full_name = "CIE 1960 uvY"
+
+ @staticmethod
+ def from_XYZ(C):
+ X, Y, Z = C
+ u = 4 * X / (X + 15 * Y + 3 * Z)
+ v = 6 * Y / (X + 15 * Y + 3 * Z)
+ return uvY(u, v, Y)
+
+ def to_UVW(self):
+ u, v, V = self
+ try:
+ U = V * u / v
+ W = -(V * v + V * u - V) / v
+ except ZeroDivisionError:
+ U = 0
+ W = 0
+ return UVW(U, V, W)
+
+ def to_XYZ(self):
+ return self.to_UVW().to_XYZ()
+
+#
+# CIE 1964
+#
+
+class UVWstar(Color):
+ full_name = "CIE 1964 U*V*W*"
+
+ @classmethod
+ def from_uvY(cls, C, CAT=True):
+ u, v, Y = C
+ un, vn, _ = cls.white_uvY if CAT else (0, 0, 0)
+
+ Wstar = 25 * Y ** (1/3) - 17
+ Ustar = 13 * Wstar * (u - un)
+ Vstar = 13 * Wstar * (v - vn)
+
+ return cls(Ustar, Vstar, Wstar)
+
+ @classmethod
+ def from_XYZ(cls, C):
+ return cls.from_uvY(C.to(uvY))
+
+ def to_uvY(self, CAT=True):
+ Ustar, Vstar, Wstar = self
+ un, vn, _ = self.white_uvY if CAT else (0, 0, 0)
+
+ u = Ustar / (13 * Wstar) + un
+ v = Vstar / (13 * Wstar) + vn
+ Y = (Wstar**3 + 51 * Wstar**2 + 867 * Wstar + 4913) / 15625
+
+ return uvY(u, v, Y)
+
+ def to_XYZ(self):
+ return self.to_uvY().to(XYZ)
+
+
+
+#
+# CIE 1976
+#
+
+class uvY76(Color):
+ full_name = "CIE 1976 u'v'Y"
+
+ @staticmethod
+ def from_XYZ(C):
+ X, Y, Z = C
+ u = 4 * X / (X + 15 * Y + 3 * Z)
+ v = 9 * Y / (X + 15 * Y + 3 * Z)
+ return uvY76(u, v, Y)
+
+ def to_XYZ(self):
+ u, v, Y = self
+ X = 9 * u / (4 * v) * Y
+ Z = (12 - 3 * u - 20 * v) / (4 * v) * Y
+ return XYZ(X, Y, Z)
+
+# The magic CIE 1976 function for mapping Y to L*
+def CIE1976_f(t):
+ if t > (6/29) ** 3:
+ return t ** (1/3)
+
+ return 841 * t / 108 + 4/29
+
+# The inverse of the above
+def CIE1976_f_inverse(t):
+ if t > 6/29:
+ return t ** 3
+
+ return (3132 * t - 432) / 24389
+
+class Luvstar(Color):
+ full_name = "CIE 1976 L*u*v*"
+
+ @classmethod
+ def from_uvY76(cls, C):
+ u, v, Y = C
+ un, vn, Yn = cls.white_uvY76
+
+ Lstar = 116 * CIE1976_f(Y / Yn) - 16
+ ustar = 13 * Lstar * (u - un)
+ vstar = 13 * Lstar * (v - vn)
+
+ return cls(Lstar, ustar, vstar)
+
+ @classmethod
+ def from_XYZ(cls, C):
+ return cls.from_uvY76(C.to(uvY76))
+
+ def to_uvY76(self):
+ Lstar, ustar, vstar = self
+ un, vn, Yn = self.white_uvY76
+
+ u = ustar / (13 * Lstar) + un
+ v = vstar / (13 * Lstar) + vn
+ Y = Yn * CIE1976_f_inverse((Lstar + 16) / 116)
+
+ return uvY76(u, v, Y)
+
+ def to_XYZ(self):
+ return self.to_uvY76().to_XYZ()
+
+class Labstar(Color):
+ full_name = "CIE 1976 L*a*b*"
+
+ @classmethod
+ def from_XYZ(cls, C):
+ X, Y, Z = C
+ Xn, Yn, Zn = cls.white
+
+ Lstar = 116 * CIE1976_f(Y / Yn) - 16
+ astar = 500 * (CIE1976_f(X / Xn) - CIE1976_f(Y / Yn))
+ bstar = 200 * (CIE1976_f(Y / Yn) - CIE1976_f(Z / Zn))
+
+ return cls(Lstar, astar, bstar)
+
+ def to_XYZ(self):
+ Lstar, astar, bstar = self
+ Xn, Yn, Zn = self.white
+
+ X = Xn * CIE1976_f_inverse(astar / 500 + (Lstar + 16) / 116)
+ Y = Yn * CIE1976_f_inverse((Lstar + 16) / 116)
+ Z = Zn * CIE1976_f_inverse(-bstar / 200 + (Lstar + 16) / 116)
+
+ return XYZ(X, Y, Z)
+
+#
+# Some useful predefined color spaces
+#
+
+class UVWstarD65(UVWstar):
+ white_uvY = uvY(Illuminants.D65)
+
+class LuvstarD65(Luvstar):
+ white_uvY76 = uvY76(Illuminants.D65)
+
+class LabstarD65(Labstar):
+ white = XYZ(Illuminants.D65).normed()
+
+class UVstarNormed(Color):
+ full_name = "CIE 1964 U*V*W*, normed (plots only!)"
+
+ @staticmethod
+ def from_UVWstar(C):
+ if not issubclass(type(C), UVWstar):
+ raise TypeError(f"{type(C)} is not a subclass of UVWstar")
+
+ Ustar, Vstar, Wstar = C
+
+ Unorm = Ustar / Wstar
+ Vnorm = Vstar / Wstar
+
+ return UVstarNormed(Unorm, Vnorm, Wstar)
+
+ @staticmethod
+ def from_XYZ(C):
+ raise NotImplementedError("use the from_UVWstar method directly")
+
+ def to_UVWstar(self):
+ Unorm, Vnorm, Wstar = self
+
+ Ustar = Unorm * Wstar
+ Vstar = Vnorm * Wstar
+
+ return UVWstarD65(Ustar, Vstar, Wstar)
+
+ def to_XYZ(self):
+ return self.to_UVWstar().to_XYZ()
+
+class UVstarNormedD65(UVstarNormed):
+ @staticmethod
+ def from_XYZ(C):
+ return UVstarNormedD65.from_UVWstar(UVWstarD65(C))
+
+#
+# CCT
+#
+
+def cct(C):
+ u, v, _ = uvY(C)
+ uv = np.array([u, v])
+
+ uv_planck = scipy.interpolate.interp1d(
+ tables2.locus_uv[:, 0], tables2.locus_uv[:, 1:3], axis=0,
+ bounds_error=False, fill_value="extrapolate")
+
+ def delta(T):
+ return np.linalg.norm(uv_planck(T) - uv)
+
+ res = scipy.optimize.minimize_scalar(delta, bounds=[1000, 20000])
+ return res.x, res.fun
+
+def cct_McCamy(C):
+ x, y, _ = xyY(C)
+
+ n = (x - 0.3320) / (y - 0.1858)
+ return -449 * n**3 + 3525 * n**2 - 6823.3 * n + 5520.33
generated by cgit (git 2.34.1) at 2024-05-02 10:51:44 +0000