Clarify the resampleLineTo "midpoint close to an end" condition

[stbuehler]

d3-geo/src/projection/resample.js

Line 41 in d10a651

|| abs((dx * dx2 + dy * dy2) / d2 - 0.5) > 0.3 // midpoint close to an end

Where (dx, dy) = (x1-x0, y1-y0) (as vectors: d = p1 - p0), and d2 = dx*dx + dy*dy, i.e. d2 = square(norm(d)), and (dx2, dy2) = (x2-x0, y2-y0) (as vectors: d' = p2 - p0).

Now the "progress" of p2 on the line from p0 to p1 is (p2-p0)*(p1-p0)/norm(p1-p0) = d*d'/sqrt(d2) = (dx*dx2+dy*dy2)/sqrt(d2), but the code uses d2 instead of sqrt(d2).

Which raises another question: is the condition even useful to trigger another recursion?

[jrus]

Specifically what is being checked is the ratio of the dot products |v₂·v₁ / v₁² − 0.5| < 0.3, not the "distance" to the proposed midpoint per se. Where v₂ is the displacement vector from the start point to the proposed midpoint, and v₁ is the displacement vector from start point to end point.

Checking this ratio is equivalent to checking that the foot of the perpendicular from p₂ onto the segment p₀p₁ lies between 0.2–0.8 of the way along the segment.

If you want to play with how this works, the green areas in this interactive plot would be allowed (i.e. would stop further bisections) under this criterion: https://www.desmos.com/calculator/1gugvfmy9v

When you also add the primary area check, what you get for allowed midpoints is a segment-aligned rectangle whose length depends on the segment length, and whose width depends on the output resolution / specified tolerance.

---

That check was added by @mbostock in d3/d3@c7d4fab#diff-66cfd10649bb3fd0d57a09c64fcd22cd0f2ac091ef218ca04046417c3ab1b69f

Handle distortions near poles.

The threshhold was changed by @jasondavies from 0.4 to 0.3 in d3/d3@1a7810b#diff-66cfd10649bb3fd0d57a09c64fcd22cd0f2ac091ef218ca04046417c3ab1b69f

Adjust threshold for distortion resampling. Includes test case.

Then that part of the check was removed by @mbostock in d3/d3@6b3c198#diff-66cfd10649bb3fd0d57a09c64fcd22cd0f2ac091ef218ca04046417c3ab1b69f

Simplify the minimum sample distance check. Rather than compute the distance from the start to the new middle, use the already computed distance from the start to the end. Also, we no longer need the additional check that the middle is close to either the start or the end because this is covered by the minimum sample distance check.

It was subsequently added back by @jasondavies in d3/d3@3f8097d#diff-66cfd10649bb3fd0d57a09c64fcd22cd0f2ac091ef218ca04046417c3ab1b69f

Resample using a great-circle distance threshold. Projections will typically exhibit less curvature at larger scales, with geodesics becoming closer and closer to straight lines. Using a minimum angular distance for resampling therefore makes more sense than using an arbitrary scale-dependent threshold.

It's not clear whether this new threshold should depend on the precision parameter, or be configurable separately.

(Which is talking about the other heuristic condition added at the same time.)

---

It’s not clear from the code itself which specific cases the “midpoint close to an end” condition ever occurs or occurred. By sleuthing in the history we discover that it is targeted at fixing an issue that cropped up near the singularity at the pole for conic projections.

Presumably the specific thing that happens in such cases is the midpoint sometimes gets projected to very close to one end of the segment, and therefore looks like it is on the line, even though the derivatives (d(projection) / d(arclength)) is very steep at the midpoint. The resampling ends too soon, and you get part of the arc chopped straight. Adding this check eliminated the specific rendering bugs noticed.

This whole set of heuristics can of course still be fooled. You can get a projected S-shaped curve whose midpoint is near the midpoint of the chord, even though the curve itself diverges dramatically. These heuristics are a balance between performance and accuracy: They work in most cases, but they can also fail for various extreme map projections, and if you intentionally try to construct a failing case it isn’t too hard.

It would be nice if each of the heuristics came with specific test cases / explanations, citations to prior sources, etc., ideally even linked from comments in-line in the code.

---

Unfortunately when D3 v4 split the repos up, old git history became a bit harder to track across the split. If you want the full history of resample.js you have to look in a few places:

June 2016 – present:
https://github.com/d3/d3-geo/commits/a5558b83343702d2a5e352e19ff107dc7e71199c/src/projection/resample.js

December 2012 – November 2015:
https://github.com/d3/d3/commits/d05692b4c7629a05eb178a70cdfa669694ceca4c/src/geo/resample.js

September 2012 – December 2012 (this list also includes changes to other functions):
https://github.com/d3/d3/commits/6091249215e8687544482460ff83dcc76563a6b4/src/geo/projection.js

---

Here is the 26 September 2012 commit where @jasondavies first added adaptive resampling: d3/d3@45d4b31#diff-0dfd24d6cf92b0917f4bee9fcd9490230dfa74c2fd9d347b78be78cd7e067758

Adaptive resampling. Can still be optimised further by caching trig. results.

[jrus]

In case LaTeX is clearer:

[stbuehler]

Uff.. obviously this was a long time ago for me. I probably should have also written down what the original problem was (i.e. why I was looking at the code in the first place), because I don't remember...

I think you are right: what I named "progress" in my first post is not the "relative" position (from 0 to 1), but the "length" on the segment p₀p₁. To get the relative progress (d·d')/(d·d) seems to be the correct formula. No idea why I thought different 2 years ago.

I'm quite sure I'm not qualified to comment on anything else you said, but thanks anyway :)

Want me to close this?

[jrus]

I went to the trouble to sleuth in the history because I also was curious about how these heuristic criteria came about, and found the code lacking inline explanation.

Once we understand how the criteria work, one thing we can try to do is implement the inverse form of adaptive resampling: starting from a straight segment on a map, we can inverse-project to get a chain of great-circle arcs on the sphere. That’s what I started writing up at https://observablehq.com/d/cf5950e1f8bc3070; explaining the planar d3-geo heuristic is a means to an end.

(This is useful for example for inverse-projecting global grids onto the sphere, or for inverse-projecting the shape of a flat raster map. Or could be helpful when projecting a geojson file whose author was strictly following the geojson spec and anticipating linear interpolation in terms of an equirectangular projection, or for projecting political boundaries defined to lie along a latitude line or property boundaries defined in a particular planar coordinate system. Or could be helpful for approximating an arbitrary small circle on the sphere e.g. as part of the graticule, or approximating a rhumb line via inverse Mercator projection.)

[jrus]

@mbostock How attached are you to the current test for stopping bisection? My suspicion is that an ellipse will generally make for a better test than a rectangle, and should cost about the same, especially since we are already squaring everything in the current version to avoid square roots from calculating distances.

Here are some possible alternative criteria, in a Desmos plot: https://www.desmos.com/calculator/ntnlzzutn0

Any time the projected midpoint is outside the shape, we trigger an additional bisection.

Red rectangle = current D3-geo test.
Green lens = angle-based test: the width at the midpoint defines a maximum turning angle for the intermediate point (the locus of points of some fixed turning angle relative to some fixed chord is a circle)
Purple ellipse = make the current test into an ellipse by combining the width/height tests (more or less turn a test of the form a² < c² && b² < c² into a test of the form a² + b² < c²), but change the width constant to widen it a bit, because we no longer have to worry about the corners.

To elaborate about the issue with corners that stick out past the lens shape: imagine our projected segment has roughly constant curvature on the map (locally looks like a circular arc), but is not parametrized in an arclength-uniform way. With the current rectangle we could end up with a circular arc which sticks out to a bit over 1.5 times our δ threshold, but parametrized so that its projected midpoint lands in the corner area, and we will stop bisecting.

It may also be worth including one extra test for minimum distance. That is, any time the segment is longer than, say, 16 times the tolerance (or could be 100x or ...), always bisect even if the projected midpoint is right in the center. This should generally cause few extra bisections compared to not having it, but should prevent the most egregious current failures where some very large large visible curve gets rendered straight because the projected midpoint happened to fall close to the chord center. (Other possible hacks to prevent such behaviors will generally require using more samples across the board, which could bog things down for very expensive projections.)

---

It would be nice if there were (a) a good (ideally already implemented) error metric for total error between two paths, and (b) some kind of benchmark example collections of points and map projections, against which we could test a few possibilities to see what kind of "mistakes" they make, how many segments, what the total error is, etc., as the tolerance parameter gets tweaked.