So I’ve started putting some code together. For all of the base polynomials used by lensfun the coefficients can be determined exactly through the solution of a suitably framed linear least squares problem (this is important, as non-linear least squares is a real pain and seldom suitable for real-time applications).
I’ve done a first pass of almost everything except for the distortion molestation routine. This is responsible from converting the coefficients from their natural form to that used by the lensfun database (the irony being that not too long after reading these coefficients in lensfun undoes these transforms to return the polynomial back to its natural form). Finishing this routine will require me to write a full multi-variate implementation of Newton’s method. That and a huge amount of testing and debugging.
#include <stdio.h>
#include <tgmath.h>
static inline
void givens(double f, double g, double *c, double *s)
{
if (g == 0.0)
{
*c = 1.0;
*s = 0.0;
}
else if (f == 0.0)
{
*c = 0.0;
*s = copysign(1.0, g);
}
else
{
double l = sqrt(f*f + g*g);
*c = fabs(f) / l;
*s = copysign(1.0, f)*g / l;
}
}
static inline
void givens_qr(double *A, double *b, int m, int n)
{
for (int j = 0; j < n; j++)
{
for (int i = m - 1; i > j; i--)
{
int k = i - 1;
double c, s;
givens(A[k*n + j], A[i*n + j], &c, &s);
// Apply to A
for (int l = j; l < n; l++)
{
double p = A[k*n + l], q = A[i*n + l];
A[k*n + l] = c*p + s*q;
A[i*n + l] = -s*p + c*q;
}
// Apply to b
double p = b[k], q = b[i];
b[k] = c*p + s*q;
b[i] = -s*p + c*q;
}
}
}
static inline
void solve_bwd(double *R, double *b, double *x, int n)
{
for (int i = n - 1; i >= 0; i--)
{
double t = b[i];
for (int j = i + 1; j < n; j++)
t -= R[i*n + j]*x[j];
x[i] = t / R[i*n + i];
}
}
static inline
void interpolate(double *spcoeffs, int nc, double *xv, double *yv, int m)
{
for (int i = 0; i < m; i++)
{
xv[i] = (double) i / (m - 1);
for (int j = 1; j < nc; j++)
{
double l = (j - 1.0) / (nc - 1), r = (double) j / (nc - 1);
if (xv[i] >= l && xv[i] <= r)
{
double dydx = (nc - 1)*(spcoeffs[j] - spcoeffs[j - 1]);
yv[i] = spcoeffs[j - 1] + (xv[i] - l)*dydx;
break;
}
}
}
}
int
train_sony_dist(const short *coeffs, int nc, double tcoeffs[3])
{
enum { M = 128, N = 3 };
// Check the coefficient count
if (nc != 16)
return -1;
// Pre-process the coefficients
double spcoeffs[16];
for (int i = 0; i < nc; i++)
spcoeffs[i] = coeffs[i]*pow(2.0, -14.0) + 1.0;
// Interpolate the resulting spline onto a fine grid
double xv[M], yv[M];
interpolate(spcoeffs, 16, xv, yv, M);
// Form the A matrix and b vector
double A[M*N], b[M];
for (int i = 0; i < M; i++)
{
A[i*N + 0] = pow(xv[i], 1);
A[i*N + 1] = pow(xv[i], 2);
A[i*N + 2] = pow(xv[i], 3);
b[i] = yv[i] - 1.0;
}
// Solve the linear least squares problem
givens_qr(A, b, M, N);
solve_bwd(A, b, tcoeffs, N);
return 0;
}
void
molest_sony_dist(double fl, const double tcoeffs[3], double mcoeffs[3])
{
//double hugin_scale = fl / (hypot(36.0, 24.0) / hypot(1.5, 1.0) / 2.0);
// TODO: This needs Newton's method
}
int
train_sony_tca(const short *coeffs, int nc, double tcoeffs[3])
{
enum { M = 64, N = 3 };
// Check the coefficient count
if (nc != 16)
return -1;
// Pre-process the coefficients
double spcoeffs[16];
for (int i = 0; i < nc; i++)
spcoeffs[i] = coeffs[i]*pow(2.0, -21.0) + 1.0;
// Interpolate the resulting spline onto a fine grid
double xv[M], yv[M];
interpolate(spcoeffs, 16, xv, yv, M);
// Form the A matrix and b vector
double A[M*N], b[M];
for (int i = 0; i < M; i++)
{
A[i*N + 0] = 1.0;
A[i*N + 1] = pow(xv[i], 1);
A[i*N + 2] = pow(xv[i], 2);
b[i] = yv[i];
}
// Solve the linear least squares problem
givens_qr(A, b, M, N);
solve_bwd(A, b, tcoeffs, N);
return 0;
}
void
molest_sony_tca(double fl, const double tcoeffs[3], double mcoeffs[3])
{
double hugin_scale = fl / (hypot(36.0, 24.0) / hypot(1.5, 1.0) / 2.0);
mcoeffs[0] = tcoeffs[0];
mcoeffs[1] = tcoeffs[1] / hugin_scale;
mcoeffs[2] = tcoeffs[2] / (hugin_scale*hugin_scale);
}
int
train_sony_vig(const short *coeffs, int nc, double tcoeffs[3])
{
enum { M = 128, N = 3 };
// Check the coefficient count
if (nc != 16)
return -1;
// Pre-process the coefficients
double spcoeffs[16];
for (int i = 0; i < nc; i++)
spcoeffs[i] = pow(2.0, pow(2.0, pow(2.0, -13.0)*coeffs[i]) - 1.0);
// Interpolate the resulting spline onto a fine grid
double xv[M], yv[M];
interpolate(spcoeffs, 16, xv, yv, M);
// Form the A matrix and b vector
double A[M*N], b[M];
for (int i = 0; i < M; i++)
{
/*
* Place a strong emphasis on getting the final point correct.
*/
double wgt = (i == M - 1) ? 25.0 : 1.0;
A[i*N + 0] = wgt*pow(xv[i], 2.0);
A[i*N + 1] = wgt*pow(xv[i], 4.0);
A[i*N + 2] = wgt*pow(xv[i], 6.0);
b[i] = wgt*(1.0 / yv[i] - 1.0);
}
// Solve the linear least squares problem
givens_qr(A, b, M, N);
solve_bwd(A, b, tcoeffs, N);
return 0;
}
int
main()
{
// Test TCA
double vcb[3], mvcb[3];
short coeffs[] = { 384, 384, 384, 512, 512, 512, 512, 512, 512, 512, 640, 640, 640, 768, 768, 1024 };
train_sony_tca(coeffs, 16, vcb);
molest_sony_tca(24.0, vcb, mvcb);
printf("TCA: v = %f\tc = %e\tb = %e\n", mvcb[0], mvcb[1], mvcb[2]);
// Test vignetting
double k123[3];
short vcoeffs[] = { 0, 32, 128, 288, 480, 736, 1056, 1376, 1984, 3360, 5024, 6560, 8128, 9728, 11360, 13184 };
train_sony_vig(vcoeffs, 16, k123);
printf("Vig: k1 = %f\tk2 = %f\tk3 = %f\n", k123[0], k123[1], k123[2]);
return 0;
}