A new quincuncial projection

[theohonohan]

This is obscure, but might be interesting to implement.

https://mpetroff.net/2021/08/square-equal-area-map-projection/

https://dl.acm.org/doi/10.1145/3460521

From the paper:

"[T]o improve upon existing square equal-area projections, a new projection must combine the low angular distortion of the Gringorten projection and its reduced occurrence of significant visible cusps with the Collignon quincuncial projection’s advantage of having closed-form solutions. These are the properties of the new projection presented in this manuscript."

[mbostock]

I ported the implementation to a notebook: https://observablehq.com/@d3/petroff-quincuncial

[jrus]

You may also find https://observablehq.com/@jrus/sac-quincuncial interesting, though that is properly intended as an intermediate step toward a plot on equilateral triangles.

(The spherical area coordinates version should actually be better for storing panorama textures, as long as the pixels are interpreted as hexagons, with some care at octant edges.)

---

Edit: note this is not quite equal-area. The conceptual definition of this map is: for any point X in the octant ABC, the ratio of areas of the three spherical triangles ABX:BCX:CAX is the same as the ratio of areas of the three planar triangles ABX:BCX:CAX.

Because in stereographic projection, the locus of vertices of equal-area triangles for a given base with one vertex at the origin is a straight line (see https://link.springer.com/article/10.1007/s13366-018-0426-2), this ends up being a particularly simple formula to write down.

This map was first proposed in Praun / Hoppe (2003) http://hhoppe.com/sphereparam.pdf

Later in:
Carfora (2007) https://doi.org/10.1016/j.cam.2006.10.068
Lei / Qi / Tian (2020) https://www.mdpi.com/2076-3417/10/2/655/htm

---

I was working on this about 5–6 months ago, but then got stalled / distracted by other stuff.

I wrote up a bit about how to round points to the nearest hexagon in https://observablehq.com/@jrus/hexround (but this does not discuss the proper corrections near octant boundaries, which I worked out on paper but never wrote up or implemented in code)

Here’s a 32x32 "pixel" grid on the sphere (which can be made denser for any 2n x 2n size; I recommend powers of 2): https://observablehq.com/@jrus/sphere-resample

[mbostock]

That’s super cool, @jrus!

[Fil]

petroff would need an inverse—but it seems a lot of work compared to Jacob’s :)

[mpetroff]

@Fil My projection has a closed-form inverse (that was part of the design requirements), but I never ported the Python implementation of the inverse to JavaScript, since I didn't need the inverse for making the paper figures, and the analysis was done in Python. Porting the forward projection implementation was pretty straightforward, so I wouldn't expect any difficulties with porting the inverse either.

[Fil]

What am I doing wrong? For me the JS implementation doesn't return the same values as the python implementation:

python:

print (new_projection(0.001, 0.001)) # (0.0005127294198217123, -0.99897454084992399)

JS:

projectionRaw(0.001, 0.001) // [0.0008026311120535786, 0.0008026315086362794]

[mpetroff]

There are layout and scaling differences between the two implementations. The Python implementation directly returns coordinates in a quincuncial arrangement with $x, y \in [-1, 1]$. The JavaScript implementation returns coordinates appropriate for passing to d3.geoQuincuncial, like is done for d3.geoGringortenRaw.

I probably shouldn't have recommended porting the Python inverse, since it's a branchless fully-vectorized implementation. Such conditions were necessary for allowing for auto-differentiation (or for a GPU implementation), but they're not ideal for a JavaScript implementation, since those restrictions aren't necessary for JavaScript and make the code more difficult to follow. I'll try to get to writing a JavaScript implementation of the inverse projection later this week, but it's probably best to base it off Appendix A of the paper instead of the Python implementation.

[mpetroff]

I implemented the inverse in JavaScript (also CC0): https://bl.ocks.org/mpetroff/dcf7b090eabda85d081d47e8a0c71d4a

Additionally, I updated the forward projection implementation (the hexadecant() function) to better match the paper. The copy in the supplemental information (which @mbostock's notebook is based on) is mathematically equivalent but was written prior to some changes made during the paper's review process, which clarified the derivation.

Unfortunately, there's a sign error in Appendix A of the paper, which I discovered while working on the new implementation of the inverse. In the second case of Equation (70), $\beta - \psi_0$ should read $\beta + \psi_0$ (the Python implementation was correct). The sentence after Equation (54) is also misleading and should read:

The fixed sub-triangle lengths $c$, $G$, $G'$, $F$, $a'$, and $c'$ are defined by Equations (17), (18), (19), (20), (21), and (22), respectively, except with $\theta$ replaced with $\theta'$, $\psi_0$ replaced with $\psi_0'$, and $\psi_1$ replaced with $\psi_1'$, in the case condition statements only.

I'm going to fix these mistakes in the arXiv copy and submit a corrigendum to ACM (in a couple weeks, in case other issues are found).

Update (2021-12-22): There's also a mistake in Equation (47). The 0 case should be $x_m >= 0, y_m <= 0$.

[Fil]

Fantastic! Here's a direct copy of your code in observable, linked to a world-map function which uses the projection's inverse to drag the globe with the mouse pointer: https://observablehq.com/@fil/petroff-quincuncial-w-inverse

(It seems to work perfectly)

[jrus]

The round trip through the inverse still has a few issues at edges (I should maybe try to understand and then rewrite the code for reflecting and unreflecting inputs/outputs, instead of just adding kludges on top), but you can still roughly see how a square "pixel" grid in the plane projects onto the sphere:

https://observablehq.com/d/331ce3a332c89c15