calculate temperature with unit Kelvin

.gitignore

@@ -6,9 +6,6 @@
 *.ipynb
 *.o
 *.so
-*.py
-*.sh
-*.txt
 .vscode/
 build/
 dist/

include/color.h

@@ -73,6 +73,7 @@ constexpr double rad2Deg(const double rad) { return ((180.0 / M_PI) * rad); }
 
 RGB XYZtoRGB(const XYZ &xyz);
 XYZ LABtoXYZ(const LAB &lab);
+LAB XYZtoLAB(const XYZ &xyz);
 CMY RGBtoCMY(const RGB &rgb);
 
 inline RGB LABtoRGB(const LAB &lab) { return XYZtoRGB(LABtoXYZ(lab)); }
@@ -89,6 +90,14 @@ inline LAB LCHtoLAB(const LCH &lch) {
   return LAB{lch.l, lch.c * cos(theta), lch.c * sin(theta)};
 }
 
+// @param T temperature in Kelvin 4000 < T < 25000
+// @param y luminocity in XYZ
+XYZ IlluminantDKelvin(double T, double y);
+
+// @param cx chromaticity x in xyY
+// @param y luminocity in XYZ
+XYZ IlluminantDChromaticity(double cx, double y);
+
 } // namespace color
 
 std::ostream &operator<<(std::ostream &s, const color::LAB &lab);

python/__init__.py

@@ -1,2 +1,2 @@
-from .native import LABtoRGB, LABtoLCH, LCHtoLAB, maxChroma, perception
+from .native import *
 from .version import __version__

python/cli.py

@@ -60,11 +60,11 @@ def argparser():
     ps.add_argument(
         '-T',
         '--temperature',
-        default=0,
-        dest='dK',
-        metavar='dK',
+        default=5000,
+        dest='T',
+        metavar='T',
         type=float,
-        help='temperature hint, -1(warm) < dK < 1(cold)')
+        help='temperature 4000 < T < 25000')
     ps.add_argument(
         '-p',
         '--profile',

python/theme.py

@@ -1,8 +1,7 @@
 from itertools import product
-from .native import LABtoRGB, LCHtoLAB, maxChroma, perception
+from .native import LABtoRGB, LCHtoLAB, maxChroma, perception, IlluminantD
 
-D65_A = -1.65 / 100  # a component of D65 in LAB
-D65_B = -19.33 / 100  # b component of D65 in LAB
+# D50_x = 0.3457048
 SEMANTIC_HUE = {'red': 30, 'yellow': 80, 'green': 120, 'blue': 260}
 
 
@@ -22,15 +21,11 @@ def interpolate2D(x1, y1, z1, x2, y2, z2, x3, y3, z3):
     return lambda x, y: (c - y * b - x * a) / d
 
 
-def tempered_gray(L, dK):
-    return (L, L * D65_A * dK, L * D65_B * dK)
-
-
 # @param args.Lf
 # @param args.Lb
 # @param args.L
 # @param args.maxC
-# @param args.dK
+# @param args.T
 # @param args.verbose
 def update_context(args, ctx):
     if args.Lb < 0:
@@ -40,7 +35,7 @@ def update_context(args, ctx):
     elif args.Lf < 0:
         args.Lf = 100 - args.Lb
 
-    ctx['fg-hex'] = rgb_hex(LABtoRGB((tempered_gray(args.Lf, args.dK))))
+    ctx['fg-hex'] = rgb_hex(LABtoRGB((IlluminantD(args.T, args.Lf))))
     ctx['ui-theme'] = 'vs-dark' if args.Lb < 50 else 'vs'
 
     # L3: text color, very close to Lf
@@ -63,14 +58,14 @@ def update_context(args, ctx):
     if (args.verbose):
         print(
             f'Init parameters: Lf={args.Lf} L3={args.L3} args.L2={args.L2} args.L1={args.L1} '
-            + f'args.L0={args.L0} Lb={args.Lb} dK={args.dK}')
+            + f'args.L0={args.L0} Lb={args.Lb} T={args.T}')
 
     sgn = 1 if args.Lf > args.Lb else -1
     palette3 = perception(
         7,
         L=args.L3,
         maxC=args.maxC,
-        fixed=[tempered_gray(args.L3, args.dK)],
+        fixed=[IlluminantD(args.T, args.L3)],
         quiet=not args.verbose)
 
     for i, rgb in enumerate(palette3['rgb']):
