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:

• ChamberlinAfrica
• ModifiedStereographicAlaska
• ModifiedStereographicGs48
• ModifiedStereographicMiller
• ModifiedStereographicLee
• ModifiedStereographicGs50 (the part over Europe is broken)
• TwoPointAzimuthalUsa

It doesn't work with

• InterruptedSinuMollweide (problem of definition of the invert in interrupted.js)

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 <a href='http://bl.ocks.org/Fil/'>Fil</a>'s block: <a href='http://bl.ocks.org/Fil/cb3718040261307ef04b4031615e0c5c'>Interrupted Homolosine</a>