[babl] babl: code handling ICC v4, CIE TRC formulas

[Øyvind "pippin" Kolås <ok src gnome org>]

commit f25c510ae26866b2a7261428b373e38ffaa8dda5
Author: Øyvind Kolås <pippin gimp org>
Date:   Fri Aug 7 23:12:15 2020 +0200

    babl: code handling ICC v4, CIE TRC formulas

 babl/babl-icc.c      | 24 ++++++++++++---
 babl/babl-internal.h |  2 ++
 babl/babl-trc.c      | 84 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 babl/babl-trc.h      |  3 +-
 4 files changed, 108 insertions(+), 5 deletions(-)
---
diff --git a/babl/babl-icc.c b/babl/babl-icc.c
index aa5a3862d..7d1fd6705 100644
--- a/babl/babl-icc.c
+++ b/babl/babl-icc.c
@@ -384,9 +384,27 @@ babl_trc_from_icc (ICC         *state,
               g = icc_read (s15f16, offset + 12 + 4 * 0);
               return babl_trc_gamma (g);
               break;
+            case 1:
+              {
+                float a,b,c;
+                g = icc_read (s15f16, offset + 12 + 4 * 0);
+                a = icc_read (s15f16, offset + 12 + 4 * 1);
+                b = icc_read (s15f16, offset + 12 + 4 * 2);
+                c = 0;
+                return babl_trc_formula_cie (g, a, b, c);
+              }
+            case 2:
+              {
+                float a,b,c;
+                g = icc_read (s15f16, offset + 12 + 4 * 0);
+                a = icc_read (s15f16, offset + 12 + 4 * 1);
+                b = icc_read (s15f16, offset + 12 + 4 * 2);
+                c = icc_read (s15f16, offset + 12 + 4 * 3);
+                return babl_trc_formula_cie (g, a, b, c);
+              }
             case 3:
               {
-                float a,b,c,d, e, f;
+                float a,b,c,d,e,f;
                 g = icc_read (s15f16, offset + 12 + 4 * 0);
                 a = icc_read (s15f16, offset + 12 + 4 * 1);
                 b = icc_read (s15f16, offset + 12 + 4 * 2);
@@ -408,8 +426,6 @@ babl_trc_from_icc (ICC         *state,
                 f = icc_read (s15f16, offset + 12 + 4 * 6);
                 return babl_trc_formula_srgb (g, a, b, c, d, e, f);
               }
-            case 1: // NYI
-            case 2: // NYI - can share code, like srgb formulas
             default:
               *error = "unhandled parametric TRC";
               fprintf (stderr, "unhandled parametric TRC type %i\n", function_type);
@@ -534,6 +550,7 @@ switch (trc->type)
       break;
     }
   case BABL_TRC_FORMULA_SRGB:
+  case BABL_TRC_FORMULA_CIE:
     {
       int lut_size = 512;
       if (flags == BABL_ICC_COMPACT_TRC_LUT)
@@ -671,7 +688,6 @@ babl_space_to_icc_rgb (const Babl  *babl,
         icc_write (u8, state->o + 12 + i, description[i]);
     }
 
-
     icc_write (u32, 0, state->no + 0);
     length = state->no + 0;
   }
diff --git a/babl/babl-internal.h b/babl/babl-internal.h
index 1efdb0f2f..56e95e4a2 100644
--- a/babl/babl-internal.h
+++ b/babl/babl-internal.h
@@ -386,6 +386,8 @@ babl_conversion_create_name (Babl *source, Babl *destination, int type,
 void _babl_space_add_universal_rgb (const Babl *space);
 const Babl *
 babl_trc_formula_srgb (double gamma, double a, double b, double c, double d, double e, double f);
+const Babl *
+babl_trc_formula_cie (double gamma, double a, double b, double c);
 
 
 const Babl *babl_space_match_trc_matrix (const Babl *trc_red,
diff --git a/babl/babl-trc.c b/babl/babl-trc.c
index b1ba8a81b..d8a69dd21 100644
--- a/babl/babl-trc.c
+++ b/babl/babl-trc.c
@@ -200,6 +200,44 @@ _babl_trc_formula_srgb_to_linear (const Babl *trc_,
   }
   return c * x + f;
 }
+static inline float 
+_babl_trc_formula_cie_from_linear (const Babl *trc_, 
+                                   float       value)
+{
+  BablTRC *trc = (void*)trc_;
+  float x= value;
+  float a = trc->lut[1];
+  float b = trc->lut[2];
+  float c = trc->lut[3];
+
+  if (x > c)
+  {
+    float v = _babl_trc_gamma_from_linear ((Babl *) trc, x - c);
+    v = (v-b)/a;
+    if (v < 0.0 || v >= 0.0)
+      return v;
+  }
+  return 0.0;
+}
+
+static inline float 
+_babl_trc_formula_cie_to_linear (const Babl *trc_, 
+                                 float       value)
+{
+  BablTRC *trc = (void*)trc_;
+  float x= value;
+  float a = trc->lut[1];
+  float b = trc->lut[2];
+  float c = trc->lut[3];
+
+  if (x >= -b / a)
+  {
+    return _babl_trc_gamma_to_linear ((Babl *) trc, a * x + b) + c;
+  }
+  return c;
+}
+
+
 
 static inline float 
 _babl_trc_srgb_to_linear (const Babl *trc_, 
@@ -421,6 +459,33 @@ babl_trc_new (const char *name,
                                          trc_db[i].poly_gamma_from_linear_x1,
                                          POLY_GAMMA_DEGREE, POLY_GAMMA_SCALE);
       break;
+    case BABL_TRC_FORMULA_CIE:
+      trc_db[i].lut = babl_calloc (sizeof (float), 4);
+      {
+        int j;
+        for (j = 0; j < 4; j++)
+          trc_db[i].lut[j] = lut[j];
+      }
+      trc_db[i].fun_to_linear = _babl_trc_formula_cie_to_linear;
+      trc_db[i].fun_from_linear = _babl_trc_formula_cie_from_linear;
+
+      trc_db[i].poly_gamma_to_linear_x0 = lut[4];
+      trc_db[i].poly_gamma_to_linear_x1 = POLY_GAMMA_X1;
+      babl_polynomial_approximate_gamma (&trc_db[i].poly_gamma_to_linear,
+                                         trc_db[i].gamma,
+                                         trc_db[i].poly_gamma_to_linear_x0,
+                                         trc_db[i].poly_gamma_to_linear_x1,
+                                         POLY_GAMMA_DEGREE, POLY_GAMMA_SCALE);
+
+      trc_db[i].poly_gamma_from_linear_x0 = lut[3] * lut[4];
+      trc_db[i].poly_gamma_from_linear_x1 = POLY_GAMMA_X1;
+      babl_polynomial_approximate_gamma (&trc_db[i].poly_gamma_from_linear,
+                                         trc_db[i].rgamma,
+                                         trc_db[i].poly_gamma_from_linear_x0,
+                                         trc_db[i].poly_gamma_from_linear_x1,
+                                         POLY_GAMMA_DEGREE, POLY_GAMMA_SCALE);
+      break;
+
     case BABL_TRC_FORMULA_SRGB:
       trc_db[i].lut = babl_calloc (sizeof (float), 7);
       {
@@ -510,6 +575,25 @@ babl_trc_formula_srgb (double g,
   return babl_trc_new (name, BABL_TRC_FORMULA_SRGB, g, 0, params);
 }
 
+const Babl *
+babl_trc_formula_cie (double g, 
+                      double a, 
+                      double b, 
+                      double c)
+{
+  char name[128];
+  int i;
+  float params[4]={g, a, b, c};
+
+  snprintf (name, sizeof (name), "%.6f %.6f %.4f %.4f", g, a, b, c);
+  for (i = 0; name[i]; i++)
+    if (name[i] == ',') name[i] = '.';
+  while (name[strlen(name)-1]=='0')
+    name[strlen(name)-1]='\0';
+  return babl_trc_new (name, BABL_TRC_FORMULA_CIE, g, 0, params);
+}
+
+
 const Babl *
 babl_trc_gamma (double gamma)
 {
diff --git a/babl/babl-trc.h b/babl/babl-trc.h
index 380a73610..cfe742c60 100644
--- a/babl/babl-trc.h
+++ b/babl/babl-trc.h
@@ -30,7 +30,8 @@ typedef enum {BABL_TRC_LINEAR,
               BABL_TRC_FORMULA_GAMMA,
               BABL_TRC_SRGB,
               BABL_TRC_FORMULA_SRGB,
-              BABL_TRC_LUT}
+              BABL_TRC_LUT,
+              BABL_TRC_FORMULA_CIE}
 BablTRCType;
 
 typedef struct