@@ -84,8 +79,8 @@ def update_context(args, ctx):
 
     for i, delta in enumerate([15, 25, 60]):
         ctx[f'line-{i}-hex'] = rgb_hex(
-            LABtoRGB(tempered_gray(args.Lb + sgn * delta, args.dK)))
+            LABtoRGB(IlluminantD(args.T, args.Lb + sgn * delta)))
 
     for i, delta in enumerate([0, 5, 12]):
         ctx[f'bg-{i}-hex'] = rgb_hex(
-            LABtoRGB(tempered_gray(args.Lb + sgn * delta, args.dK)))
+            LABtoRGB(IlluminantD(args.T, args.Lb + sgn * delta)))

src/colorspace.cpp

@@ -3,10 +3,9 @@
 #include <constant.h>
 
 namespace color {
-double sgn(double val);
-constexpr inline double pow3(double x);
 
 double sgn(double val) { return (0. < val) - (val < 0.); }
+constexpr inline double pow2(double x) { return x * x; }
 constexpr inline double pow3(double x) { return x * x * x; }
 
 RGB XYZtoRGB(const XYZ &xyz) {
@@ -40,6 +39,22 @@ XYZ LABtoXYZ(const LAB &lab) {
   return XYZ{x, y, z};
 }
 
+LAB XYZtoLAB(const XYZ &xyz) {
+  double x = xyz.x / PCS_X;
+  double y = xyz.y / PCS_Y;
+  double z = xyz.z / PCS_Z;
+
+  double fx = x > 0.008856 ? pow(x, 1. / 3.) : (7.787 * x + 16. / 116.);
+  double fy = y > 0.008856 ? pow(y, 1. / 3.) : (7.787 * y + 16. / 116.);
+  double fz = z > 0.008856 ? pow(z, 1. / 3.) : (7.787 * z + 16. / 116.);
+
+  double l = 116 * fy - 16;
+  double a = 500 * (fx - fy);
+  double b = 200 * (fy - fz);
+
+  return LAB{l, a, b};
+}
+
 CMY RGBtoCMY(const RGB &rgb) {
   const double cr = 1 - rgb.r;
   const double cg = 1 - rgb.g;
@@ -51,6 +66,24 @@ CMY RGBtoCMY(const RGB &rgb) {
   return CMY{c, m, y};
 }
 
+// http://www.brucelindbloom.com/index.html?Eqn_DIlluminant.html
+XYZ IlluminantDKelvin(double T, double y) {
+  double cx = T < 7000 ? (-4.6070e9 / pow3(T) + 2.9678e6 / pow2(T) +
+                          0.09911e3 / T + 0.244063)
+                       : (-2.0064e9 / pow3(T) + 1.9018e6 / pow2(T) +
+                          0.24748e3 / T + 0.237040);
+  return IlluminantDChromaticity(cx, y);
+}
+
+XYZ IlluminantDChromaticity(double cx, double y) {
+  double cy = -3 * pow2(cx) + 2.87 * cx - 0.275;
+
+  double x = cx * y / cy;
+  double z = (1 - cx - cy) * y / cy;
+
+  return XYZ{x, y, z};
+}
+
 } // namespace color
 
 std::ostream &operator<<(std::ostream &s, const color::LAB &lab) {

src/python.cpp

@@ -35,6 +35,16 @@ static PyObject *LABtoRGB(PyObject *self, PyObject *args) {
   return Py_BuildValue("(ddd)", rgb.r, rgb.g, rgb.b);
 }
 
+static PyObject *LABtoXYZ(PyObject *self, PyObject *args) {
+  color::LAB lab;
+
+  if (!PyArg_ParseTuple(args, "(ddd)", &lab.l, &lab.a, &lab.b))
+    return nullptr;
+
+  const auto xyz = color::LABtoXYZ(lab);
+  return Py_BuildValue("(ddd)", xyz.x, xyz.y, xyz.z);
+}
+
 static PyObject *LABtoLCH(PyObject *self, PyObject *args) {
   color::LAB lab;
 
@@ -55,6 +65,30 @@ static PyObject *LCHtoLAB(PyObject *self, PyObject *args) {
   return Py_BuildValue("(ddd)", lab.l, lab.a, lab.b);
 }
 
+static PyObject *IlluminantDKelvin(PyObject *self, PyObject *args) {
+  double T;
+  double l;
+
+  if (!PyArg_ParseTuple(args, "dd", &T, &l))
+    return nullptr;
+
+  const auto lab = color::XYZtoLAB(
+      color::IlluminantDKelvin(T, color::LABtoXYZ({l, 0, 0}).y));
+  return Py_BuildValue("(ddd)", lab.l, lab.a, lab.b);
+}
+
+static PyObject *IlluminantDChromaticity(PyObject *self, PyObject *args) {
+  double cx;
+  double l;
+
+  if (!PyArg_ParseTuple(args, "dd", &cx, &l))
+    return nullptr;
+
+  const auto lab = color::XYZtoLAB(
+      color::IlluminantDChromaticity(cx, color::LABtoXYZ({l, 0, 0}).y));
+  return Py_BuildValue("(ddd)", lab.l, lab.a, lab.b);
+}
+
 static PyObject *pyperception(PyObject *self, PyObject *args, PyObject *kw) {
   // color::LAB fg, bg;
   unsigned long M;
@@ -160,6 +194,12 @@ static PyMethodDef methods[] = {
      "\n"
      "@param lab (l,a,b)\n"},
 
+    {"LABtoXYZ", LABtoXYZ, METH_VARARGS,
+     "convert LAB to XYZ\n"
+     "------------------\n"
+     "\n"
+     "@param lab (l,a,b)\n"},
+
     {"LABtoLCH", LABtoLCH, METH_VARARGS,
      "convert LAB to LCH\n"
      "------------------\n"
@@ -172,6 +212,13 @@ static PyMethodDef methods[] = {
      "\n"
      "@param lch (l,c,h)\n"},
 
+    {"IlluminantD", IlluminantDKelvin, METH_VARARGS,
+     "Calculate Illuminant D\n"
+     "----------------------\n"
+     "\n"
+     "@param T temperature in Kelvin, 4000 < T < 25000\n"
+     "@param L luminocity in LAB"},
+
     {nullptr, nullptr, 0, nullptr}};
 
 static struct PyModuleDef module = {

test/testcolorspace.cpp

@@ -66,6 +66,12 @@ int testcolor() {
   check(rgb5.g < 0, "RGB limit");
   check(rgb5.b > 1, "RGB limit");
 
+  color::LAB lab3{50, 10, 20};
+  auto lab3a = XYZtoLAB(LABtoXYZ(lab3));
+  check(lab3a.l, lab3.l, "LAB->XYZ->LAB");
+  check(lab3a.a, lab3.a, "LAB->XYZ->LAB");
+  check(lab3a.b, lab3.b, "LAB->XYZ->LAB");
+
   check(color::LABtoLCH(color::LAB{50, 20, 0}).h, 0, "LCH hue");
   check(color::LABtoLCH(color::LAB{50, 0, 20}).h, 90, "LCH hue");
   check(color::LABtoLCH(color::LAB{50, -20, 0}).h, 180, "LCH hue");
@@ -78,6 +84,11 @@ int testcolor() {
   check(color::LCHtoLAB(color::LCH{50, 20, 225}).b, -20. / sqrt(2),
         "LCH -> LAB");
 
+  auto d65 = color::IlluminantDKelvin(6500, 1);
+  check(d65.x, 0.95047, "D65");
+  check(d65.y, 1, "D65");
+  check(d65.z, 1.08833, "D65");
+
   std::cout << std::endl;
   int ret = EXIT_SUCCESS;
   for (unsigned int i = 0; i < passFail.size(); i++) {

test/testilluminantd.py

@@ -0,0 +1,6 @@
+#! /usr/bin/env python3
+
+from perception import IlluminantD
+
+print(IlluminantD(5000, 100))
+print(IlluminantD(6500, 100))