Approximate 2-d Newton-Raphson method: inverts Bertin1953, Chamberlin…

… and interruptedSinuMollweide

The general Jacobian inversion method is detailed in Ipbuker & Bildirici 2002

In this implementation, I added a dampening factor when the determinant is too small — which helps with convergence near the singularities of the projection.

Inspired by Robert B. Schmunk’s work on the Bertin1953 inverse for G.Projector, which showed it was possible.

test all inverse projections

closes d3#85 (at long last)

src/bertin.js

@@ -1,14 +1,15 @@
 import {geoProjection as projection} from "d3-geo";
 import {hammerRaw} from "./hammer.js";
 import {cos, pi, sin} from "./math.js";
+import {solve2d} from "./newton.js";
 
 // Bertin 1953 as a modified Briesemeister
 // https://bl.ocks.org/Fil/5b9ee9636dfb6ffa53443c9006beb642
 export function bertin1953Raw() {
   var hammer = hammerRaw(1.68, 2),
       fu = 1.4, k = 12;
 
-  return function(lambda, phi) {
+  function forward(lambda, phi) {
 
     if (lambda + phi < -fu) {
       var u = (lambda - phi + 1.6) * (lambda + phi + fu) / 8;
@@ -28,7 +29,10 @@ export function bertin1953Raw() {
     }
 
     return r;
-  };
+  }
+
+  forward.invert = solve2d(forward);
+  return forward;
 }
 
 export default function() {

src/chamberlin.js

@@ -1,5 +1,6 @@
 import {geoCentroid as centroid, geoProjection as projection, geoRotation as rotation} from "d3-geo";
 import {abs, acos, asin, atan2, cos, epsilon, floor, pi, radians, sin, sqrt} from "./math.js";
+import {solve2d} from "./newton.js";
 
 // Azimuthal distance.
 function distance(dPhi, c1, s1, c2, s2, dLambda) {
@@ -103,7 +104,9 @@ export default function chamberlin(p0, p1, p2) { // TODO order matters!
   var c = centroid({type: "MultiPoint", coordinates: [p0, p1, p2]}),
       R = [-c[0], -c[1]],
       r = rotation(R),
-      p = projection(chamberlinRaw(pointRadians(r(p0)), pointRadians(r(p1)), pointRadians(r(p2)))).rotate(R),
+      f = chamberlinRaw(pointRadians(r(p0)), pointRadians(r(p1)), pointRadians(r(p2)));
+  f.invert = solve2d(f);
+  var p = projection(f).rotate(R),
       center = p.center;
 
   delete p.rotate;

src/interrupted/index.js

@@ -64,7 +64,7 @@ function interpolateSphere(lobes) {
   };
 }
 
-export default function(project, lobes) {
+export default function(project, lobes, inverse) {
   var sphere, bounds;
 
   function forward(lambda, phi) {
@@ -75,18 +75,21 @@ export default function(project, lobes) {
     return p;
   }
 
-  // Assumes mutually exclusive bounding boxes for lobes.
-  if (project.invert) forward.invert = function(x, y) {
-    var bound = bounds[+(y < 0)], lobe = lobes[+(y < 0)];
-    for (var i = 0, n = bound.length; i < n; ++i) {
-      var b = bound[i];
-      if (b[0][0] <= x && x < b[1][0] && b[0][1] <= y && y < b[1][1]) {
-        var p = project.invert(x - project(lobe[i][1][0], 0)[0], y);
-        p[0] += lobe[i][1][0];
-        return pointEqual(forward(p[0], p[1]), [x, y]) ? p : null;
+  if (inverse) {
+    forward.invert = inverse(forward);
+  } else if (project.invert) {
+    forward.invert = function(x, y) {
+      var bound = bounds[+(y < 0)], lobe = lobes[+(y < 0)];
+      for (var i = 0, n = bound.length; i < n; ++i) {
+        var b = bound[i];
+        if (b[0][0] <= x && x < b[1][0] && b[0][1] <= y && y < b[1][1]) {
+          var p = project.invert(x - project(lobe[i][1][0], 0)[0], y);
+          p[0] += lobe[i][1][0];
+          return pointEqual(forward(p[0], p[1]), [x, y]) ? p : null;
+        }
       }
-    }
-  };
+    };
+  }
 
   var p = projection(forward),
       stream_ = p.stream;

src/interrupted/sinuMollweide.js

@@ -1,5 +1,6 @@
 import {sinuMollweideRaw} from "../sinuMollweide.js";
 import interrupt from "./index.js";
+import {solve2d} from "../newton.js";
 
 var lobes = [[ // northern hemisphere
   [[-180,  35], [ -30,  90], [   0,  35]],
@@ -11,7 +12,7 @@ var lobes = [[ // northern hemisphere
 ]];
 
 export default function() {
-  return interrupt(sinuMollweideRaw, lobes)
+  return interrupt(sinuMollweideRaw, lobes, solve2d)
       .rotate([-20, -55])
       .scale(164.263)
       .center([0, -5.4036]);

src/newton.js

@@ -1,6 +1,6 @@
-import {abs, epsilon} from "./math.js";
+import {abs, epsilon, epsilon2} from "./math.js";
 
-// Newton-Raphson
+// Approximate Newton-Raphson
 // Solve f(x) = y, start from x
 export function solve(f, y, x) {
   var steps = 100, delta, f0, f1;
@@ -14,3 +14,38 @@ export function solve(f, y, x) {
   } while (steps-- > 0 && abs(delta) > epsilon);
   return steps < 0 ? NaN : x;
 }
+
+// Approximate Newton-Raphson in 2D
+// Solve f(a,b) = [x,y]
+export function solve2d(f, MAX_ITERATIONS = 40, eps = epsilon2) {
+  return function(x, y, a = 0, b = 0) {
+    var tx, ty;
+    for (var i = 0; i < MAX_ITERATIONS; i++) {
+      var p = f(a, b);
+      // diffs
+      tx = p[0] - x;
+      ty = p[1] - y;
+      if (abs(tx) < eps && abs(ty) < eps) break; // we're there!
+
+      // partial derivatives
+      var ea = (a > 0 ? -1 : 1) * eps,
+        eb = (b > 0 ? -1 : 1) * eps,
+        pa = f(a + ea, b),
+        pb = f(a, b + eb),
+        dxa = (pa[0] - p[0]) / ea,
+        dya = (pa[1] - p[1]) / ea,
+        dxb = (pb[0] - p[0]) / eb,
+        dyb = (pb[1] - p[1]) / eb,
+        // determinant
+        D = dyb * dxa - dya * dxb,
+        // newton step — or half-step for small D
+        l = (abs(D) < 0.5 ? 0.5 : 1) / D,
+        da = (ty * dxb - tx * dyb) * l,
+        db = (tx * dya - ty * dxa) * l;
+      a += da;
+      b += db;
+      if (abs(da) < eps && abs(db) < eps) break; // we're crawling
+    }
+    return [a, b];
+  };
+}

test/invert-test.js

@@ -12,13 +12,13 @@ var points = [[0, 0], [30.3, 24.1], [-10, 42], [-2, -5], [0,-55]];
   d3.geoAugust,
   d3.geoBaker,
   d3.geoBerghaus,
-  // d3.geoBertin1953,
+  d3.geoBertin1953,
   d3.geoBoggs,
   d3.geoBonne,
   d3.geoBottomley,
   d3.geoBromley,
-  // d3.geoChamberlin,
-  // d3.geoChamberlinAfrica,
+  // d3.geoChamberlin, // factory
+  d3.geoChamberlinAfrica,
   d3.geoCollignon,
   d3.geoCraig,
   d3.geoCraster,
@@ -55,7 +55,7 @@ var points = [[0, 0], [30.3, 24.1], [-10, 42], [-2, -5], [0,-55]];
   d3.geoInterruptedHomolosine,
   d3.geoInterruptedMollweide,
   d3.geoInterruptedMollweideHemispheres,
-  // d3.geoInterruptedSinuMollweide,
+  d3.geoInterruptedSinuMollweide,
   d3.geoInterruptedSinusoidal,
   d3.geoInterruptedQuarticAuthalic,
   d3.geoKavrayskiy7,
@@ -66,11 +66,11 @@ var points = [[0, 0], [30.3, 24.1], [-10, 42], [-2, -5], [0,-55]];
   d3.geoLoximuthal,
   d3.geoMiller,
   // d3.geoModifiedStereographic, // factory of factory
-  // d3.geoModifiedStereographicAlaska,
-  // d3.geoModifiedStereographicGs48,
-  // d3.geoModifiedStereographicGs50,
-  // d3.geoModifiedStereographicMiller,
-  // d3.geoModifiedStereographicLee,
+  // d3.geoModifiedStereographicAlaska, // see below
+  // d3.geoModifiedStereographicGs48, // see below
+  // d3.geoModifiedStereographicGs50, // see below
+  // d3.geoModifiedStereographicMiller, // see below
+  // d3.geoModifiedStereographicLee, // see below
   d3.geoMollweide,
   d3.geoMtFlatPolarParabolic,
   d3.geoMtFlatPolarQuartic,
@@ -90,7 +90,7 @@ var points = [[0, 0], [30.3, 24.1], [-10, 42], [-2, -5], [0,-55]];
   d3.geoSinusoidal,
   d3.geoTimes,
   // d3.geoTwoPointAzimuthal,
-  // d3.geoTwoPointAzimuthalUsa,
+  // d3.geoTwoPointAzimuthalUsa, // see below
   // d3.geoTwoPointEquidistant,
   d3.geoTwoPointEquidistantUsa,
   d3.geoVanDerGrinten,
@@ -112,3 +112,22 @@ var points = [[0, 0], [30.3, 24.1], [-10, 42], [-2, -5], [0,-55]];
     test.end();
   });
 });
+
+
+[
+{ factory: d3.geoModifiedStereographicAlaska, points: [[-149.9025632,61.2150138]] },
+{ factory: d3.geoModifiedStereographicGs48, points: [[-104.9833053, 39.7309179]] },
+{ factory: d3.geoModifiedStereographicGs50, points: [[-104.9833053, 39.7309179]] },
+{ factory: d3.geoModifiedStereographicMiller, points: [[0, 10]] },
+{ factory: d3.geoModifiedStereographicLee, points: [[179, 10]] },
+{ factory: d3.geoTwoPointAzimuthalUsa, points: [[-104.9833053, 39.7309179]] },
+].forEach(function(p) {
+  var factory = p.factory, points = p.points;
+  var name = factory.name, projection = factory();
+  tape(name + "(point) and " + name + ".invert(point) are symmetric YO", function(test) {
+    points.forEach(function(point) {
+      test.projectionEqual(projection, point, projection(point));
+    });
+    test.end();
+  });
+});