MatrixM-for-standard-white-A.py

import numpy as np

# Source: http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html

# chromaticity coordinates
# values are for sRGB (http://www.brucelindbloom.com/index.html?WorkingSpaceInfo.html)
xr, xg, xb = 0.6400, 0.3000, 0.1500
yr, yg, yb = 0.3300, 0.6000, 0.0600

# Reference Whites in XYZ, source: https://de.mathworks.com/help/images/ref/whitepoint.html
# https://en.wikipedia.org/wiki/Standard_illuminant#White_points_of_standard_illuminants
# A: Simulates typical, domestic, tungsten-filament lighting with correlated color temperature of 2856 K.
A = np.array([1.0985, 1.0000, 0.3558])
# C: Simulates average or north sky daylight with correlated color temperature of 6774 K.
C = np.array([0.9807, 1.0000, 1.1822])
# D50: Simulates warm daylight at sunrise or sunset with correlated color temperature of 5003 K. Also known as horizon light.
D50 = np.array([0.9642, 1.0, 0.8251])
# D55: Simulates mid-morning or mid-afternoon daylight with correlated color temperature of 5500 K.
D55 = np.array([0.9568, 1.0, 0.9214])
# D65: Simulates noon daylight with correlated color temperature of 6504 K.
D65 = np.array([0.95047, 1.0, 1.08883])

# Ideal would be 2700 Kelvin as reference white, but I didn't find a value for this online
# So instead we take standard illuminant A, which has 2856K and is closest to 2700K
W = A

X = np.array([xr/yr, xg/yg, xb/yb])
Y = np.array([1, 1, 1])
Z = np.array([(1-xr-yr)/yr, (1-xg-yg)/yg, (1-xb-yb)/yb])

print("X", X)
print("Y", Y)
print("Z", Z)
print("W", W)

XYZ = np.vstack([X, Y, Z])
print("XZY", XYZ)

XYZinv = np.linalg.inv(XYZ)
print("XYZ-1", XYZinv)

S = np.dot(XYZinv, W)
print("S", S)

S3 = np.transpose(np.repeat(S, repeats=3).reshape(3,3))
print("S3", S3)

M = np.multiply(S3, XYZ)
print("M", M)