Computing the inverse of the ChamberlinAfrica projection for issue #85.
When a projection does not have an invert function, we apply the Newton-Raphson method (see gist 8903368) on a function of the cartesian coordinates (x,y,z). This avoids lots of problems at the poles.
This approach seems to work well with:
It doesn’t work with
The main code is the following:
function solve(project, precision) {
var n = 100,
step = (halfPi - epsilon) / n,
i,
j,
grid = [];
for (i = 0; i <= 4 * n; i++) {
grid[i] = [];
for (j = 0; j <= 2 * n; j++) {
grid[i][j] = project((i - 2 * n) * step, (j - n) * step);
}
}
function invert (x, y) {
// find a start point c "close enough" to x,y
var i, j,
c,
m, min = +Infinity;
// d3.scan
for (i = 0; i <= 4 * n; i++) {
for (j = 0; j <= 2 * n; j++) {
m = abs(x - grid[i][j][0]) + abs(y - grid[i][j][1]);
if (m < min) {
c = [i,j];
min = m;
}
}
}
c = [ step * (c[0] - 2 * n), step * (c[1] - n) ];
c = cartesian(c);
// solve for x,y
var solution = Newton.Solve(g => {
var norm = g[0] * g[0] + g[1] * g[1] + g[2] * g[2];
cartesianScale(g, 1 / sqrt(norm));
var s = spherical(g),
p = project(s[0], s[1]);
return [p[0], p[1], norm];
}, [x, y, 1], { start: c, acc: precision });
var norm = solution[0] * solution[0] + solution[1] * solution[1] + solution[2] * solution[2];
cartesianScale(solution, 1 / sqrt(norm));
return spherical(solution);
}
return invert;
}
/*
* Newton's method for finding roots
*
* code adapted from D.V. Fedorov,
* “Introduction to Numerical Methods with examples in Javascript”
* http://owww.phys.au.dk/~fedorov/nucltheo/Numeric/11/book.pdf
* (licensed under the GPL)
* by Philippe Riviere <philippe.riviere@illisible.net> March 2014
* modified for compatibility with Chrome/Safari
* added a max iterations parameter
*
* Usage: Newton.Solve(Math.exp, 2)); => ~ log(2)
* Newton.Solve(d3.geo.chamberlin(), [200,240])
*/
var Newton = {version: "1.0.0"}; // semver
Newton.Norm = function (v) {
return Math.sqrt(v.reduce(function (s, e) {
return s + e * e
}, 0));
}
Newton.Dot = function (a, b) {
var s = 0;
for (var i in a) s += a[i] * b[i];
return s;
}
// QR - decomposition A=QR of matrix A
Newton.QRDec = function(A){
var m = A.length, R = [], i, j, k;
for (j in A) {
R[j] = [];
for (i in A) R[j][i] = 0;
}
var Q = [];
for (i in A) {
Q[i] = [];
for (j in A[0]) Q[i][j] = A[i][j];
}
// Q is a copy of A
for (i = 0; i < m; i++) {
var e = Q[i],
r = Math.sqrt(Newton.Dot(e, e));
if (r == 0) throw "Newton.QRDec: singular matrix"
R[i][i] = r;
for (k in e) e[k] /= r;
// normalization
for (j = i + 1; j < m; j++) {
var q = Q[j],
s = Newton.Dot(e, q);
for (k in q) q[k] -= s * e[k];
// orthogonalization
R[j][i] = s;
}
}
return [Q, R];
};
// QR - backsubstitution
// input: matrices Q,R, array b; output: array x such that QRx=b
Newton.QRBack = function(Q, R, b) {
var m = Q.length,
c = new Array(m),
x = new Array(m),
i, k, s;
for (i in Q) {
// c = QˆT b
c[i] = 0;
for (k in b) c[i] += Q[i][k] * b[k];
}
for (i = m - 1; i >= 0; i--) {
// back substitution
for (s = 0, k = i + 1; k < m; k++) s += R[k][i] * x[k];
x[i] = (c[i] - s) / R[i][i];
}
return x;
}
// calculates inverse of matrix A
Newton.Inverse = function(A){
var t = Newton.QRDec(A),
Q = t[0],
R = t[1];
var m = [], i, k, n;
for (i in A) {
n = [];
for (k in A) {
n[k] = (k == i ? 1 : 0)
}
m[i] = Newton.QRBack(Q, R, n);
}
return m;
};
Newton.Zero = function (fs, x, opts={} /* acc, dx, max */) {
// Newton's root-finding method
var i, j, k;
if (opts.acc == undefined) opts.acc = 1e-6
if (opts.dx == undefined) opts.dx = 1e-3
if (opts.max == undefined) opts.max = 40 // max iterations
var J = [];
for (i in x) {
J[i] = [];
for (j in x) J[i][j] = 0;
}
var minusfx = [];
var v = fs(x);
if (v == null) throw "unable to compute fs at "+JSON.stringify(x)
for (i in x) minusfx[i] = -v[i];
do {
if (opts.max-- < 0) return null;
for (i in x)
for (k in x) {
// calculate Jacobian
x[k] += opts.dx
v = fs(x);
if (v == null) throw "unable to compute fs at "+JSON.stringify(x)
J[k][i] = (v[i] + minusfx[i]) / opts.dx
x[k] -= opts.dx
}
var t = Newton.QRDec(J),
Q = t[0],
R = t[1],
Dx = Newton.QRBack(Q, R, minusfx)
// Newton's step
var s = 2
do {
// simple backtracking line search
s = s / 2;
var z = [];
for (i in x) {
z[i] = x[i] + s * Dx[i];
}
var minusfz = [];
v = fs(z);
if (v == null) throw "unable to compute fs at "+JSON.stringify(z)
for (i in x) {
minusfz[i] = -v[i];
}
}
while (Newton.Norm(minusfz) > (1 - s / 2) * Newton.Norm(minusfx) && s > 1. / 128)
minusfx = minusfz;
x = z;
// step done
}
while (Newton.Norm(minusfx) > opts.acc)
return x;
}
Newton.Solve = function(fs, res, opts={}){
if (typeof res != 'object') res = [ typeof res == 'number'
? + res
: 0
];
var _fs = fs;
fs = function(x) {
var r = _fs(x);
if (typeof r == 'number') r = [ r ];
for (var i in r) r[i] -= res[i];
return r;
}
var start = [];
if (opts.start) {
start = opts.start;
} else {
for (var i in res) start[i]=0;
}
var n = Newton.Zero(fs, start, opts);
if (n && n.length == 1) n = n[0];
return n;
};forked from Fil‘s block: Interrupted Homolosine