diff options
Diffstat (limited to 'third_party')
-rw-r--r-- | third_party/skcms/src/GaussNewton.c | 41 | ||||
-rw-r--r-- | third_party/skcms/src/GaussNewton.h | 10 | ||||
-rw-r--r-- | third_party/skcms/src/LinearAlgebra.c | 97 | ||||
-rw-r--r-- | third_party/skcms/src/LinearAlgebra.h | 10 | ||||
-rw-r--r-- | third_party/skcms/src/PolyTF.c | 6 | ||||
-rw-r--r-- | third_party/skcms/src/TransferFunction.c | 6 | ||||
-rwxr-xr-x | third_party/skcms/version.sha1 | 2 |
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 |