Initial implementation of a SkColorSpace_A2B xform

src/core/SkColorSpaceXformPriv.h

Index: src/core/SkColorSpaceXformPriv.h
diff --git a/src/core/SkColorSpaceXformPriv.h b/src/core/SkColorSpaceXformPriv.h
new file mode 100644
index 0000000000000000000000000000000000000000..6938198a83ac14a7a52af5b60abbcfd5211008eb
--- /dev/null
+++ b/src/core/SkColorSpaceXformPriv.h
@@ -0,0 +1,354 @@
+/*
+ * Copyright 2016 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#ifndef SkColorSpaceXformPriv_DEFINED
+#define SkColorSpaceXformPriv_DEFINED
+
+#include "SkColorSpace_Base.h"
+#include "SkHalf.h"
+#include "SkSRGB.h"
+
+#define AI SK_ALWAYS_INLINE
+
+#define SkCSXformPrintfDefined 0

[msarett, 2016/10/26 20:43:29]

What's this for?

[raftias, 2016/10/27 18:39:38]

To block out chunks of code (variable declarations, etc) that exist just for
printing.

+#define SkCSXformPrintf(...)
+
+static AI void interp_3d_clut(float dst[3], float src[3], const SkColorLookUpTable* colorLUT) {

[msarett, 2016/10/26 20:43:29]

Can this be moved to SkColorSpaceXform_A2B.cpp?

[raftias, 2016/10/27 18:39:38]

Done.

+    // Call the src components x, y, and z.
+    uint8_t maxX = colorLUT->fGridPoints[0] - 1;
+    uint8_t maxY = colorLUT->fGridPoints[1] - 1;
+    uint8_t maxZ = colorLUT->fGridPoints[2] - 1;
+
+    // An approximate index into each of the three dimensions of the table.
+    float x = src[0] * maxX;
+    float y = src[1] * maxY;
+    float z = src[2] * maxZ;
+
+    // This gives us the low index for our interpolation.
+    int ix = sk_float_floor2int(x);
+    int iy = sk_float_floor2int(y);
+    int iz = sk_float_floor2int(z);
+
+    // Make sure the low index is not also the max index.
+    ix = (maxX == ix) ? ix - 1 : ix;
+    iy = (maxY == iy) ? iy - 1 : iy;
+    iz = (maxZ == iz) ? iz - 1 : iz;
+
+    // Weighting factors for the interpolation.
+    float diffX = x - ix;
+    float diffY = y - iy;
+    float diffZ = z - iz;
+
+    // Constants to help us navigate the 3D table.
+    // Ex: Assume x = a, y = b, z = c.
+    //     table[a * n001 + b * n010 + c * n100] logically equals table[a][b][c].
+    const int n000 = 0;
+    const int n001 = 3 * colorLUT->fGridPoints[1] * colorLUT->fGridPoints[2];
+    const int n010 = 3 * colorLUT->fGridPoints[2];
+    const int n011 = n001 + n010;
+    const int n100 = 3;
+    const int n101 = n100 + n001;
+    const int n110 = n100 + n010;
+    const int n111 = n110 + n001;
+
+    // Base ptr into the table.
+    const float* ptr = &(colorLUT->table()[ix*n001 + iy*n010 + iz*n100]);
+
+    // The code below performs a tetrahedral interpolation for each of the three
+    // dst components.  Once the tetrahedron containing the interpolation point is
+    // identified, the interpolation is a weighted sum of grid values at the
+    // vertices of the tetrahedron.  The claim is that tetrahedral interpolation
+    // provides a more accurate color conversion.
+    // blogs.mathworks.com/steve/2006/11/24/tetrahedral-interpolation-for-colorspace-conversion/
+    //
+    // I have one test image, and visually I can't tell the difference between
+    // tetrahedral and trilinear interpolation.  In terms of computation, the
+    // tetrahedral code requires more branches but less computation.  The
+    // SampleICC library provides an option for the client to choose either
+    // tetrahedral or trilinear.
+    for (int i = 0; i < 3; i++) {
+        if (diffZ < diffY) {
+            if (diffZ < diffX) {
+                dst[i] = (ptr[n000] + diffZ * (ptr[n110] - ptr[n010]) +
+                                      diffY * (ptr[n010] - ptr[n000]) +
+                                      diffX * (ptr[n111] - ptr[n110]));
+            } else if (diffY < diffX) {
+                dst[i] = (ptr[n000] + diffZ * (ptr[n111] - ptr[n011]) +
+                                      diffY * (ptr[n011] - ptr[n001]) +
+                                      diffX * (ptr[n001] - ptr[n000]));
+            } else {
+                dst[i] = (ptr[n000] + diffZ * (ptr[n111] - ptr[n011]) +
+                                      diffY * (ptr[n010] - ptr[n000]) +
+                                      diffX * (ptr[n011] - ptr[n010]));
+            }
+        } else {
+            if (diffZ < diffX) {
+                dst[i] = (ptr[n000] + diffZ * (ptr[n101] - ptr[n001]) +
+                                      diffY * (ptr[n111] - ptr[n101]) +
+                                      diffX * (ptr[n001] - ptr[n000]));
+            } else if (diffY < diffX) {
+                dst[i] = (ptr[n000] + diffZ * (ptr[n100] - ptr[n000]) +
+                                      diffY * (ptr[n111] - ptr[n101]) +
+                                      diffX * (ptr[n101] - ptr[n100]));
+            } else {
+                dst[i] = (ptr[n000] + diffZ * (ptr[n100] - ptr[n000]) +
+                                      diffY * (ptr[n110] - ptr[n100]) +
+                                      diffX * (ptr[n111] - ptr[n110]));
+            }
+        }
+
+        // Increment the table ptr in order to handle the next component.
+        // Note that this is the how table is designed: all of nXXX
+        // variables are multiples of 3 because there are 3 output
+        // components.
+        ptr++;
+    }
+}
+
+// Interpolating lookup in a variably sized table.
+static AI float interp_lut(float input, const float* table, int tableSize) {
+    float index = input * (tableSize - 1);
+    float diff = index - sk_float_floor2int(index);
+    return table[(int) sk_float_floor2int(index)] * (1.0f - diff) +
+           table[(int) sk_float_ceil2int(index)] * diff;
+}
+
+
+static AI void transform_gamut(const Sk4f& r, const Sk4f& g, const Sk4f& b, const Sk4f& a,
+                                   const Sk4f& rXgXbX, const Sk4f& rYgYbY, const Sk4f& rZgZbZ,

[msarett, 2016/10/26 20:43:29]

nit: spacing

[raftias, 2016/10/27 18:39:38]

Will fix in next patchset

+                                   Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f& da) {
+    dr = rXgXbX[0]*r + rYgYbY[0]*g + rZgZbZ[0]*b;
+    dg = rXgXbX[1]*r + rYgYbY[1]*g + rZgZbZ[1]*b;
+    db = rXgXbX[2]*r + rYgYbY[2]*g + rZgZbZ[2]*b;
+    da = a;
+}
+
+static AI void translate_gamut(const Sk4f& rTgTbT, Sk4f& dr, Sk4f& dg, Sk4f& db) {
+    dr = dr + rTgTbT[0];
+    dg = dg + rTgTbT[1];
+    db = db + rTgTbT[2];
+}
+
+static AI void load_matrix(const float matrix[16],
+                               Sk4f& rXgXbX, Sk4f& rYgYbY, Sk4f& rZgZbZ, Sk4f& rTgTbT) {

[msarett, 2016/10/26 20:43:29]

nit: spacing

+    rXgXbX = Sk4f::Load(matrix +  0);
+    rYgYbY = Sk4f::Load(matrix +  4);
+    rZgZbZ = Sk4f::Load(matrix +  8);
+    rTgTbT = Sk4f::Load(matrix + 12);
+}
+
+
+enum Order {
+    kRGBA_Order,
+    kBGRA_Order,
+};
+
+static AI void set_rb_shifts(Order kOrder, int* kRShift, int* kBShift) {
+    if (kRGBA_Order == kOrder) {
+        *kRShift = 0;
+        *kBShift = 16;
+    } else {
+        *kRShift = 16;
+        *kBShift = 0;
+    }
+}
+
+template <Order kOrder>
+static AI void load_rgb_linear(const uint32_t* src,
+                                   Sk4f& r, Sk4f& g, Sk4f& b, Sk4f& a,
+                                   const float* const[3]) {
+    int kRShift, kGShift = 8, kBShift;
+    set_rb_shifts(kOrder, &kRShift, &kBShift);
+    r = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> kRShift) & 0xFF);
+    g = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> kGShift) & 0xFF);
+    b = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> kBShift) & 0xFF);
+    a = 0.0f; // Don't let the compiler complain that |a| is uninitialized.
+}
+
+template <Order kOrder>
+static AI void load_rgba_linear(const uint32_t* src,
+                                    Sk4f& r, Sk4f& g, Sk4f& b, Sk4f& a,
+                                    const float* const[3]) {
+    int kRShift, kGShift = 8, kBShift;
+    set_rb_shifts(kOrder, &kRShift, &kBShift);
+    r = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> kRShift) & 0xFF);
+    g = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> kGShift) & 0xFF);
+    b = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> kBShift) & 0xFF);
+    a = (1.0f / 255.0f) * SkNx_cast<float>((Sk4u::Load(src) >> 24));
+}
+
+static AI void premultiply(Sk4f& dr, Sk4f& dg, Sk4f& db, const Sk4f& da) {
+    dr = da * dr;
+    dg = da * dg;
+    db = da * db;
+}
+
+template <Order kOrder>
+static AI void store_srgb(void* dst, const uint32_t* src,
+                              Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f&,
+                              const uint8_t* const[3]) {
+    int kRShift, kGShift = 8, kBShift;
+    set_rb_shifts(kOrder, &kRShift, &kBShift);
+    dr = sk_linear_to_srgb_needs_trunc(dr);
+    dg = sk_linear_to_srgb_needs_trunc(dg);
+    db = sk_linear_to_srgb_needs_trunc(db);
+
+    dr = sk_clamp_0_255(dr);
+    dg = sk_clamp_0_255(dg);
+    db = sk_clamp_0_255(db);
+
+    Sk4i da = Sk4i::Load(src) & 0xFF000000;
+
+    Sk4i rgba = (SkNx_cast<int>(dr) << kRShift)
+              | (SkNx_cast<int>(dg) << kGShift)
+              | (SkNx_cast<int>(db) << kBShift)
+              | (da                           );
+    rgba.store(dst);
+}
+
+static AI Sk4f linear_to_2dot2(const Sk4f& x) {
+    // x^(29/64) is a very good approximation of the true value, x^(1/2.2).
+    auto x2  = x.rsqrt(),                            // x^(-1/2)
+         x32 = x2.rsqrt().rsqrt().rsqrt().rsqrt(),   // x^(-1/32)
+         x64 = x32.rsqrt();                          // x^(+1/64)
+
+    // 29 = 32 - 2 - 1
+    return 255.0f * x2.invert() * x32 * x64.invert();
+}
+
+template <Order kOrder>
+static AI void store_2dot2(void* dst, const uint32_t* src,
+                               Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f&,
+                               const uint8_t* const[3]) {
+    int kRShift, kGShift = 8, kBShift;
+    set_rb_shifts(kOrder, &kRShift, &kBShift);
+    dr = linear_to_2dot2(dr);
+    dg = linear_to_2dot2(dg);
+    db = linear_to_2dot2(db);
+
+    dr = sk_clamp_0_255(dr);
+    dg = sk_clamp_0_255(dg);
+    db = sk_clamp_0_255(db);
+
+    Sk4i da = Sk4i::Load(src) & 0xFF000000;
+
+    Sk4i rgba = (Sk4f_round(dr) << kRShift)
+              | (Sk4f_round(dg) << kGShift)
+              | (Sk4f_round(db) << kBShift)
+              | (da                       );
+    rgba.store(dst);
+}
+
+template <Order kOrder>
+static AI void store_linear(void* dst, const uint32_t* src,
+                               Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f&,
+                               const uint8_t* const[3]) {
+    int kRShift, kGShift = 8, kBShift;
+    set_rb_shifts(kOrder, &kRShift, &kBShift);
+    dr = sk_clamp_0_255(255.0f * dr);
+    dg = sk_clamp_0_255(255.0f * dg);
+    db = sk_clamp_0_255(255.0f * db);
+
+    Sk4i da = Sk4i::Load(src) & 0xFF000000;
+
+    Sk4i rgba = (Sk4f_round(dr) << kRShift)
+              | (Sk4f_round(dg) << kGShift)
+              | (Sk4f_round(db) << kBShift)
+              | (da                       );
+    rgba.store(dst);
+}
+
+template <Order kOrder>
+static AI void store_f16(void* dst, const uint32_t* src,
+                             Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f& da,
+                             const uint8_t* const[3]) {
+    Sk4h::Store4(dst, SkFloatToHalf_finite_ftz(dr),
+                      SkFloatToHalf_finite_ftz(dg),
+                      SkFloatToHalf_finite_ftz(db),
+                      SkFloatToHalf_finite_ftz(da));
+}
+
+template <Order kOrder>
+static AI void store_f32(void* dst, const uint32_t* src,
+                             Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f& da,
+                             const uint8_t* const[3]) {
+    Sk4f::Store4(dst, dr, dg, db, da);
+}
+
+template <Order kOrder>
+static AI void store_f16_opaque(void* dst, const uint32_t* src,
+                                    Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f&,
+                                    const uint8_t* const[3]) {
+    Sk4h::Store4(dst, SkFloatToHalf_finite_ftz(dr),
+                      SkFloatToHalf_finite_ftz(dg),
+                      SkFloatToHalf_finite_ftz(db),
+                      SK_Half1);
+}
+
+template <Order kOrder>
+static AI void store_generic(void* dst, const uint32_t* src,
+                                 Sk4f& dr, Sk4f& dg, Sk4f& db, Sk4f&,
+                                 const uint8_t* const dstTables[3]) {
+    int kRShift, kGShift = 8, kBShift;
+    set_rb_shifts(kOrder, &kRShift, &kBShift);
+    dr = Sk4f::Min(Sk4f::Max(1023.0f * dr, 0.0f), 1023.0f);
+    dg = Sk4f::Min(Sk4f::Max(1023.0f * dg, 0.0f), 1023.0f);
+    db = Sk4f::Min(Sk4f::Max(1023.0f * db, 0.0f), 1023.0f);
+
+    Sk4i ir = Sk4f_round(dr);
+    Sk4i ig = Sk4f_round(dg);
+    Sk4i ib = Sk4f_round(db);
+
+    Sk4i da = Sk4i::Load(src) & 0xFF000000;
+
+    uint32_t* dst32 = (uint32_t*) dst;
+    dst32[0] = dstTables[0][ir[0]] << kRShift
+             | dstTables[1][ig[0]] << kGShift
+             | dstTables[2][ib[0]] << kBShift
+             | da[0];
+    dst32[1] = dstTables[0][ir[1]] << kRShift
+             | dstTables[1][ig[1]] << kGShift
+             | dstTables[2][ib[1]] << kBShift
+             | da[1];
+    dst32[2] = dstTables[0][ir[2]] << kRShift
+             | dstTables[1][ig[2]] << kGShift
+             | dstTables[2][ib[2]] << kBShift
+             | da[2];
+    dst32[3] = dstTables[0][ir[3]] << kRShift
+             | dstTables[1][ig[3]] << kGShift
+             | dstTables[2][ib[3]] << kBShift
+             | da[3];
+}
+
+typedef decltype(load_rgb_linear<kRGBA_Order>  )* LoadFn;
+typedef decltype(store_linear<kRGBA_Order>     )* StoreFn;
+
+static AI int num_tables(SkColorSpace_XYZ* space) {
+    switch (space->gammaNamed()) {
+        case kSRGB_SkGammaNamed:
+        case k2Dot2Curve_SkGammaNamed:
+        case kLinear_SkGammaNamed:
+            return 0;
+        default: {
+            const SkGammas* gammas = space->gammas();
+            SkASSERT(gammas);
+
+            bool gammasAreMatching = (gammas->type(0) == gammas->type(1)) &&
+                                     (gammas->data(0) == gammas->data(1)) &&
+                                     (gammas->type(0) == gammas->type(2)) &&
+                                     (gammas->data(0) == gammas->data(2));
+
+            // It's likely that each component will have the same gamma.  In this case,
+            // we only need to build one table.
+            return gammasAreMatching ? 1 : 3;
+        }
+    }
+}
+
+#undef AI
+
+#endif