diff --git a/library/colorspaces.py b/library/colorspaces.py new file mode 100644 index 0000000..94a2168 --- /dev/null +++ b/library/colorspaces.py @@ -0,0 +1,250 @@ +import numpy as np + +# UTILITY {{{1 + +F = lambda a: np.array(a, np.float64) +A = lambda *args: F(args) +and1 = lambda a: A(*a.T, np.ones_like(a.T[0])).T # appends a column of 1s +# (works for both vectors and matrices thanks to some transposition magic) + + +def either1d2d(fun): + from functools import wraps + + @wraps(fun) + def wrap(a, *args, **kwargs): + a = np.asfarray(a) + assert a.ndim in (1, 2), f"expected a 1D or 2D array, got {a.ndim}" + return fun(a, *args, **kwargs) + + return wrap + + +# DATA {{{1 + + +d65 = A(95.0489, 100.0000, 108.8840) +d50 = A(96.4212, 100.0000, +82.5188) + +mat_alt_xyz = A( # for SRLAB2 + (0.320530, 0.636920, 0.042560), + (0.161987, 0.756636, 0.081376), + (0.017228, 0.108660, 0.874112), +) + +mat_alt_lab = A( # for SRLAB2 + (37.0950, +62.9054, -0.0008), + (663.4684, -750.5078, +87.0328), + (63.9569, +108.4576, -172.4152), +) + +lab_magic = A( + [0, 116, 0], + [500, -500, 0], + [0, 200, -200], +) + +lab_offset = A(-16, 0, 0) + +bradford = A( + [+0.8951, +0.2664, -0.1614], + [-0.7502, +1.7135, +0.0367], + [+0.0389, -0.0685, +1.0296], +) + +mat_BFD = bradford + +mat_97 = A( + [+0.8562, +0.3372, -0.1934], + [-0.8360, +1.8327, +0.0033], + [+0.0357, -0.0469, +1.0112], +) + +mat_02 = A( + [+0.7328, +0.4296, -0.1624], + [-0.7036, +1.6975, +0.0061], + [+0.0030, +0.0136, +0.9834], +) + +mat_16 = A( + [+0.401288, +0.650173, -0.051461], + [-0.250268, +1.204414, +0.045854], + [-0.002079, +0.048952, +0.953127], +) + +# REC709_PRI +sRGB_primaries = A( # strangely, these seem to be missing a digit of significance. + # actually, this might be ITU-R BT.709-5 which is naturally less accurate? + [+0.64000, +0.33000], # +0.03000 + [+0.30000, +0.60000], # +0.10000 + [+0.15000, +0.06000], # +0.79000 + [+0.31270, +0.32900], # +0.35830 +) + +# AP0 +AP0_primaries = A( + [+0.73470, +0.26530], # +0.00000 + [+0.00000, +1.00000], # +0.00000 + [+0.00010, -0.07700], # +1.07690 + [+0.32168, +0.33767], # +0.34065 +) + +# AP1 +AP1_primaries = A( + [+0.71300, +0.29300], # -0.00600 + [+0.16500, +0.83000], # +0.00500 + [+0.12800, +0.04400], # +0.82800 + [+0.32168, +0.33767], # +0.34065 +) + +ref_XYZ = A( # sRGB -> XYZ (still needs to be converted from D65) + [506752 / 1228815, 87881 / 245763, 12673 / 70218], # sum: 3127/3290 + [87098 / 409605, 175762 / 245763, 12673 / 175545], # sum: 1/1 + [7918 / 409605, 87881 / 737289, 1001167 / 1053270], # sum: 3583/3290 +) + +_Zify = A( # goes *after* input vector, not before. + [1, 0, -1], + [0, 1, -1], + [0, 0, +1], +) + +# DATA (INVERSES) {{{1 + +inv_alt_xyz = np.linalg.inv(mat_alt_xyz) +inv_alt_lab = np.linalg.inv(mat_alt_lab) +inv_lab_magic = np.linalg.inv(lab_magic) + +inv_BFD = np.linalg.inv(mat_BFD) # TODO: use this! +inv_97 = np.linalg.inv(mat_97) # TODO: rename to _CAT97? +inv_02 = np.linalg.inv(mat_02) # TODO: rename to _CAT02? +inv_16 = np.linalg.inv(mat_16) + + +# FUNCTIONS {{{1 + + +@either1d2d +def cartesian_to_polar(a): + return A(a.T[0], np.hypot(a.T[1], a.T[2]), np.arctan2(a.T[2], a.T[1])).T + + +@either1d2d +def polar_to_cartesian(a): + return A(a.T[0], a.T[1] * np.cos(a.T[2]), a.T[1] * np.sin(a.T[2])).T + + +def rgb2lin(a): + a = np.asfarray(a) + return np.where(a > 0.04045, ((a + 0.055) / 1.055) ** 2.4, a / 12.92) + + +def lin2rgb(a): + a = np.asfarray(a) + return np.where(a > 0.04045 / 12.92, (a ** (1 / 2.4) * 1.055) - 0.055, a * 12.92) + + +@either1d2d +def xyY_2_XYZ(xyY): + # transposition allows this function to be vectorized (over vectors). + x, y, Y = xyY.T + y = np.maximum(y, 1e-10) + return (Y / y * A(x, y, 1.0 - x - y)).T + + +@either1d2d +def lab_to_xyz(lab, *, ref): + ijk = (lab - lab_offset) @ inv_lab_magic.T + xyz = np.where(ijk > np.cbrt(0.008856), ijk**3, (ijk - 16 / 116) / 7.787) + xyz = xyz * ref + return xyz + + +@either1d2d +def xyz_to_lab(xyz, *, ref): + xyz = xyz / ref + # absolute value just to shut up numpy: + ijk = np.where(xyz > 0.008856, np.cbrt(np.abs(xyz)), 7.787 * xyz + 16 / 116) + lab = ijk @ lab_magic.T + lab_offset + return lab + + +# print(xyz_to_lab((30, 50, 20), ref=d50)) +# print("[ 76.06926101 -58.04284385 34.04319485]") +# print(lab_to_xyz(( 76.06926101, -58.04284385, 34.04319485), ref=d50)) +# stop + + +@either1d2d +def xyz_to_srlab2(xyz): + rgb = xyz @ inv_XYZ.T + xyz = rgb @ mat_alt_xyz.T + ijk = np.where( + xyz <= 216.0 / 24389.0, 24389.0 / 2700.0 * xyz, 1.16 * np.cbrt(xyz) - 0.16 + ) + lab = ijk @ mat_alt_lab.T + return lab + + +@either1d2d +def srlab2_to_xyz(lab): + ijk = lab @ inv_alt_lab.T + xyz = np.where(ijk <= 0.08, 2700.0 / 24389.0 * ijk, ((ijk + 0.16) / 1.16) ** 3) + rgb = xyz @ inv_alt_xyz.T + xyz = rgb @ mat_XYZ.T + return xyz + + +# print(srlab2_to_xyz(xyz_to_srlab2((0.3, 0.4, 0.5)))) +# print(srlab2_to_xyz(xyz_to_srlab2((0.3, 0.4, 0.5)))) +# stop + + +def ref_to_uv(ref): + # NOTE: this assumes that ref is valid i.e. won't cause a division by zero. + xt, yt, zt = ref # tristimulus values (CIE 1931 XYZ) + den = xt + yt + zt # denominator + xc, yc = xt / den, yt / den # chromaticity values (still CIE 1931) + den = xc + 15 * yc + 3 * (1 - xc - yc) + u, v = 4 * xc / den, 6 * yc / den # u, v coordinates (CIE 1960 UCS) + return u, v + + +@either1d2d +def xyz_to_uvw(xyz, *, ref): + u0, v0 = ref_to_uv(ref) + x, y, z = xyz.T + den = np.maximum(x + 15 * y + 3 * z, 1e-12) + u = 4 * x / den + v = 6 * y / den + W = 25 * np.cbrt(y) - 17 + U = 13 * W * (u - u0) + V = 13 * W * (v - v0) + return A(U, V, W).T + + +@either1d2d +def uvw_to_xyz(uvw, *, ref): + u0, v0 = ref_to_uv(ref) + U, V, W = uvw.T + Y = ((W + 17) / 25) ** 3 + u = U / (13 * W) + u0 + v = V / (13 * W) + v0 + den = Y * 6 / np.maximum(v, 1e-12) + X = u / 4 * den + Z = (den - X - 15 * Y) / 3 + return A(X, Y, Z).T + + +def makemat(prim, normalize=True): + primz, white = and1(prim[:3]) @ _Zify, xyY_2_XYZ(and1(prim[3])) + rescalant = np.cross(primz[((1, 2, 0),)], primz[((2, 0, 1),)]) @ white + mat = (primz * rescalant[:, None]).T + # normalize such that max inputs (1.0, 1.0, 1.0) yields (any, 1.0, any) + return mat / np.sum(mat[1, :]) if normalize else mat + + +# GENERATED DATA {{{1 + +mat_XYZ = makemat(sRGB_primaries) +inv_XYZ = np.linalg.inv(mat_XYZ)