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>
This commit is contained in:
skcms-skia-autoroll@skia-buildbots.google.com.iam.gserviceaccount.com 2018-05-17 19:25:40 +00:00 committed by Skia Commit-Bot
parent 8f288d9399
commit 53ea91139a
6 changed files with 91 additions and 138 deletions

View File

@ -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"

View File

@ -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",

View File

@ -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 <assert.h>
#include <string.h>
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]);
}

View File

@ -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 <stdbool.h>
// 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);

View File

@ -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;
}
}

View File

@ -1 +1 @@
3e527c628503e17ba8f4756b00b52b1b17f24cfd
5bfec7723e0fefd2a257e8be4710282cd033eeb2