Update jedypod

using oklab, a little better than cie-lab, but still not mathematically correct.

Author: natural-harmonia-gropius

gamut-mapping/jedypod.glsl

@@ -126,84 +126,63 @@ vec3 XYZ_to_RGB(vec3 XYZ) {
     return XYZ * M;
 }
 
-vec3 XYZD65_to_XYZD50(vec3 XYZ) {
+vec3 XYZ_to_LMS(vec3 XYZ) {
     mat3 M = mat3(
-         1.0479298208405488,   0.022946793341019088, -0.05019222954313557,
-         0.029627815688159344, 0.990434484573249,    -0.01707382502938514,
-        -0.009243058152591178, 0.015055144896577895,  0.7518742899580008);
+        0.8190224432164319,   0.3619062562801221,  -0.12887378261216414,
+        0.0329836671980271,   0.9292868468965546,   0.03614466816999844,
+        0.048177199566046255, 0.26423952494422764,  0.6335478258136937);
     return XYZ * M;
 }
 
-vec3 XYZD50_to_XYZD65(vec3 XYZ) {
+vec3 LMS_to_XYZ(vec3 LMS) {
     mat3 M = mat3(
-         0.9554734527042182,   -0.023098536874261423, 0.0632593086610217,
-        -0.028369706963208136,  1.0099954580058226,   0.021041398966943008,
-         0.012314001688319899, -0.020507696433477912, 1.3303659366080753);
-    return XYZ * M;
-}
-
-float delta = 6.0 / 29.0;
-float deltac = delta * 2.0 / 3.0;
-
-float f1(float x, float delta) {
-    return x > pow(delta, 3.0) ?
-        cbrt(x) :
-        deltac + x / (3.0 * pow(delta, 2.0));
+         1.2268798733741557,  -0.5578149965554813,  0.28139105017721583,
+        -0.04057576262431372,  1.1122868293970594, -0.07171106666151701,
+        -0.07637294974672142, -0.4214933239627914,  1.5869240244272418);
+    return LMS * M;
 }
 
-float f2(float x, float delta) {
-    return x > delta ?
-        pow(x, 3.0) :
-        (x - deltac) * (3.0 * pow(delta, 2.0));
-}
-
-vec3 XYZ_ref = RGB_to_XYZ(vec3(L_sdr));
-
-vec3 XYZ_to_Lab(vec3 XYZ) {
-    float X = XYZ.x;
-    float Y = XYZ.y;
-    float Z = XYZ.z;
-
-    X = f1(X / XYZ_ref.x, delta);
-    Y = f1(Y / XYZ_ref.y, delta);
-    Z = f1(Z / XYZ_ref.z, delta);
-
-    float L = 116.0 * Y - 16.0;
-    float a = 500.0 * (X - Y);
-    float b = 200.0 * (Y - Z);
+vec3 LMS_to_Lab(vec3 LMS) {
+    mat3 M = mat3(
+        0.2104542553,  0.7936177850, -0.0040720468,
+        1.9779984951, -2.4285922050,  0.4505937099,
+        0.0259040371,  0.7827717662, -0.8086757660);
+
+    LMS = vec3(
+        cbrt(LMS.x),
+        cbrt(LMS.y),
+        cbrt(LMS.z)
+    );
 
-    return vec3(L, a, b);
+    return LMS * M;
 }
 
-vec3 Lab_to_XYZ(vec3 Lab) {
-    float L = Lab.x;
-    float a = Lab.y;
-    float b = Lab.z;
-
-    float Y = (L + 16.0) / 116.0;
-    float X = Y + a / 500.0;
-    float Z = Y - b / 200.0;
+vec3 Lab_to_LMS(vec3 Lab) {
+    mat3 M = mat3(
+        0.99999999845051981432,  0.39633779217376785678,   0.21580375806075880339,
+        1.0000000088817607767,  -0.1055613423236563494,   -0.063854174771705903402,
+        1.0000000546724109177,  -0.089484182094965759684, -1.2914855378640917399);
 
-    X = f2(X, delta) * XYZ_ref.x;
-    Y = f2(Y, delta) * XYZ_ref.y;
-    Z = f2(Z, delta) * XYZ_ref.z;
+    Lab = Lab * M;
 
-    return vec3(X, Y, Z);
+    return vec3(
+        pow(Lab.x, 3.0),
+        pow(Lab.y, 3.0),
+        pow(Lab.z, 3.0)
+    );
 }
 
 vec3 RGB_to_Lab(vec3 color) {
-    color *= L_sdr;
     color  = RGB_to_XYZ(color);
-    color  = XYZD65_to_XYZD50(color);
-    color  = XYZ_to_Lab(color);
+    color  = XYZ_to_LMS(color);
+    color  = LMS_to_Lab(color);
     return color;
 }
 
 vec3 Lab_to_RGB(vec3 color) {
-    color  = Lab_to_XYZ(color);
-    color  = XYZD50_to_XYZD65(color);
+    color  = Lab_to_LMS(color);
+    color  = LMS_to_XYZ(color);
     color  = XYZ_to_RGB(color);
-    color /= L_sdr;
     return color;
 }