Roll skia/third_party/skcms c3b186a..5a327ce (1 commits)

https://skia.googlesource.com/skcms.git/+log/c3b186a..5a327ce

2018-04-26 mtklein@chromium.org drop GaussNewton to max 3 parameters


The AutoRoll server is located here: https://skcms-skia-roll.skia.org

Documentation for the AutoRoller is here:
https://skia.googlesource.com/buildbot/+/master/autoroll/README.md

If the roll is causing failures, please contact the current sheriff, who should
be CC'd on the roll, and stop the roller if necessary.



TBR=stani@google.com

Change-Id: Ic242cd2f2d2e3ba599362b6f19725c831b1e9f61
Reviewed-on: https://skia-review.googlesource.com/124143
Reviewed-by: skcms-skia-autoroll <skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com>

7 files changed, 39 insertions, 133 deletions

diff --git a/third_party/skcms/src/GaussNewton.c b/third_party/skcms/src/GaussNewton.c
index 842f7ed804..5bbc64475d 100644
--- a/third_party/skcms/src/GaussNewton.c
+++ b/third_party/skcms/src/GaussNewton.c
@@ -14,11 +14,11 @@
 
 bool skcms_gauss_newton_step(float (*     t)(float x, const void*),
                              const void* t_ctx,
-                             float (*     f)(float x, const void*, const float P[4]),
+                             float (*     f)(float x, const void*, const float P[3]),
                              const void* f_ctx,
-                             void  (*grad_f)(float x, const void*, const float P[4], float dfdP[4]),
+                             void  (*grad_f)(float x, const void*, const float P[3], float dfdP[3]),
                              const void* g_ctx,
-                             float P[4],
+                             float P[3],
                              float x0, float x1, int N) {
     // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing.
     //
@@ -28,12 +28,12 @@ bool skcms_gauss_newton_step(float (*     t)(float x, const void*),
     //
     // Let's review the shape of each of these expressions:
     //   r(P)   is [N x 1], a column vector with one entry per value of x tested
-    //   Jf     is [N x 4], a matrix with an entry for each (x,P) pair
-    //   Jf^T   is [4 x N], the transpose of Jf
+    //   Jf     is [N x 3], a matrix with an entry for each (x,P) pair
+    //   Jf^T   is [3 x N], the transpose of Jf
     //
-    //   Jf^T Jf   is [4 x N] * [N x 4] == [4 x 4], a 4x4 matrix,
+    //   Jf^T Jf   is [3 x N] * [N x 3] == [3 x 3], a 3x3 matrix,
     //                                              and so is its inverse (Jf^T Jf)^-1
-    //   Jf^T r(P) is [4 x N] * [N x 1] == [4 x 1], a column vector with the same shape as P
+    //   Jf^T r(P) is [3 x N] * [N x 1] == [3 x 1], a column vector with the same shape as P
     //
     // Our implementation strategy to get to the final ∆P is
     //   1) evaluate Jf^T Jf,   call that lhs
@@ -48,13 +48,13 @@ bool skcms_gauss_newton_step(float (*     t)(float x, const void*),
     // Other implementation strategies could trade this off, e.g. evaluating the
     // pseudoinverse of Jf ( (Jf^T Jf)^-1 Jf^T ) directly, then multiplying that by
     // the residuals.  That would probably require implementing singular value
-    // decomposition, and would create a [4 x N] matrix to be multiplied by the
+    // decomposition, and would create a [3 x N] matrix to be multiplied by the
     // [N x 1] residual vector, but on the upside I think that'd eliminate the
     // possibility of this skcms_gauss_newton_step() function ever failing.
 
     // 0) start off with lhs and rhs safely zeroed.
-    skcms_Matrix4x4 lhs = {{ {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} }};
-    skcms_Vector4   rhs = {  {0,0,0,0} };
+    skcms_Matrix3x3 lhs = {{ {0,0,0}, {0,0,0}, {0,0,0} }};
+    skcms_Vector3   rhs = {  {0,0,0} };
 
     // 1,2) evaluate lhs and evaluate rhs
     //   We want to evaluate Jf only once, but both lhs and rhs involve Jf^T,
@@ -65,38 +65,37 @@ bool skcms_gauss_newton_step(float (*     t)(float x, const void*),
 
         float resid = t(x,t_ctx) - f(x,f_ctx,P);
 
-        float dfdP[4] = {0,0,0,0};
+        float dfdP[3] = {0,0,0};
         grad_f(x,g_ctx,P, dfdP);
 
-        for (int r = 0; r < 4; r++) {
-            for (int c = 0; c < 4; c++) {
+        for (int r = 0; r < 3; r++) {
+            for (int c = 0; c < 3; c++) {
                 lhs.vals[r][c] += dfdP[r] * dfdP[c];
             }
             rhs.vals[r] += dfdP[r] * resid;
         }
     }
 
-    // If any of the 4 P parameters are unused, this matrix will be singular.
+    // If any of the 3 P parameters are unused, this matrix will be singular.
     // Detect those cases and fix them up to indentity instead, so we can invert.
-    for (int k = 0; k < 4; k++) {
-        if (lhs.vals[0][k]==0 && lhs.vals[1][k]==0 && lhs.vals[2][k]==0 && lhs.vals[3][k]==0 &&
-            lhs.vals[k][0]==0 && lhs.vals[k][1]==0 && lhs.vals[k][2]==0 && lhs.vals[k][3]==0) {
+    for (int k = 0; k < 3; k++) {
+        if (lhs.vals[0][k]==0 && lhs.vals[1][k]==0 && lhs.vals[2][k]==0 &&
+            lhs.vals[k][0]==0 && lhs.vals[k][1]==0 && lhs.vals[k][2]==0) {
             lhs.vals[k][k] = 1;
         }
     }
 
     // 3) invert lhs
-    skcms_Matrix4x4 lhs_inv;
-    if (!skcms_Matrix4x4_invert(&lhs, &lhs_inv)) {
+    skcms_Matrix3x3 lhs_inv;
+    if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) {
         return false;
     }
 
     // 4) multiply inverse lhs by rhs
-    skcms_Vector4 dP = skcms_Matrix4x4_Vector4_mul(&lhs_inv, &rhs);
+    skcms_Vector3 dP = skcms_MV_mul(&lhs_inv, &rhs);
     P[0] += dP.vals[0];
     P[1] += dP.vals[1];
     P[2] += dP.vals[2];
-    P[3] += dP.vals[3];
     return true;
 }
 
diff --git a/third_party/skcms/src/GaussNewton.h b/third_party/skcms/src/GaussNewton.h
index 16ad80bcf7..2cabbd0075 100644
--- a/third_party/skcms/src/GaussNewton.h
+++ b/third_party/skcms/src/GaussNewton.h
@@ -9,7 +9,7 @@
 
 #include <stdbool.h>
 
-// One Gauss-Newton step, tuning up to 4 parameters P to minimize [ t(x,ctx) - f(x,P) ]^2.
+// One Gauss-Newton step, tuning up to 3 parameters P to minimize [ t(x,ctx) - f(x,P) ]^2.
 //
 //   t:        target function of x to approximate
 //   t_ctx:    any context needed for t, passed blindly into calls to t()
@@ -18,16 +18,16 @@
 //   P:        in-out, both your initial guess for parameters of f(), and our updated values
 //   x0,x1,N:  N x-values to test in [x0,x1] (both inclusive) with even spacing
 //
-// If you have fewer than 4 parameters, set the unused P to zero, don't touch their dfdP.
+// If you have fewer than 3 parameters, set the unused P to zero, don't touch their dfdP.
 //
 // Returns true and updates P on success, or returns false on failure.
 bool skcms_gauss_newton_step(float (*     t)(float x, const void*),
                              const void* t_ctx,
-                             float (*     f)(float x, const void*, const float P[4]),
+                             float (*     f)(float x, const void*, const float P[3]),
                              const void* f_ctx,
-                             void  (*grad_f)(float x, const void*, const float P[4], float dfdP[4]),
+                             void  (*grad_f)(float x, const void*, const float P[3], float dfdP[3]),
                              const void* g_ctx,
-                             float P[4],
+                             float P[3],
                              float x0, float x1, int N);
 
 // A target function for skcms_Curve, passed as ctx.
diff --git a/third_party/skcms/src/LinearAlgebra.c b/third_party/skcms/src/LinearAlgebra.c
index 9ef7079fdc..51127e6d1d 100644
--- a/third_party/skcms/src/LinearAlgebra.c
+++ b/third_party/skcms/src/LinearAlgebra.c
@@ -10,92 +10,6 @@
 #include "PortableMath.h"
 #include <float.h>
 
-bool skcms_Matrix4x4_invert(const skcms_Matrix4x4* src, skcms_Matrix4x4* dst) {
-    double a00 = src->vals[0][0],
-           a01 = src->vals[1][0],
-           a02 = src->vals[2][0],
-           a03 = src->vals[3][0],
-           a10 = src->vals[0][1],
-           a11 = src->vals[1][1],
-           a12 = src->vals[2][1],
-           a13 = src->vals[3][1],
-           a20 = src->vals[0][2],
-           a21 = src->vals[1][2],
-           a22 = src->vals[2][2],
-           a23 = src->vals[3][2],
-           a30 = src->vals[0][3],
-           a31 = src->vals[1][3],
-           a32 = src->vals[2][3],
-           a33 = src->vals[3][3];
-
-    double b00 = a00*a11 - a01*a10,
-           b01 = a00*a12 - a02*a10,
-           b02 = a00*a13 - a03*a10,
-           b03 = a01*a12 - a02*a11,
-           b04 = a01*a13 - a03*a11,
-           b05 = a02*a13 - a03*a12,
-           b06 = a20*a31 - a21*a30,
-           b07 = a20*a32 - a22*a30,
-           b08 = a20*a33 - a23*a30,
-           b09 = a21*a32 - a22*a31,
-           b10 = a21*a33 - a23*a31,
-           b11 = a22*a33 - a23*a32;
-
-    double determinant = b00*b11
-                       - b01*b10
-                       + b02*b09
-                       + b03*b08
-                       - b04*b07
-                       + b05*b06;
-
-    if (determinant == 0) {
-        return false;
-    }
-
-    double invdet = 1.0 / determinant;
-    if (invdet > +FLT_MAX || invdet < -FLT_MAX || !isfinitef_((float)invdet)) {
-        return false;
-    }
-
-    b00 *= invdet;
-    b01 *= invdet;
-    b02 *= invdet;
-    b03 *= invdet;
-    b04 *= invdet;
-    b05 *= invdet;
-    b06 *= invdet;
-    b07 *= invdet;
-    b08 *= invdet;
-    b09 *= invdet;
-    b10 *= invdet;
-    b11 *= invdet;
-
-    dst->vals[0][0] = (float)( a11*b11 - a12*b10 + a13*b09 );
-    dst->vals[1][0] = (float)( a02*b10 - a01*b11 - a03*b09 );
-    dst->vals[2][0] = (float)( a31*b05 - a32*b04 + a33*b03 );
-    dst->vals[3][0] = (float)( a22*b04 - a21*b05 - a23*b03 );
-    dst->vals[0][1] = (float)( a12*b08 - a10*b11 - a13*b07 );
-    dst->vals[1][1] = (float)( a00*b11 - a02*b08 + a03*b07 );
-    dst->vals[2][1] = (float)( a32*b02 - a30*b05 - a33*b01 );
-    dst->vals[3][1] = (float)( a20*b05 - a22*b02 + a23*b01 );
-    dst->vals[0][2] = (float)( a10*b10 - a11*b08 + a13*b06 );
-    dst->vals[1][2] = (float)( a01*b08 - a00*b10 - a03*b06 );
-    dst->vals[2][2] = (float)( a30*b04 - a31*b02 + a33*b00 );
-    dst->vals[3][2] = (float)( a21*b02 - a20*b04 - a23*b00 );
-    dst->vals[0][3] = (float)( a11*b07 - a10*b09 - a12*b06 );
-    dst->vals[1][3] = (float)( a00*b09 - a01*b07 + a02*b06 );
-    dst->vals[2][3] = (float)( a31*b01 - a30*b03 - a32*b00 );
-    dst->vals[3][3] = (float)( a20*b03 - a21*b01 + a22*b00 );
-
-    for (int r = 0; r < 4; ++r)
-    for (int c = 0; c < 4; ++c) {
-        if (!isfinitef_(dst->vals[r][c])) {
-            return false;
-        }
-    }
-    return true;
-}
-
 bool skcms_Matrix3x3_invert(const skcms_Matrix3x3* src, skcms_Matrix3x3* dst) {
     double a00 = src->vals[0][0],
            a01 = src->vals[1][0],
@@ -123,7 +37,7 @@ bool skcms_Matrix3x3_invert(const skcms_Matrix3x3* src, skcms_Matrix3x3* dst) {
     }
 
     double invdet = 1.0 / determinant;
-    if (!isfinitef_((float)invdet)) {
+    if (invdet > +FLT_MAX || invdet < -FLT_MAX || !isfinitef_((float)invdet)) {
         return false;
     }
 
@@ -153,13 +67,12 @@ bool skcms_Matrix3x3_invert(const skcms_Matrix3x3* src, skcms_Matrix3x3* dst) {
     return true;
 }
 
-skcms_Vector4 skcms_Matrix4x4_Vector4_mul(const skcms_Matrix4x4* m, const skcms_Vector4* v) {
-    skcms_Vector4 dst = {{0,0,0,0}};
-    for (int row = 0; row < 4; ++row) {
+skcms_Vector3 skcms_MV_mul(const skcms_Matrix3x3* m, const skcms_Vector3* v) {
+    skcms_Vector3 dst = {{0,0,0}};
+    for (int row = 0; row < 3; ++row) {
         dst.vals[row] = m->vals[row][0] * v->vals[0]
                       + m->vals[row][1] * v->vals[1]
-                      + m->vals[row][2] * v->vals[2]
-                      + m->vals[row][3] * v->vals[3];
+                      + m->vals[row][2] * v->vals[2];
     }
     return dst;
 }
diff --git a/third_party/skcms/src/LinearAlgebra.h b/third_party/skcms/src/LinearAlgebra.h
index 95a2cf9bf9..47f53d9230 100644
--- a/third_party/skcms/src/LinearAlgebra.h
+++ b/third_party/skcms/src/LinearAlgebra.h
@@ -9,15 +9,9 @@
 
 #include <stdbool.h>
 
-// A row-major 4x4 matrix (ie vals[row][col])
-typedef struct { float vals[4][4]; } skcms_Matrix4x4;
-// (skcms_Matrix3x3 is in skcms.h, exactly analogous to skcms_Matrix4x4.)
-
-typedef struct { float vals[4]; } skcms_Vector4;
-
+typedef struct { float vals[3]; } skcms_Vector3;
 
 // It is _not_ safe to alias the pointers to invert in-place.
-bool skcms_Matrix4x4_invert(const skcms_Matrix4x4*, skcms_Matrix4x4*);
 bool skcms_Matrix3x3_invert(const skcms_Matrix3x3*, skcms_Matrix3x3*);
 
-skcms_Vector4 skcms_Matrix4x4_Vector4_mul(const skcms_Matrix4x4*, const skcms_Vector4*);
+skcms_Vector3 skcms_MV_mul(const skcms_Matrix3x3*, const skcms_Vector3*);
diff --git a/third_party/skcms/src/PolyTF.c b/third_party/skcms/src/PolyTF.c
index f7f9a3fb3d..0137a3ca8b 100644
--- a/third_party/skcms/src/PolyTF.c
+++ b/third_party/skcms/src/PolyTF.c
@@ -41,7 +41,7 @@
 // It's important to evaluate as f(x) as A(x^3-1) + B(x^2-1) + 1
 // and not Ax^3 + Bx^2 + (1-A-B) to ensure that f(1.0f) == 1.0f.
 
-static float eval_poly_tf(float x, const void* ctx, const float P[4]) {
+static float eval_poly_tf(float x, const void* ctx, const float P[3]) {
     const skcms_PolyTF* tf = (const skcms_PolyTF*)ctx;
 
     float A = P[0],
@@ -53,7 +53,7 @@ static float eval_poly_tf(float x, const void* ctx, const float P[4]) {
                  : A*(x*x*x-1) + B*(x*x-1) + 1;
 }
 
-static void grad_poly_tf(float x, const void* ctx, const float P[4], float dfdP[4]) {
+static void grad_poly_tf(float x, const void* ctx, const float P[3], float dfdP[3]) {
     const skcms_PolyTF* tf = (const skcms_PolyTF*)ctx;
     (void)P;
     float D = tf->D;
@@ -114,7 +114,7 @@ static bool fit_poly_tf(const skcms_Curve* curve, skcms_PolyTF* tf) {
     }
 
     // Start with guess A = 0, i.e. f(x) ≈ x^2.
-    float P[4] = {0, 0,0,0};
+    float P[3] = {0, 0,0};
     for (int i = 0; i < 3; i++) {
         if (!skcms_gauss_newton_step(skcms_eval_curve, curve,
                                      eval_poly_tf, tf,
diff --git a/third_party/skcms/src/TransferFunction.c b/third_party/skcms/src/TransferFunction.c
index c32d1204a3..6daf0cf4d3 100644
--- a/third_party/skcms/src/TransferFunction.c
+++ b/third_party/skcms/src/TransferFunction.c
@@ -143,7 +143,7 @@ bool skcms_TransferFunction_invert(const skcms_TransferFunction* src, skcms_Tran
 //    ∂tf/∂b =  g(ax + b)^(g-1)
 //           -  g(ad + b)^(g-1)
 
-static float eval_nonlinear(float x, const void* ctx, const float P[4]) {
+static float eval_nonlinear(float x, const void* ctx, const float P[3]) {
     const skcms_TransferFunction* tf = (const skcms_TransferFunction*)ctx;
 
     const float g = P[0],  a = P[1],  b = P[2],
@@ -158,7 +158,7 @@ static float eval_nonlinear(float x, const void* ctx, const float P[4]) {
          - powf_(D, g)
          + c*d + f;
 }
-static void grad_nonlinear(float x, const void* ctx, const float P[4], float dfdP[4]) {
+static void grad_nonlinear(float x, const void* ctx, const float P[3], float dfdP[3]) {
     const skcms_TransferFunction* tf = (const skcms_TransferFunction*)ctx;
 
     const float g = P[0], a = P[1], b = P[2],
@@ -222,7 +222,7 @@ int skcms_fit_linear(const skcms_Curve* curve, int N, float tol, float* c, float
 
 // Fit the points in [start,N) to the non-linear piece of tf, or return false if we can't.
 static bool fit_nonlinear(const skcms_Curve* curve, int start, int N, skcms_TransferFunction* tf) {
-    float P[4] = { tf->g, tf->a, tf->b, 0 };
+    float P[3] = { tf->g, tf->a, tf->b };
 
     // No matter where we start, x_scale should always represent N even steps from 0 to 1.
     const float x_scale = 1.0f / (N-1);
diff --git a/third_party/skcms/version.sha1 b/third_party/skcms/version.sha1
index 6cdb91df73..7c2fcaddc8 100755
--- a/third_party/skcms/version.sha1
+++ b/third_party/skcms/version.sha1
@@ -1 +1 @@
-c3b186a912acceb3dec3f49b01cfa02977cbde88
\ No newline at end of file
+5a327ce9eb094efdb75893d1ac6d46637e006f14
\ No newline at end of file