sketch PQ and HLG APIs

In the end I'd like to support 3 new forms to transfer function,

   - PQ / inv PQ
   - HLG
   - inv HLG

PQ's inverse is in the same form, but no luck for HLG.

I've written the HLG and inv HLG factories to take the same parameters,
hoping that makes things easier to work with.  I figure most people come
at this from the perspective of the encoding inverse HLG (hence the
"log" in the name... it's not HEG hybrid-exponent-gamma), so those
shared parameters are in the form most natural to describe the inverse,
making the description of the decoding HLG forward direction just a bit
messier.

It's unclear to me if we really want to be the arbiters of canonical
"PQ" or "HLG" functions or if we just provide the PQ-ish and HLG-ish
factories.  There's so much nuance and opinion involved in what the
output range of these functions should be.   PQ 0->1 or 0->10000?  HLG
0->1 or 0->12?  Or some other range!?

I'm also uncertain whether HLG and its inverse will need a parameter
to control the gamma/log crossover point, or if it's always best at 1.
Will probably find this out as I implement and test.

I learned even the beautiful PQ function needs a clamp up to zero
in its numerator or the math will go complex.  Crushed my heart.

So far I'm roughly at:
   [x] draft some APIs
   [x] impl. and test scalar PQ
   [ ] impl. and test scalar HLG
   [ ] impl. and test vector PQ
   [ ] impl. and test vector HLG
   [ ] make sure we're happy with APIs

I think this might be a good first spot to land now that
the unit tests do their first vaguely interesting work.

Bug: chromium:960620
Change-Id: If22584f8ad24158fdaf3c30e084dd6754455d8fb
Reviewed-on: https://skia-review.googlesource.com/c/skcms/+/246807
Commit-Queue: Mike Klein <mtklein@google.com>
Reviewed-by: Brian Osman <brianosman@google.com>
diff --git a/skcms.cc b/skcms.cc
index 2a86753..587671b 100644
--- a/skcms.cc
+++ b/skcms.cc
@@ -61,6 +61,80 @@
     return x;
 }
 
