Add support for 3D colorLUTs to SkColorXform

src/core/SkColorSpaceXform.cpp

 /*
  * Copyright 2016 Google Inc.
  *
  * Use of this source code is governed by a BSD-style license that can be
  * found in the LICENSE file.
  */
 
 #include "SkColorPriv.h"
 #include "SkColorSpace_Base.h"
 #include "SkColorSpaceXform.h"
 #include "SkOpts.h"
 
 static inline bool compute_gamut_xform(SkMatrix44* srcToDst, const SkMatrix44& srcToXYZ,
                                        const SkMatrix44& dstToXYZ) {
     if (!dstToXYZ.invert(srcToDst)) {
         return false;
     }
 
     srcToDst->postConcat(srcToXYZ);
     return true;
 }
 
 std::unique_ptr<SkColorSpaceXform> SkColorSpaceXform::New(const sk_sp<SkColorSpace>& srcSpace,
                                                           const sk_sp<SkColorSpace>& dstSpace) {
     if (!srcSpace || !dstSpace) {
         // Invalid input
         return nullptr;
     }
 
-    if (as_CSB(srcSpace)->colorLUT() || as_CSB(dstSpace)->colorLUT()) {
-        // Unimplemented
+    if (as_CSB(dstSpace)->colorLUT()) {
+        // It would be really weird for a dst profile to have a color LUT.  I don't think
+        // we need to support this.
         return nullptr;
     }
 
     SkMatrix44 srcToDst(SkMatrix44::kUninitialized_Constructor);
     if (!compute_gamut_xform(&srcToDst, srcSpace->xyz(), dstSpace->xyz())) {
         return nullptr;
     }
 
     if (0.0f == srcToDst.getFloat(3, 0) &&
         0.0f == srcToDst.getFloat(3, 1) &&
-        0.0f == srcToDst.getFloat(3, 2))
+        0.0f == srcToDst.getFloat(3, 2) &&
+        !as_CSB(srcSpace)->colorLUT())
     {
         switch (srcSpace->gammaNamed()) {
             case SkColorSpace::kSRGB_GammaNamed:
                 if (SkColorSpace::kSRGB_GammaNamed == dstSpace->gammaNamed()) {
                     return std::unique_ptr<SkColorSpaceXform>(
                             new SkFastXform<SkColorSpace::kSRGB_GammaNamed,
                                             SkColorSpace::kSRGB_GammaNamed>(srcToDst));
                 } else if (SkColorSpace::k2Dot2Curve_GammaNamed == dstSpace->gammaNamed()) {
                     return std::unique_ptr<SkColorSpaceXform>(
                             new SkFastXform<SkColorSpace::kSRGB_GammaNamed,
                                             SkColorSpace::k2Dot2Curve_GammaNamed>(srcToDst));
                 }
                 break;
             case SkColorSpace::k2Dot2Curve_GammaNamed:
                 if (SkColorSpace::kSRGB_GammaNamed == dstSpace->gammaNamed()) {
                     return std::unique_ptr<SkColorSpaceXform>(
                             new SkFastXform<SkColorSpace::k2Dot2Curve_GammaNamed,
                                             SkColorSpace::kSRGB_GammaNamed>(srcToDst));
                 } else if (SkColorSpace::k2Dot2Curve_GammaNamed == dstSpace->gammaNamed()) {
                     return std::unique_ptr<SkColorSpaceXform>(
                             new SkFastXform<SkColorSpace::k2Dot2Curve_GammaNamed,
                                             SkColorSpace::k2Dot2Curve_GammaNamed>(srcToDst));
                 }
                 break;
             default:
                 break;
         }
     }
 
-    return std::unique_ptr<SkColorSpaceXform>(
-            new SkDefaultXform(srcSpace, srcToDst, dstSpace));
+    return std::unique_ptr<SkColorSpaceXform>(new SkDefaultXform(srcSpace, srcToDst, dstSpace));
 }
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 
 static void build_src_to_dst(float srcToDstArray[12], const SkMatrix44& srcToDstMatrix) {
     // Build the following row major matrix:
     //   rX gX bX 0
     //   rY gY bY 0
     //   rZ gZ bZ 0
     // Swap R and B if necessary to make sure that we output SkPMColor order.
@@ skipping 449 matching lines @@
                                         float b, float c, float d, float e, float f) {
     for (int i = 0; i < outTableSize; i++) {
         float x = ((float) i) * (1.0f / ((float) (outTableSize - 1)));
         float y = inverse_parametric(x, g, a, b, c, d, e, f);
         outTable[i] = clamp_normalized_float_to_byte(y);
     }
 }
 
 SkDefaultXform::SkDefaultXform(const sk_sp<SkColorSpace>& srcSpace, const SkMatrix44& srcToDst,
                                const sk_sp<SkColorSpace>& dstSpace)
