Skip to content

Commit 79dbf56

Browse files
mblochmbloch
committed
Switch to better-performing robust segment geometry
1 parent 41a6dd3 commit 79dbf56

12 files changed

+576
-164
lines changed

‎src/geojson/geojson-export.mjs‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,7 @@ GeoJSON.exportPointGeom = function(points, arcs) {
255255
return geom;
256256
};
257257

258+
// ids: shape data (array of array of arc ids)
258259
GeoJSON.exportLineGeom = function(ids, arcs) {
259260
var obj = exportPathData(ids, arcs, "polyline");
260261
if (obj.pointCount === 0) return null;
@@ -270,6 +271,7 @@ GeoJSON.exportLineGeom = function(ids, arcs) {
270271
};
271272
};
272273

274+
// ids: shape data (array of array of arc ids)
273275
GeoJSON.exportPolygonGeom = function(ids, arcs, opts) {
274276
var obj = exportPathData(ids, arcs, "polygon");
275277
if (obj.pointCount === 0) return null;

‎src/geom/mapshaper-bigint-utils.mjs‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
2+
export function toScaledStr(num, k) {
3+
var abs = Math.abs(num);
4+
var signed = num != abs;
5+
// String() matches behavior of big.js
6+
// (as opposed to toFixed() or toPrecision())
7+
var s = abs < 1e-6 ? abs.toFixed(k) : String(abs);
8+
var parts = s.split('.');
9+
var integer = parts[0] == '0' ? '' : parts[0];
10+
var decimal = parts.length == 2 ? parts[1] : '';
11+
if (decimal.length < k) {
12+
decimal = decimal.padEnd(k, '0');
13+
} else if (decimal.length > k) {
14+
decimal = decimal.slice(0, k);
15+
}
16+
var s2 = integer + decimal;
17+
// remove any leading 0s
18+
while (s2[0] == 0 && s2.length > 1) {
19+
s2 = s2.slice(1);
20+
}
21+
if (signed) {
22+
s2 = '-' + s2;
23+
}
24+
return s2;
25+
}
26+
27+
// returns Number
28+
export function fromScaledStr(s, decimals) {
29+
var signed = s[0] == '-';
30+
var uns = signed ? s.slice(1) : s;
31+
var s2, len;
32+
len = uns.length;
33+
if (len > decimals) {
34+
s2 = uns.slice(0, len - decimals) + '.' + uns.slice(-decimals);
35+
} else if (len == decimals) {
36+
s2 = '0.' + uns;
37+
} else {
38+
s2 = '0.' + uns.padStart(decimals, '0');
39+
}
40+
if (signed) {
41+
s2 = '-' + s2;
42+
}
43+
return Number(s2);
44+
}
45+
46+
export function findBigIntScaleFactor() {
47+
var minVal = Infinity;
48+
var intLen = 0, s;
49+
for (var i=0, n=arguments.length; i<n; i++) {
50+
minVal = Math.min(minVal, Math.abs(arguments[i]));
51+
}
52+
if (minVal >= 1) {
53+
s = minVal.toFixed(1);
54+
intLen = s.indexOf('.');
55+
} else if (minVal !== 0) {
56+
s = minVal.toFixed(10).slice(2); // decimal part, up to _ decimals
57+
while (s[0] === '0') {
58+
intLen--;
59+
s = s.slice(1);
60+
}
61+
}
62+
return 17 - intLen;
63+
}

‎src/geom/mapshaper-segment-geom-big.mjs‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
import Big from 'big.js';
2+
import { debug } from '../utils/mapshaper-logging';
3+
4+
// Find the intersection point of two segments that cross each other,
5+
// or return null if the segments do not cross.
6+
// Assumes endpoint intersections have already been detected
7+
export function findCrossIntersection_big(ax, ay, bx, by, cx, cy, dx, dy) {
8+
var ax_big = Big(ax);
9+
var ay_big = Big(ay);
10+
var bx_big = Big(bx);
11+
var by_big = Big(by);
12+
var cx_big = Big(cx);
13+
var cy_big = Big(cy);
14+
var dx_big = Big(dx);
15+
var dy_big = Big(dy);
16+
var v1x = bx_big.minus(ax_big);
17+
var v1y = by_big.minus(ay_big);
18+
var den_big = determinant2D_big(v1x, v1y, dx_big.minus(cx_big), dy_big.minus(cy_big));
19+
if (den_big.eq(0)) {
20+
debug('DIV0 error - should have been identified as collinear "touch" intersection.');
21+
// console.log("hit?", segmentHit_big(ax, ay, bx, by, cx, cy, dx, dy))
22+
// console.log('Seg 1', getSegFeature(ax, ay, bx, by, true))
23+
// console.log('Seg 2', getSegFeature(cx, cy, dx, dy, false))
24+
return null;
25+
}
26+
// perform division using regular math, which does not reduce overall
27+
// precision in test data (big.js division is very slow)
28+
// tests show identical result to:
29+
// var m_big = orient2D_big(cx_big, cy_big, dx_big, dy_big, ax_big, ay_big).div(den_big)
30+
var m = orient2D_big(cx_big, cy_big, dx_big, dy_big, ax_big,
31+
ay_big).toNumber() / den_big.toNumber();
32+
var m_big = Big(m);
33+
// console.log("big m:", m_big.toString())
34+
var x_big = ax_big.plus(m_big.times(v1x).round(16));
35+
var y_big = ay_big.plus(m_big.times(v1y).round(16));
36+
var p = [x_big.toNumber(), y_big.toNumber()];
37+
return p;
38+
}
39+
40+
export function orient2D_big(ax, ay, bx, by, cx, cy) {
41+
var a = (ax.minus(cx)).times(by.minus(cy));
42+
var b = (ay.minus(cy)).times(bx.minus(cx));
43+
return a.minus(b);
44+
}
45+
46+
export function orient2D_big2(ax, ay, bx, by, cx, cy) {
47+
return orient2D_big(Big(ax), Big(ay), Big(bx), Big(by), Big(cx), Big(cy));
48+
}
49+
50+
export function segmentHit_big(ax, ay, bx, by, cx, cy, dx, dy) {
51+
return orient2D_big(ax, ay, bx, by, cx, cy).times(
52+
orient2D_big(ax, ay, bx, by, dx, dy)).lte(0) &&
53+
orient2D_big(cx, cy, dx, dy, ax, ay).times(
54+
orient2D_big(cx, cy, dx, dy, bx, by)).lte(0);
55+
}
56+
57+
export function segmentHit_big2(ax, ay, bx, by, cx, cy, dx, dy) {
58+
return segmentHit_big(Big(ax), Big(ay), Big(bx), Big(by),
59+
Big(cx), Big(cy), Big(dx), Big(dy));
60+
}
61+
62+
function determinant2D_big(a, b, c, d) {
63+
return a.times(d).minus(b.times(c));
64+
}

‎src/geom/mapshaper-segment-geom.mjs‎

Lines changed: 70 additions & 90 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,8 @@
22
import { getHighPrecisionSnapInterval } from '../paths/mapshaper-snapping';
33
import { debug } from '../utils/mapshaper-logging';
44
import { distance2D, distanceSq, pointSegDistSq2 } from '../geom/mapshaper-basic-geom';
5-
import Big from "big.js";
5+
import { fromScaledStr, toScaledStr, findBigIntScaleFactor } from '../geom/mapshaper-bigint-utils';
6+
//import { findCrossIntersection_big } from '../geom/mapshaper-segment-geom-big';
67

78
// Find the intersection between two 2D segments
89
// Returns 0, 1 or 2 [x, y] locations as null, [x, y], or [x1, y1, x2, y2]
@@ -38,18 +39,31 @@ export function segmentIntersection(ax, ay, bx, by, cx, cy, dx, dy, epsArg) {
3839
return touches || cross || null;
3940
}
4041

41-
4242
function findCrossIntersection(ax, ay, bx, by, cx, cy, dx, dy, eps) {
4343
var p;
44-
if (eps > 0 && !segmentHit_fast(ax, ay, bx, by, cx, cy, dx, dy)) return null;
45-
else if (eps === 0 && !segmentHit_big(ax, ay, bx, by, cx, cy, dx, dy)) return null;
44+
// The normal-precision hit function works for all inputs when eps > 0 because
45+
// the geometries that cause the ordinary function fails are detected as
46+
// 'touches' or endpoint hits (at least this was true in all the real-world
47+
// data samples that were tested).
48+
//
49+
if (eps > 0 && !segmentHit_fast(ax, ay, bx, by, cx, cy, dx, dy)) {
50+
return null;
51+
} else if (eps === 0 && !segmentHit_robust(ax, ay, bx, by, cx, cy, dx, dy)) {
52+
return null;
53+
}
4654

47-
// in a typical layer with many intersections, robust is preferred in
48-
// most (>90%) segment intersections in order to keep the positional
49-
// error within a small interval (e.g. 50% of eps)
55+
// in a typical layer with many intersections along shared polygon boundaries,
56+
// robust is preferred in most (>90%) segment intersections in order to keep
57+
// the positional error within a small interval (e.g. 50% of eps)
5058
//
5159
if (useRobustCross(ax, ay, bx, by, cx, cy, dx, dy)) {
5260
p = findCrossIntersection_robust(ax, ay, bx, by, cx, cy, dx, dy);
61+
// var p2 = findCrossIntersection_big(ax, ay, bx, by, cx, cy, dx, dy);
62+
// var dx = p[0] - p2[0];
63+
// var dy = p[1] - p2[1];
64+
// if (dx != 0 || dy != 0) {
65+
// console.log(dx, dy)
66+
// }
5367
} else {
5468
p = findCrossIntersection_fast(ax, ay, bx, by, cx, cy, dx, dy);
5569
}
@@ -69,41 +83,6 @@ function findCrossIntersection(ax, ay, bx, by, cx, cy, dx, dy, eps) {
6983
return p;
7084
}
7185

72-
// Find the intersection point of two segments that cross each other,
73-
// or return null if the segments do not cross.
74-
// Assumes endpoint intersections have already been detected
75-
function findCrossIntersection_robust(ax, ay, bx, by, cx, cy, dx, dy) {
76-
var ax_big = Big(ax);
77-
var ay_big = Big(ay);
78-
var bx_big = Big(bx);
79-
var by_big = Big(by);
80-
var cx_big = Big(cx);
81-
var cy_big = Big(cy);
82-
var dx_big = Big(dx);
83-
var dy_big = Big(dy);
84-
var v1x = bx_big.minus(ax_big);
85-
var v1y = by_big.minus(ay_big);
86-
var den_big = determinant2D_big(v1x, v1y, dx_big.minus(cx), dy_big.minus(cy));
87-
if (den_big.eq(0)) {
88-
debug("DIV0 error (should have been caught upstream)");
89-
// console.log("hit?", segmentHit_big(ax, ay, bx, by, cx, cy, dx, dy))
90-
// console.log('Seg 1', getSegFeature(ax, ay, bx, by, true))
91-
// console.log('Seg 2', getSegFeature(cx, cy, dx, dy, false))
92-
return null;
93-
}
94-
// perform division using regular math, which does not reduce overall
95-
// precision in test data (big.js division is very slow)
96-
// tests show identical result to:
97-
// orient2D_big(cx_big, cy_big, dx_big, dy_big, ax_big, ay_big).div(den_big)
98-
var m = orient2D_big(cx_big, cy_big, dx_big, dy_big, ax_big,
99-
ay_big).toNumber() / den_big.toNumber();
100-
var m_big = Big(m);
101-
var x_big = ax_big.plus(m_big.times(v1x).round(16));
102-
var y_big = ay_big.plus(m_big.times(v1y).round(16));
103-
var p = [x_big.toNumber(), y_big.toNumber()];
104-
return p;
105-
}
106-
10786
function findCrossIntersection_fast(ax, ay, bx, by, cx, cy, dx, dy) {
10887
var den = determinant2D(bx - ax, by - ay, dx - cx, dy - cy);
10988
var m = orient2D(cx, cy, dx, dy, ax, ay) / den;
@@ -118,6 +97,33 @@ function findCrossIntersection_fast(ax, ay, bx, by, cx, cy, dx, dy) {
11897
return p;
11998
}
12099

100+
// this function, using BigInt, is 3-4x faster than the version using big.js
101+
export function findCrossIntersection_robust(ax, ay, bx, by, cx, cy, dx, dy) {
102+
var d = findBigIntScaleFactor(ax, ay, bx, by, cx, cy, dx, dy);
103+
var d2 = 16; // scale numerator of integer division by this many decimal digits
104+
var k_bi = 10000000000000000n; // matches d2
105+
var ax_bi = BigInt(toScaledStr(ax, d));
106+
var ay_bi = BigInt(toScaledStr(ay, d));
107+
var bx_bi = BigInt(toScaledStr(bx, d));
108+
var by_bi = BigInt(toScaledStr(by, d));
109+
var cx_bi = BigInt(toScaledStr(cx, d));
110+
var cy_bi = BigInt(toScaledStr(cy, d));
111+
var dx_bi = BigInt(toScaledStr(dx, d));
112+
var dy_bi = BigInt(toScaledStr(dy, d));
113+
var den = determinant2D(bx_bi - ax_bi, by_bi - ay_bi, dx_bi - cx_bi, dy_bi - cy_bi);
114+
if (den === 0n) {
115+
debug('DIV0 error - should have been identified as collinear "touch" intersection.');
116+
return null;
117+
}
118+
var num = orient2D(cx_bi, cy_bi, dx_bi, dy_bi, ax_bi, ay_bi) * k_bi;
119+
var m_bi = num / den;
120+
var x_bi = ax_bi * k_bi + m_bi * (bx_bi - ax_bi);
121+
var y_bi = ay_bi * k_bi + m_bi * (by_bi - ay_bi);
122+
var x = fromScaledStr(x_bi.toString(), d + d2);
123+
var y = fromScaledStr(y_bi.toString(), d + d2);
124+
return [x, y];
125+
}
126+
121127
function useRobustCross(ax, ay, bx, by, cx, cy, dx, dy) {
122128
// angle and seg length ratio thresholds were found by comparing
123129
// fast and robust outputs on sample data
@@ -129,7 +135,7 @@ function useRobustCross(ax, ay, bx, by, cx, cy, dx, dy) {
129135
return false;
130136
}
131137

132-
// Returns smaller unsigned angle between two segments
138+
// Returns smaller of two angles between two segments (unsigned)
133139
function innerAngle(ax, ay, bx, by, cx, cy, dx, dy) {
134140
var v1x = bx - ax;
135141
var v1y = by - ay;
@@ -278,18 +284,6 @@ function clampIntersectionPoint(p, ax, ay, bx, by, cx, cy, dx, dy) {
278284
p[1] = y;
279285
}
280286

281-
// function getSegFeature(x1, y1, x2, y2, hot) {
282-
// return JSON.stringify({
283-
// type: "Feature",
284-
// properties: {
285-
// stroke: hot ? "orange" : "blue"
286-
// },
287-
// geometry: {
288-
// type: "LineString",
289-
// coordinates: [[x1, y1], [x2, y2]]
290-
// }
291-
// });
292-
// }
293287

294288
// a: coordinate of point
295289
// b: endpoint coordinate of segment
@@ -322,10 +316,6 @@ function determinant2D(a, b, c, d) {
322316
return a * d - b * c;
323317
}
324318

325-
function determinant2D_big(a, b, c, d) {
326-
return a.times(d).minus(b.times(c));
327-
}
328-
329319
// returns a positive value if the points a, b, and c are arranged in
330320
// counterclockwise order, a negative value if the points are in clockwise
331321
// order, and zero if the points are collinear.
@@ -334,22 +324,31 @@ export function orient2D(ax, ay, bx, by, cx, cy) {
334324
return determinant2D(ax - cx, ay - cy, bx - cx, by - cy);
335325
}
336326

337-
export function orient2D_big(ax, ay, bx, by, cx, cy) {
338-
var a = (ax.minus(cx)).times(by.minus(cy));
339-
var b = (ay.minus(cy)).times(bx.minus(cx));
340-
return a.minus(b);
327+
export function orient2D_robust(ax, ay, bx, by, cx, cy) {
328+
var d = findBigIntScaleFactor(ax, ay, bx, by, cx, cy);
329+
var ax_bi = BigInt(toScaledStr(ax, d));
330+
var ay_bi = BigInt(toScaledStr(ay, d));
331+
var bx_bi = BigInt(toScaledStr(bx, d));
332+
var by_bi = BigInt(toScaledStr(by, d));
333+
var cx_bi = BigInt(toScaledStr(cx, d));
334+
var cy_bi = BigInt(toScaledStr(cy, d));
335+
var o2d_bi = orient2D(ax_bi, ay_bi, bx_bi, by_bi, cx_bi, cy_bi);
336+
return fromScaledStr(o2d_bi.toString(), d);
341337
}
342338

343-
export function orient2D_big2(ax, ay, bx, by, cx, cy) {
344-
return orient2D_big(Big(ax), Big(ay), Big(bx), Big(by), Big(cx), Big(cy));
339+
function segmentHit_robust(ax, ay, bx, by, cx, cy, dx, dy) {
340+
var d = findBigIntScaleFactor(ax, ay, bx, by, cx, cy, dx, dy);
341+
var ax_bi = BigInt(toScaledStr(ax, d));
342+
var ay_bi = BigInt(toScaledStr(ay, d));
343+
var bx_bi = BigInt(toScaledStr(bx, d));
344+
var by_bi = BigInt(toScaledStr(by, d));
345+
var cx_bi = BigInt(toScaledStr(cx, d));
346+
var cy_bi = BigInt(toScaledStr(cy, d));
347+
var dx_bi = BigInt(toScaledStr(dx, d));
348+
var dy_bi = BigInt(toScaledStr(dy, d));
349+
return segmentHit_fast(ax_bi, ay_bi, bx_bi, by_bi, cx_bi, cy_bi, dx_bi, dy_bi);
345350
}
346351

347-
348-
// export function orient2D_v2(ax, ay, bx, by, cx, cy) {
349-
// return -orient2D_robust(ax, ay, bx, by, cx, cy);
350-
// }
351-
352-
353352
// Source: Sedgewick, _Algorithms in C_
354353
// (Other functions were tried that were more sensitive to floating point errors
355354
// than this function)
@@ -360,25 +359,6 @@ export function segmentHit_fast(ax, ay, bx, by, cx, cy, dx, dy) {
360359
orient2D(cx, cy, dx, dy, bx, by) <= 0;
361360
}
362361

363-
export function segmentHit_big(ax, ay, bx, by, cx, cy, dx, dy) {
364-
return orient2D_big(ax, ay, bx, by, cx, cy).times(
365-
orient2D_big(ax, ay, bx, by, dx, dy)).lte(0) &&
366-
orient2D_big(cx, cy, dx, dy, ax, ay).times(
367-
orient2D_big(cx, cy, dx, dy, bx, by)).lte(0);
368-
}
369-
370-
export function segmentHit_big2(ax, ay, bx, by, cx, cy, dx, dy) {
371-
return segmentHit_big(Big(ax), Big(ay), Big(bx), Big(by),
372-
Big(cx), Big(cy), Big(dx), Big(dy));
373-
}
374-
375-
// export function segmentHit_robust(ax, ay, bx, by, cx, cy, dx, dy) {
376-
// return -orient2D_robust(ax, ay, bx, by, cx, cy) *
377-
// -orient2D_robust(ax, ay, bx, by, dx, dy) <= 0 &&
378-
// -orient2D_robust(cx, cy, dx, dy, ax, ay) *
379-
// -orient2D_robust(cx, cy, dx, dy, bx, by) <= 0;
380-
// }
381-
382362
// Useful for determining if a segment that intersects another segment is
383363
// entering or leaving an enclosed buffer area
384364
// returns -1 if angle of p1p2 -> p3p4 is counter-clockwise (left turn)

0 commit comments

Comments
 (0)