From 53ea91139ad9f5ec32fba2cfd167277a5cc3f443 Mon Sep 17 00:00:00 2001 From: "skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com" Date: Thu, 17 May 2018 19:25:40 +0000 Subject: [PATCH] Roll skia/third_party/skcms 3e527c6..5bfec77 (1 commits) https://skia.googlesource.com/skcms.git/+log/3e527c6..5bfec77 2018-05-17 mtklein@chromium.org rm GaussNewton.[ch] 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. CQ_INCLUDE_TRYBOTS=master.tryserver.blink:linux_trusty_blink_rel TBR=herb@google.com Change-Id: Idf7e24d60750f69db9d09a71e9665073380b8912 Reviewed-on: https://skia-review.googlesource.com/128987 Reviewed-by: skcms-skia-autoroll Commit-Queue: skcms-skia-autoroll --- third_party/skcms/skcms.c | 1 - third_party/skcms/skcms.gni | 2 - third_party/skcms/src/GaussNewton.c | 94 -------------------- third_party/skcms/src/GaussNewton.h | 25 ------ third_party/skcms/src/TransferFunction.c | 105 +++++++++++++++++++---- third_party/skcms/version.sha1 | 2 +- 6 files changed, 91 insertions(+), 138 deletions(-) delete mode 100644 third_party/skcms/src/GaussNewton.c delete mode 100644 third_party/skcms/src/GaussNewton.h diff --git a/third_party/skcms/skcms.c b/third_party/skcms/skcms.c index 80ea143cb0..6b1abf59dd 100644 --- a/third_party/skcms/skcms.c +++ b/third_party/skcms/skcms.c @@ -8,7 +8,6 @@ // skcms.c is a unity build target for skcms, #including every other C source file. #include "src/Curve.c" -#include "src/GaussNewton.c" #include "src/ICCProfile.c" #include "src/LinearAlgebra.c" #include "src/PortableMath.c" diff --git a/third_party/skcms/skcms.gni b/third_party/skcms/skcms.gni index 64eeb16111..8577de3465 100644 --- a/third_party/skcms/skcms.gni +++ b/third_party/skcms/skcms.gni @@ -6,8 +6,6 @@ skcms_sources = [ "src/Curve.c", "src/Curve.h", - "src/GaussNewton.c", - "src/GaussNewton.h", "src/ICCProfile.c", "src/LinearAlgebra.c", "src/LinearAlgebra.h", diff --git a/third_party/skcms/src/GaussNewton.c b/third_party/skcms/src/GaussNewton.c deleted file mode 100644 index 61aeddbcdb..0000000000 --- a/third_party/skcms/src/GaussNewton.c +++ /dev/null @@ -1,94 +0,0 @@ -/* - * Copyright 2018 Google Inc. - * - * Use of this source code is governed by a BSD-style license that can be - * found in the LICENSE file. - */ - -#include "../skcms.h" -#include "GaussNewton.h" -#include "LinearAlgebra.h" -#include "PortableMath.h" -#include "TransferFunction.h" -#include -#include - -bool skcms_gauss_newton_step(float (*rg)(float x, const void*, const float P[3], float dfdP[3]), - const void* ctx, - float P[3], - float x0, float dx, int N) { - // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing. - // - // We want to do P' = P + (Jf^T Jf)^-1 Jf^T r(P), - // where r(P) is the residual vector - // and Jf is the Jacobian matrix of f(), ∂r/∂P. - // - // 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 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 [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 [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 - // 2) evaluate Jf^T r(P), call that rhs - // 3) invert lhs - // 4) multiply inverse lhs by rhs - // - // This is a friendly implementation strategy because we don't have to have any - // buffers that scale with N, and equally nice don't have to perform any matrix - // operations that are variable size. - // - // 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 [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_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, - // so we'll have to update lhs and rhs at the same time. - for (int i = 0; i < N; i++) { - float x = x0 + i*dx; - - float dfdP[3] = {0,0,0}; - float resid = rg(x,ctx,P, dfdP); - - 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 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 < 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_Matrix3x3 lhs_inv; - if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) { - return false; - } - - // 4) multiply inverse lhs by 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]; - return isfinitef_(P[0]) && isfinitef_(P[1]) && isfinitef_(P[2]); -} diff --git a/third_party/skcms/src/GaussNewton.h b/third_party/skcms/src/GaussNewton.h deleted file mode 100644 index 5de9927ed3..0000000000 --- a/third_party/skcms/src/GaussNewton.h +++ /dev/null @@ -1,25 +0,0 @@ -/* - * Copyright 2018 Google Inc. - * - * Use of this source code is governed by a BSD-style license that can be - * found in the LICENSE file. - */ - -#pragma once - -#include - -// One Gauss-Newton step, tuning up to 3 parameters P to minimize [ r(x,ctx) ]^2. -// -// rg: residual function r(x,P) to minimize, and gradient at x in dfdP -// ctx: arbitrary context argument passed to rg -// P: in-out, both your initial guess for parameters of r(), and our updated values -// x0,dx,N: N x-values to test with even dx spacing, [x0, x0+dx, x0+2dx, ...] -// -// 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 (*rg)(float x, const void*, const float P[3], float dfdP[3]), - const void* ctx, - float P[3], - float x0, float dx, int N); diff --git a/third_party/skcms/src/TransferFunction.c b/third_party/skcms/src/TransferFunction.c index 8f8fb36c9d..0e3467c13f 100644 --- a/third_party/skcms/src/TransferFunction.c +++ b/third_party/skcms/src/TransferFunction.c @@ -7,7 +7,6 @@ #include "../skcms.h" #include "Curve.h" -#include "GaussNewton.h" #include "LinearAlgebra.h" #include "Macros.h" #include "PortableMath.h" @@ -151,19 +150,15 @@ bool skcms_TransferFunction_invert(const skcms_TransferFunction* src, skcms_Tran // ∂r/∂b = g(ay + b)^(g-1) // - g(ad + b)^(g-1) -typedef struct { - const skcms_Curve* curve; - const skcms_TransferFunction* tf; -} rg_nonlinear_arg; - // Return the residual of roundtripping skcms_Curve(x) through f_inv(y) with parameters P, // and fill out the gradient of the residual into dfdP. -static float rg_nonlinear(float x, const void* ctx, const float P[3], float dfdP[3]) { - const rg_nonlinear_arg* arg = (const rg_nonlinear_arg*)ctx; +static float rg_nonlinear(float x, + const skcms_Curve* curve, + const skcms_TransferFunction* tf, + const float P[3], + float dfdP[3]) { + const float y = skcms_eval_curve(curve, x); - const float y = skcms_eval_curve(arg->curve, x); - - const skcms_TransferFunction* tf = arg->tf; const float g = P[0], a = P[1], b = P[2], c = tf->c, d = tf->d, f = tf->f; @@ -230,6 +225,87 @@ int skcms_fit_linear(const skcms_Curve* curve, int N, float tol, float* c, float return lin_points; } +static bool gauss_newton_step(const skcms_Curve* curve, + const skcms_TransferFunction* tf, + float P[3], + float x0, float dx, int N) { + // We'll sample x from the range [x0,x1] (both inclusive) N times with even spacing. + // + // We want to do P' = P + (Jf^T Jf)^-1 Jf^T r(P), + // where r(P) is the residual vector + // and Jf is the Jacobian matrix of f(), ∂r/∂P. + // + // 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 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 [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 [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 + // 2) evaluate Jf^T r(P), call that rhs + // 3) invert lhs + // 4) multiply inverse lhs by rhs + // + // This is a friendly implementation strategy because we don't have to have any + // buffers that scale with N, and equally nice don't have to perform any matrix + // operations that are variable size. + // + // 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 [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 gauss_newton_step() function ever failing. + + // 0) start off with lhs and rhs safely zeroed. + 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, + // so we'll have to update lhs and rhs at the same time. + for (int i = 0; i < N; i++) { + float x = x0 + i*dx; + + float dfdP[3] = {0,0,0}; + float resid = rg_nonlinear(x,curve,tf,P, dfdP); + + 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 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 < 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_Matrix3x3 lhs_inv; + if (!skcms_Matrix3x3_invert(&lhs, &lhs_inv)) { + return false; + } + + // 4) multiply inverse lhs by 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]; + return isfinitef_(P[0]) && isfinitef_(P[1]) && isfinitef_(P[2]); +} + + // Fit the points in [L,N) to the non-linear piece of tf, or return false if we can't. static bool fit_nonlinear(const skcms_Curve* curve, int L, int N, skcms_TransferFunction* tf) { float P[3] = { tf->g, tf->a, tf->b }; @@ -250,10 +326,9 @@ static bool fit_nonlinear(const skcms_Curve* curve, int L, int N, skcms_Transfer assert (P[1] >= 0 && P[1] * tf->d + P[2] >= 0); - rg_nonlinear_arg arg = { curve, tf}; - if (!skcms_gauss_newton_step(rg_nonlinear, &arg, - P, - L*dx, dx, N-L)) { + if (!gauss_newton_step(curve, tf, + P, + L*dx, dx, N-L)) { return false; } } diff --git a/third_party/skcms/version.sha1 b/third_party/skcms/version.sha1 index facaf94c06..cfc5d9610e 100755 --- a/third_party/skcms/version.sha1 +++ b/third_party/skcms/version.sha1 @@ -1 +1 @@ -3e527c628503e17ba8f4756b00b52b1b17f24cfd \ No newline at end of file +5bfec7723e0fefd2a257e8be4710282cd033eeb2 \ No newline at end of file