Cox Projection by fil

index.html

<!DOCTYPE html>
<canvas width="960" height="500"></canvas>
<script src="https://d3js.org/d3.v4.js"></script>
<script src="https://d3js.org/d3-geo-projection.v2.min.js"></script>
<script src="https://d3js.org/topojson.v2.min.js"></script>
<script src="https://unpkg.com/complex.js"></script>
<script src="geoInverseNewton.js"></script>
<script src="versor.js"></script>

<script>
  
  // Approximate \int _0 ^sm(z)  dt / (1 - t^3)^(2/3)
  // sm maps a triangle to a disc, sm^-1 does the opposite
  function sm_1(s){
    var z = Complex(s),
        w = Complex([-1/2, Math.sqrt(3)/2]),
        k = Complex(0);

    // rotate to have s ~= 1
    var rot = w.clone().pow(d3.scan([0,1,2].map(
      i => -(z.clone().mul(w.clone().pow(i))).re
    )));
      
    var y = rot.clone().mul(z).mul(-1).add(1);
    
    // var n = 3
    // w1 = gamma(1/n) * gamma(1 - 2/n) / n / gamma(1 - 1/n)
    // https://bl.ocks.org/Fil/852557838117687bbd985e4b38ff77d4
    var w1 = 1.7666387502854533;

    // McIlroy formula 5 p6 and table for F3 page 16
    var F0 = [
      1.44224957030741,
      0.240374928384568,
      0.0686785509670194,
      0.0178055502507087,
      0.00228276285265497,
      -1.48379585422573e-3,
      -1.64287728109203e-3,
      -1.02583417082273e-3,
      -4.83607537673571e-4,
      -1.67030822094781e-4,
      -2.45024395166263e-5,
      2.14092375450951e-5,
      2.55897270486771e-5,
      1.73086854400834e-5,
      8.72756299984649e-6,
      3.18304486798473e-6,
      4.79323894565283e-7
      -4.58968389565456e-7,
      -5.62970586787826e-7,
      -3.92135372833465e-7
    ];
    
    var F = Complex(0);
    for (var i = F0.length; i--;) {
      F = Complex(F0[i]).add(F.mul(y));
    }

    // /!\ mutates y
    var k = Complex(w1).add(y.pow(1-2/3).mul(-1).mul(F)).mul(rot.pow(2))

    
    // when we are close to [0,0] we switch to another approximation:
    // https://www.wolframalpha.com/input/?i=(-2%2F3+choose+k)++*+(-1)%5Ek++%2F+(k%2B1)+with+k%3D0,1,2,3,4
    // the difference is _very_ tiny but necessary
    // if we want projection(0,0) === [0,0]
    var n = z.abs(), m = 0.3;
    if (n < m) {
      var H0 = [
        1, 1/3, 5/27, 10/81, 22/243 //…
       ];
      var z3 = z.clone().pow(3);
      var h = Complex(0);
      for (i = H0.length; i--;) {
        h = Complex(H0[i]).add(h.mul(z3));
      }
      h = h.mul(z);
      k.mul(n / m).add(h.mul(1 - n / m));
    }
    
    return k.toVector();
  }
  

  var canvas = d3.select("canvas"),
  width = canvas.property("width"),
  height = canvas.property("height"),
  context = canvas.node().getContext("2d");

  // retina display
  var devicePixelRatio = window.devicePixelRatio || 1;
  canvas.style('width', canvas.attr('width')+'px');
  canvas.style('height', canvas.attr('height')+'px');
  canvas.attr('width', canvas.attr('width') * devicePixelRatio);
  canvas.attr('height', canvas.attr('height') * devicePixelRatio);
  context.scale(devicePixelRatio,devicePixelRatio);


  function coxRaw(lambda, phi){
    var s = d3.geoLagrangeRaw(0.5)(lambda, phi);
    var t = sm_1([s[1] / 2, s[0] / 2]);
    return [t[1], t[0]]
  }

  if (typeof geoInverse == 'function')
    coxRaw.invert = geoInverse(coxRaw);

projection = d3.geoProjection(coxRaw)
  .rotate([-10.8,0])
  .scale(188.68216)
  .translate([480, 500 * 2/3])
  .precision(0.01);

  