+// Most transfer functions we work with are sRGBish.
+// For exotic HDR transfer functions, we encode them using a tf.g that makes no sense,
+// and repurpose the other fields to hold the parameters of the HDR functions.
+enum TFKind { Bad, sRGBish, PQish, HLGish, HLGinvish };
+struct TF_PQish  { float A,B,C,D,E,F; };
+struct TF_HLGish { float R,G,a,b,c; };
+
+static float TFKind_marker(TFKind kind) {
+    // We'd use different NaNs, but those aren't guaranteed to be preserved by WASM.
+    return -(float)kind;
+}
+
+static TFKind classify(const skcms_TransferFunction& tf, TF_PQish*   pq = nullptr
+                                                       , TF_HLGish* hlg = nullptr) {
+    if (tf.g < 0 && (int)tf.g == tf.g) {
+        // TODO: sanity checks for PQ/HLG like we do for sRGBish.
+        switch (-(int)tf.g) {
+            case PQish:     if (pq ) { memcpy(pq , &tf.a, sizeof(*pq )); } return PQish;
+            case HLGish:    if (hlg) { memcpy(hlg, &tf.a, sizeof(*hlg)); } return HLGish;
+            case HLGinvish: if (hlg) { memcpy(hlg, &tf.a, sizeof(*hlg)); } return HLGinvish;
+        }
+        return Bad;
+    }
+
+    // Basic sanity checks for sRGBish transfer functions.
+    if (isfinitef_(tf.a + tf.b + tf.c + tf.d + tf.e + tf.f + tf.g)
+            // a,c,d,g should be non-negative to make any sense.
+            && tf.a >= 0
+            && tf.c >= 0
+            && tf.d >= 0
+            && tf.g >= 0
+            // Raising a negative value to a fractional tf->g produces complex numbers.
+            && tf.a * tf.d + tf.b >= 0) {
+        return sRGBish;
+    }
+
+    return Bad;
+}
+
+// TODO: temporary shim for old call sites
+static bool tf_is_valid(const skcms_TransferFunction* tf) {
+    return classify(*tf) == sRGBish;
+}
+
+
+bool skcms_TransferFunction_makePQish(skcms_TransferFunction* tf,
+                                      float A, float B, float C,
+                                      float D, float E, float F) {
+    *tf = { TFKind_marker(PQish), A,B,C,D,E,F };
+    assert(classify(*tf) == PQish);
+    return true;
+}
+
+float skcms_TransferFunction_eval(const skcms_TransferFunction* tf, float x) {
+    float sign = x < 0 ? -1.0f : 1.0f;
+    x *= sign;
+
+    TF_PQish pq;
+    switch (classify(*tf, &pq)) {
+        case Bad:       break;
+        case HLGish:    break;
+        case HLGinvish: break;
+
+        case sRGBish: return sign * (x < tf->d ?       tf->c * x + tf->f
+                                               : powf_(tf->a * x + tf->b, tf->g) + tf->e);
+
+        case PQish: return sign * powf_(fmaxf_(pq.A + pq.B * powf_(x, pq.C), 0)
+                                            / (pq.D + pq.E * powf_(x, pq.C)),
+                                        pq.F);
+    }
+    return 0;
+}
+
+
 static float eval_curve(const skcms_Curve* curve, float x) {
     if (curve->table_entries == 0) {
         return skcms_TransferFunction_eval(&curve->parametric, x);
@@ -255,25 +329,6 @@
            read_tag_xyz(bXYZ, &toXYZ->vals[0][2], &toXYZ->vals[1][2], &toXYZ->vals[2][2]);
 }
 
-static bool tf_is_valid(const skcms_TransferFunction* tf) {
-    // Reject obviously malformed inputs
-    if (!isfinitef_(tf->a + tf->b + tf->c + tf->d + tf->e + tf->f + tf->g)) {
-        return false;
-    }
-
-    // All of these parameters should be non-negative
-    if (tf->a < 0 || tf->c < 0 || tf->d < 0 || tf->g < 0) {
-        return false;
-    }
-
-    // It's rather _complex_ to raise a negative number to a fractional power tf->g.
-    if (tf->a * tf->d + tf->b < 0) {
-        return false;
-    }
-
-    return true;
-}
-
 typedef struct {
     uint8_t type          [4];
     uint8_t reserved_a    [4];
@@ -1419,14 +1474,6 @@
                                 : exp2f_(log2f_(x) * y);
 }
 
-float skcms_TransferFunction_eval(const skcms_TransferFunction* tf, float x) {
-    float sign = x < 0 ? -1.0f : 1.0f;
-    x *= sign;
-
-    return sign * (x < tf->d ? tf->c * x + tf->f
-                             : powf_(tf->a * x + tf->b, tf->g) + tf->e);
-}
-
 #if defined(__clang__)
     [[clang::no_sanitize("float-divide-by-zero")]]  // Checked for by tf_is_valid() on the way out.
 #endif
diff --git a/skcms.h b/skcms.h
index c316385..bd3a96e 100644
--- a/skcms.h
+++ b/skcms.h
@@ -51,6 +51,50 @@
 SKCMS_API bool  skcms_TransferFunction_invert(const skcms_TransferFunction*,
                                               skcms_TransferFunction*);
 
+// We can jam a couple alternate transfer function forms into skcms_TransferFunction,
+// including those matching the general form of the SMPTE ST 2084 PQ function and its inverse:
+//
+//                       max(A + B|x|^C, 0)
+//    tf(x) = sign(x) * (------------------) ^ F
+   //                       (D + E|x|^C)
+SKCMS_API bool skcms_TransferFunction_makePQish(skcms_TransferFunction*,
+                                                float A, float B, float C,
+                                                float D, float E, float F);
+
+static inline bool skcms_TransferFunction_makePQ(skcms_TransferFunction* tf) {
+    return skcms_TransferFunction_makePQish(tf, -107/128.0f,         1.0f,   32/2523.0f
+                                              , 2413/128.0f, -2392/128.0f, 8192/1305.0f);
+}
+static inline bool skcms_TransferFunction_makePQinv(skcms_TransferFunction* tf) {
+    return skcms_TransferFunction_makePQish(tf, 107/128.0f, 2413/128.0f, 1305/8192.0f
+                                              ,       1.0f, 2392/128.0f, 2523/  32.0f);
+}
+
+// skcms_TransferFunction also supports functions of the form of HLG's inverse...
+//
+//             { sign(linear) * ( R|linear|^G          ) when 0 <= |linear| <= 1
+//   encoded = { sign(linear) * ( a ln(|linear|-b) + c ) when 1 <  |linear|
+SKCMS_API bool skcms_TransferFunction_makeHLGinvish(skcms_TransferFunction*,
+                                                    float R, float G,
+                                                    float a, float b, float c);
+
+// ... and of the form of HLG itself, factored to take the same five parameters.
+//
+//            { sign(encoded) * ( (|encoded|/R)^(1/G) )        when 0 <= |encoded| <= R
+//   linear = { sign(encoded) * ( e^( (|encoded|-c)/a ) + b )  when R <  |encoded|
+SKCMS_API bool skcms_TransferFunction_makeHLGish(skcms_TransferFunction*,
+                                                float R, float G,
+                                                float a, float b, float c);
+
+static inline bool skcms_TransferFunction_makeHLG(skcms_TransferFunction* tf) {
+    return skcms_TransferFunction_makeHLGish(tf,
+            0.5f, 0.5f , 0.17883277f, 0.28466892f, 0.55991073f);
+}
+static inline bool skcms_TransferFunction_makeHLGinv(skcms_TransferFunction* tf) {
+    return skcms_TransferFunction_makeHLGinvish(tf,
+            0.5f, 0.5f , 0.17883277f, 0.28466892f, 0.55991073f);
+}
+
 // Unified representation of 'curv' or 'para' tag data, or a 1D table from 'mft1' or 'mft2'
 typedef union skcms_Curve {
     struct {
diff --git a/tests.c b/tests.c
index 6911b0b..5fa964d 100644
--- a/tests.c
+++ b/tests.c
@@ -1373,6 +1373,51 @@
   //expect(0 == memcmp(sRGB, &sRGB2, sizeof(skcms_TransferFunction)));
 }
 
+static void test_PQ() {
+    {
+        // This PQ function maps [0,1] to [0,1].
+        skcms_TransferFunction pq;
+        expect(skcms_TransferFunction_makePQ(&pq));
+
+        expect(0.0000f == skcms_TransferFunction_eval(&pq, 0.0f));
+        expect(1.0000f == skcms_TransferFunction_eval(&pq, 1.0f));
+
+        // 100 nits is around 0.508.
+        expect(0.0099f < skcms_TransferFunction_eval(&pq, 0.508f));
+        expect(0.0101f > skcms_TransferFunction_eval(&pq, 0.508f));
+    }
+
+    {
+        // Let's see if we can get absolute 0-10000 nits.
+        skcms_TransferFunction pq_abs;
+
+        // Mathematically to get 10000 on the output, we want to
+        // scale the A and B PQ terms by R = 10000 ^ (1/F).
+        float R = powf_(10000.0f, 1305/8192.0f);   // ~= 4.33691
+        expect(skcms_TransferFunction_makePQish(&pq_abs,
+                    R*(-107/128.0f), R*       1.0f,   32/2523.0f,
+                       2413/128.0f,   -2392/128.0f, 8192/1305.0f));
+
+        // That gets us close.
+        expect(0.0f == skcms_TransferFunction_eval(&pq_abs, 0.0f));
+        expect(   99.8f < skcms_TransferFunction_eval(&pq_abs, 0.508f));
+        expect(  100.0f > skcms_TransferFunction_eval(&pq_abs, 0.508f));
+        expect( 9989.0f < skcms_TransferFunction_eval(&pq_abs, 1.0f));
+        expect( 9991.0f > skcms_TransferFunction_eval(&pq_abs, 1.0f));
+
+        // We can get a lot closer with an unprincpled tweak to that math.
+        R = powf_(10009.9f, 1305/8192.0f);  // ~= 4.33759
+        expect(skcms_TransferFunction_makePQish(&pq_abs,
+                    R*(-107/128.0f), R*       1.0f,   32/2523.0f,
+                       2413/128.0f,   -2392/128.0f, 8192/1305.0f));
+        expect(0.0f == skcms_TransferFunction_eval(&pq_abs, 0.0f));
+        expect(   99.9f < skcms_TransferFunction_eval(&pq_abs, 0.508f));
+        expect(  100.0f > skcms_TransferFunction_eval(&pq_abs, 0.508f));
+        expect( 9999.0f < skcms_TransferFunction_eval(&pq_abs, 1.0f));
+        expect(10000.0f > skcms_TransferFunction_eval(&pq_abs, 1.0f));
+    }
+}
+
 int main(int argc, char** argv) {
     bool regenTestData = false;
     for (int i = 1; i < argc; ++i) {
@@ -1407,6 +1452,7 @@
     test_TF_invert();
     test_Clamp();
     test_Premul();
+    test_PQ();
 
     // Temporarily disable some tests while getting FP16 compute working.
     if (!kFP16) {