rm GaussNewton.[ch]
Move its impl back into TransferFunction.c, rewrite internal APIs
a little to reflect that it's always calling rg_nonlinear().
Change-Id: I4a2dacbc3efaa0d7456dcdc7467c9095e498c243
Reviewed-on: https://skia-review.googlesource.com/128932
Commit-Queue: Mike Klein <mtklein@chromium.org>
Commit-Queue: Brian Osman <brianosman@google.com>
Auto-Submit: Mike Klein <mtklein@chromium.org>
Reviewed-by: Brian Osman <brianosman@google.com>
diff --git a/build/targets b/build/targets
index c25cfd1..0de391e 100644
--- a/build/targets
+++ b/build/targets
@@ -29,7 +29,6 @@
$out/fuzz/fuzz_main.o $
$out/skcms.o
-build $out/src/GaussNewton.o: compile src/GaussNewton.c
build $out/src/ICCProfile.o: compile src/ICCProfile.c
build $out/src/LinearAlgebra.o: compile src/LinearAlgebra.c
build $out/src/PortableMath.o: compile src/PortableMath.c
diff --git a/skcms.c b/skcms.c
index 80ea143..6b1abf5 100644
--- a/skcms.c
+++ b/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/skcms.gni b/skcms.gni
index 64eeb16..8577de3 100644
--- a/skcms.gni
+++ b/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/src/GaussNewton.c b/src/GaussNewton.c
deleted file mode 100644
index 61aeddb..0000000
--- a/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 <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]);
-}
diff --git a/src/GaussNewton.h b/src/GaussNewton.h
deleted file mode 100644
index 5de9927..0000000
--- a/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 <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);
diff --git a/src/TransferFunction.c b/src/TransferFunction.c
index 8f8fb36..0e3467c 100644
--- a/src/TransferFunction.c
+++ b/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 @@
// ∂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 @@
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 @@
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;
}
}