// the Sphere should go *exactly* to the vertices of the triangles
// because they are singular points
function sphere() {
  var degrees = 180 / Math.PI,
      c = 2 * Math.asin(1/Math.sqrt(5))*degrees;
  return {
    type: "Polygon",
    coordinates: [
      [ [ 0,90 ],  [ -180, -c ], [ 0, -90 ], [ 180, -c ], [ 0,90 ] ]
    ]
  };
}
  
var stream_ = projection.stream;
projection.stream = function(stream) {
    var rotate = projection.rotate(),
        rotateStream = stream_(stream),
        sphereStream = (projection.rotate([0, 0]), stream_(stream));
    projection.rotate(rotate);
    rotateStream.sphere = function() { d3.geoStream(sphere(), sphereStream); };
    return rotateStream;
  };
  
var path = d3.geoPath().projection(projection).context(context);

  canvas.call(d3.drag()
    .on("start", dragstarted)
    .on("drag", dragged));

var render = function() {},
    v0, // Mouse position in Cartesian coordinates at start of drag gesture.
    r0, // Projection rotation as Euler angles at start.
    q0; // Projection rotation as versor at start.

function dragstarted() {
  v0 = versor.cartesian(projection.invert(d3.mouse(this)));
  r0 = projection.rotate();
  q0 = versor(r0);
}

function dragged() {
  var v1 = versor.cartesian(projection.rotate(r0).invert(d3.mouse(this))),
      q1 = versor.multiply(q0, versor.delta(v0, v1)),
      r1 = versor.rotation(q1);
  projection.rotate(r1);
  render();
}

d3.json("https://unpkg.com/world-atlas@1/world/110m.json", function(
  error,
  world
) {
  if (error) throw error;

  var land = topojson.merge(world, world.objects.countries.geometries);

  render = function() {
    context.fillStyle = "#fff";
    context.fillRect(0, 0, width, height);

    context.beginPath();
    path({type:"Sphere"});
    context.strokeStyle = "black";
    context.lineWidth = 1.5;
    context.stroke(), context.closePath();

    context.beginPath();
    path(land);
    context.lineWidth = 1;
    context.strokeStyle = "#000";
    context.stroke();
    context.fillStyle = "#000";
    context.fill();
    context.closePath();
    
    context.beginPath();
    path(d3.geoGraticule()());
    context.strokeStyle = "#777";
    context.lineWidth = 0.5;
    context.stroke(), context.closePath();

  };

  render();
});

</script>

geoInverseNewton.js

// import {abs, epsilon, halfPi, sqrt} from 'math'
var abs = Math.abs, epsilon = 1e-6, halfPi = Math.PI/2, sqrt = Math.sqrt;

// import {cartesian, cartesianScale, spherical} from 'cartesian'
var asin  =Math.asin, atan2 = Math.atan2, cos = Math.cos, sin = Math.sin;

function spherical(cartesian) {
  return [atan2(cartesian[1], cartesian[0]), asin(cartesian[2])];
}

function cartesian(spherical) {
  var lambda = spherical[0], phi = spherical[1], cosPhi = cos(phi);
  return [cosPhi * cos(lambda), cosPhi * sin(lambda), sin(phi)];
}

function cartesianScale(vector, k) {
  return [vector[0] * k, vector[1] * k, vector[2] * k];
}

function initial_general(project) {
      var n = 50,
      step = (halfPi - epsilon) / n,
      i,
      j,
      grid = [];
    for (i = 0; i <= 4 * n; i++) {
      grid[i] = [];
      for (j = 0; j <= 2 * n; j++) {
        var p = [(i - 2 * n) * step, (j - n) * step]
        grid[i][j] = project(p[0], p[1]);
      }
    }

  return function(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) ];
      return c;
  }
}