-    : fSrcToDst(srcToDst)
+    : fColorLUT(sk_ref_sp((SkColorLookUpTable*) as_CSB(srcSpace)->colorLUT()))
+    , fSrcToDst(srcToDst)
 {
     // Build tables to transform src gamma to linear.
     switch (srcSpace->gammaNamed()) {
         case SkColorSpace::kSRGB_GammaNamed:
             fSrcGammaTables[0] = fSrcGammaTables[1] = fSrcGammaTables[2] = sk_linear_from_srgb;
             break;
         case SkColorSpace::k2Dot2Curve_GammaNamed:
             fSrcGammaTables[0] = fSrcGammaTables[1] = fSrcGammaTables[2] = sk_linear_from_2dot2;
             break;
         case SkColorSpace::kLinear_GammaNamed:
@@ skipping 126 matching lines @@
                     build_table_linear_to_gamma(&fDstGammaTableStorage[i * kDstGammaTableSize],
                                                 kDstGammaTableSize, curve.fG, curve.fA, curve.fB,
                                                 curve.fC, curve.fD, curve.fE, curve.fF);
                     fDstGammaTables[i] = &fDstGammaTableStorage[i * kDstGammaTableSize];
                 }
             }
         }
     }
 }
 
+static float byte_to_float(uint8_t byte) {
+    return ((float) byte) * (1.0f / 255.0f);
+}
+
 // Clamp to the 0-1 range.
 static float clamp_normalized_float(float v) {
     if (v > 1.0f) {
         return 1.0f;
     } else if ((v < 0.0f) || (v != v)) {
         return 0.0f;
     } else {
         return v;
     }
 }
 
+static void interp_3d_clut(float dst[3], float src[3], const SkColorLookUpTable* colorLUT) {
+    // 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.
+    float* ptr = &colorLUT->fTable[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++;
+    }
+}
+
 void SkDefaultXform::xform_RGB1_8888(uint32_t* dst, const uint32_t* src, uint32_t len) const {
     while (len-- > 0) {
+        uint8_t r = (*src >>  0) & 0xFF,
+                g = (*src >>  8) & 0xFF,
+                b = (*src >> 16) & 0xFF;
+
+        if (fColorLUT) {
+            float in[3];
+            float out[3];
+
+            in[0] = byte_to_float(r);
+            in[1] = byte_to_float(g);
+            in[2] = byte_to_float(b);
+
+            interp_3d_clut(out, in, fColorLUT.get());
+
+            r = sk_float_round2int(255.0f * clamp_normalized_float(out[0]));
+            g = sk_float_round2int(255.0f * clamp_normalized_float(out[1]));
+            b = sk_float_round2int(255.0f * clamp_normalized_float(out[2]));
+        }
+
         // Convert to linear.
         float srcFloats[3];
-        srcFloats[0] = fSrcGammaTables[0][(*src >>  0) & 0xFF];
-        srcFloats[1] = fSrcGammaTables[1][(*src >>  8) & 0xFF];
-        srcFloats[2] = fSrcGammaTables[2][(*src >> 16) & 0xFF];
+        srcFloats[0] = fSrcGammaTables[0][r];
+        srcFloats[1] = fSrcGammaTables[1][g];
+        srcFloats[2] = fSrcGammaTables[2][b];
 
         // Convert to dst gamut.
         float dstFloats[3];
         dstFloats[0] = srcFloats[0] * fSrcToDst.getFloat(0, 0) +
                        srcFloats[1] * fSrcToDst.getFloat(1, 0) +
                        srcFloats[2] * fSrcToDst.getFloat(2, 0) + fSrcToDst.getFloat(3, 0);
         dstFloats[1] = srcFloats[0] * fSrcToDst.getFloat(0, 1) +
                        srcFloats[1] * fSrcToDst.getFloat(1, 1) +
                        srcFloats[2] * fSrcToDst.getFloat(2, 1) + fSrcToDst.getFloat(3, 1);
         dstFloats[2] = srcFloats[0] * fSrcToDst.getFloat(0, 2) +
                        srcFloats[1] * fSrcToDst.getFloat(1, 2) +
                        srcFloats[2] * fSrcToDst.getFloat(2, 2) + fSrcToDst.getFloat(3, 2);
 
         // Clamp to 0-1.
         dstFloats[0] = clamp_normalized_float(dstFloats[0]);
         dstFloats[1] = clamp_normalized_float(dstFloats[1]);
         dstFloats[2] = clamp_normalized_float(dstFloats[2]);
 
         // Convert to dst gamma.
-        uint8_t r = fDstGammaTables[0][sk_float_round2int((kDstGammaTableSize - 1) * dstFloats[0])];
-        uint8_t g = fDstGammaTables[1][sk_float_round2int((kDstGammaTableSize - 1) * dstFloats[1])];
-        uint8_t b = fDstGammaTables[2][sk_float_round2int((kDstGammaTableSize - 1) * dstFloats[2])];
+        r = fDstGammaTables[0][sk_float_round2int((kDstGammaTableSize - 1) * dstFloats[0])];
+        g = fDstGammaTables[1][sk_float_round2int((kDstGammaTableSize - 1) * dstFloats[1])];
+        b = fDstGammaTables[2][sk_float_round2int((kDstGammaTableSize - 1) * dstFloats[2])];
 
         *dst = SkPackARGB32NoCheck(0xFF, r, g, b);
 
         dst++;
         src++;
     }
 }