diff options
author | skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com <skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com> | 2018-05-17 19:25:40 +0000 |
---|---|---|
committer | Skia Commit-Bot <skia-commit-bot@chromium.org> | 2018-05-17 19:56:43 +0000 |
commit | 53ea91139ad9f5ec32fba2cfd167277a5cc3f443 (patch) | |
tree | f7db453c32a943da059318b88387eb393684ae44 /third_party/skcms/src/TransferFunction.c | |
parent | 8f288d9399db95cd0a0994f037f6c08410a7c354 (diff) |
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 <skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com>
Commit-Queue: skcms-skia-autoroll <skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com>
Diffstat (limited to 'third_party/skcms/src/TransferFunction.c')
-rw-r--r-- | third_party/skcms/src/TransferFunction.c | 105 |
1 files changed, 90 insertions, 15 deletions
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; - - const float y = skcms_eval_curve(arg->curve, x); +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 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; } } |