function geoInverse(project, precision, initial) {
  
  _initial = initial_general(project)
  _____initial = function(x,y) {
    if (y > 0.3) return ___initial(x+0.16,y-0.035);
    if (x > 0) return ___initial(x+0.16,y+0.015);
    return ___initial(x-0.1,y);
  }
  function invert (x,y) {
    // find a start point c
    var point = [x,y],
        c = _initial(x,y);
        
    c[0] *= 0.999;
    
    c = cartesian(c);
    // solve for x,y
    try {
      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], 1/norm];
        },
        [point[0], point[1], 1],
        { start: c, acc: precision, dx: 1e-5, max: 100 }
      );
    } catch(e) {
    	console.log(e);
    }
    if (solution) 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 x /* bad approximation better that nothing ?? */;
        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),
            no = Newton.Norm(Dx)
            // 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];// * Math.min(no, 1e-2);
            }
//console.log(z, no)
            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;
  };

versor.js

// Version 0.0.0. Copyright 2017 Mike Bostock.
(function(global, factory) {
  typeof exports === 'object' && typeof module !== 'undefined' ? module.exports = factory() :
  typeof define === 'function' && define.amd ? define(factory) :
  (global.versor = factory());
}(this, (function() {'use strict';

var acos = Math.acos,
    asin = Math.asin,
    atan2 = Math.atan2,
    cos = Math.cos,
    max = Math.max,
    min = Math.min,
    PI = Math.PI,
    sin = Math.sin,
    sqrt = Math.sqrt,
    radians = PI / 180,
    degrees = 180 / PI;

// Returns the unit quaternion for the given Euler rotation angles [λ, φ, γ].
function versor(e) {
  var l = e[0] / 2 * radians, sl = sin(l), cl = cos(l), // λ / 2
      p = e[1] / 2 * radians, sp = sin(p), cp = cos(p), // φ / 2
      g = e[2] / 2 * radians, sg = sin(g), cg = cos(g); // γ / 2
  return [
    cl * cp * cg + sl * sp * sg,
    sl * cp * cg - cl * sp * sg,
    cl * sp * cg + sl * cp * sg,
    cl * cp * sg - sl * sp * cg
  ];
}

// Returns Cartesian coordinates [x, y, z] given spherical coordinates [λ, φ].
versor.cartesian = function(e) {
  var l = e[0] * radians, p = e[1] * radians;
  return [cos(p) * cos(l), cos(p) * sin(l), sin(p)];
};

// Returns the Euler rotation angles [λ, φ, γ] for the given quaternion.
versor.rotation = function(q) {
  return [
    atan2(2 * (q[0] * q[1] + q[2] * q[3]), 1 - 2 * (q[1] * q[1] + q[2] * q[2])) * degrees,
    asin(max(-1, min(1, 2 * (q[0] * q[2] - q[3] * q[1])))) * degrees,
    atan2(2 * (q[0] * q[3] + q[1] * q[2]), 1 - 2 * (q[2] * q[2] + q[3] * q[3])) * degrees
  ];
};

// Returns the quaternion to rotate between two cartesian points on the sphere.
versor.delta = function(v0, v1) {
  var w = cross(v0, v1), l = sqrt(dot(w, w));
  if (!l) return [1, 0, 0, 0];
  var t = acos(max(-1, min(1, dot(v0, v1)))) / 2, s = sin(t); // t = θ / 2
  return [cos(t), w[2] / l * s, -w[1] / l * s, w[0] / l * s];
};

// Returns the quaternion that represents q0 * q1.
versor.multiply = function(q0, q1) {
  return [
    q0[0] * q1[0] - q0[1] * q1[1] - q0[2] * q1[2] - q0[3] * q1[3],
    q0[1] * q1[0] + q0[0] * q1[1] + q0[2] * q1[3] - q0[3] * q1[2],
    q0[0] * q1[2] - q0[1] * q1[3] + q0[2] * q1[0] + q0[3] * q1[1],
    q0[0] * q1[3] + q0[1] * q1[2] - q0[2] * q1[1] + q0[3] * q1[0]
  ];
};

function cross(v0, v1) {
  return [
    v0[1] * v1[2] - v0[2] * v1[1],
    v0[2] * v1[0] - v0[0] * v1[2],
    v0[0] * v1[1] - v0[1] * v1[0]
  ];
}

function dot(v0, v1) {
  return v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2];
}

return versor;
})));