Commit 4bdd547a authored by Matthias Betz's avatar Matthias Betz
Browse files

parse coordinates without split/parseFloat, triangulate simple quads directly



- coordinates.js: scans posList/pos text directly; plain decimals up to 15
  digits are converted exactly (integer / power of ten, bit-identical to
  parseFloat), everything else goes through parseFloat
- triangulate.js (moved out of citygml.js): triangles and quadrilaterals
  without holes are split directly, concave ones at the reflex corner;
  self-intersecting, degenerate or holed polygons and corners that coincide
  in libtess's 2D projection stay with libtess

On 156 test files vertex counts and colours are identical to before; total
area differs by at most 1.7e-6 (other diagonal of non-planar quads).
274 MB LoD2 file: 7.2 s -> 5.8 s in Node.

Co-Authored-By: default avatarClaude Opus 5.5 (1M context) <noreply@anthropic.com>
parent f50dc6f6
// CityGML parsing and mesh building, independent of the DOM and WebGL.
import { libtess, glMatrix } from '../vendor.js';
import { glMatrix } from '../vendor.js';
import { createTokenizer } from './xmlTokenizer.js';
import { parseNumbers } from './coordinates.js';
import { triangulate } from './triangulate.js';
var vec3 = glMatrix.vec3;
......@@ -191,16 +193,12 @@ function createParser(opts) {
// Parses whitespace separated x y z triples into [[x, y, z], ...],
// relative to the origin.
function parseRing(coordinateString) {
var trimmed = coordinateString.trim();
if (trimmed === "") {
return [];
}
var split = trimmed.split(/\s+/);
var numbers = parseNumbers(coordinateString);
var coords = [];
for (var i = 0; i + 2 < split.length; i = i + 3) {
var x = parseFloat(split[i]);
var y = parseFloat(split[i + 1]);
var z = parseFloat(split[i + 2]);
for (var i = 0; i + 2 < numbers.length; i = i + 3) {
var x = numbers[i];
var y = numbers[i + 1];
var z = numbers[i + 2];
if (origin === null) {
origin = [x, y, z];
}
......@@ -474,95 +472,4 @@ function readFile(file, parser, opts) {
});
}
function triangulate(polygon) {
var triangleVerts = [];
tessy.gluTessBeginPolygon(triangleVerts);
var exterior = polygon[0];
var normal = calculateNormal(polygon[0]);
tessy.gluTessNormal(normal[0], normal[1], normal[2]);
tessy.gluTessBeginContour();
for (var j = 0; j < exterior.length; j++) {
tessy.gluTessVertex(exterior[j], exterior[j]);
}
tessy.gluTessEndContour();
for (var i = 1; i < polygon.length; i++) {
tessy.gluTessBeginContour();
var contour = polygon[i];
for (var k = 0; k < contour.length; k++) {
tessy.gluTessVertex(contour[k], contour[k]);
}
tessy.gluTessEndContour();
}
tessy.gluTessEndPolygon();
return {
vertices: triangleVerts,
normal: normal
};
}
function calculateNormal(ring) {
var coords = [0, 0, 0];
for (var i = 0; i < ring.length - 1; i++) {
var current = ring[i + 0];
var next = ring[i + 1];
coords[0] += (current[2] + next[2]) * (current[1] - next[1]);
coords[1] += (current[0] + next[0]) * (current[2] - next[2]);
coords[2] += (current[1] + next[1]) * (current[0] - next[0]);
}
if (coords[0] == 0 && coords[1] == 0 && coords[2] == 0) {
// no valid normal vector found
if (ring.length < 3) {
// no three points, return x-axis
return vec3.fromValues(1, 0, 0);
}
// vec3.create() takes no arguments; clone copies the points.
return calculateNormalWithCross(vec3.clone(ring[0]), vec3.clone(ring[1]), vec3.clone(ring[2]));
}
var v = vec3.fromValues(coords[0], coords[1], coords[2]);
vec3.normalize(v, v);
return v;
}
function calculateNormalWithCross(v1, v2, v3) {
var dir1 = vec3.create();
vec3.sub(dir1, v2, v1);
var dir2 = vec3.create();
vec3.sub(dir2, v3, v1);
var cross = vec3.create();
vec3.cross(cross, dir1, dir2);
vec3.normalize(cross, cross);
return cross;
}
var tessy = (function initTesselator() {
// function called for each vertex of tesselator output
function vertexCallback(data, polyVertArray) {
polyVertArray[polyVertArray.length] = data[0];
polyVertArray[polyVertArray.length] = data[1];
polyVertArray[polyVertArray.length] = data[2];
}
function begincallback(type) {
}
function errorcallback(errno) {
}
// callback for when segments intersect and must be split
function combinecallback(coords, data, weight) {
return [coords[0], coords[1], coords[2]];
}
function edgeCallback(flag) {
// don't really care about the flag, but need no-strip/no-fan behavior
}
var tessy = new libtess.GluTesselator();
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_VERTEX_DATA, vertexCallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_BEGIN, begincallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_ERROR, errorcallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_COMBINE, combinecallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_EDGE_FLAG, edgeCallback);
return tessy;
})();
export { BBox, createParser, readStream, readFile };
// Parsing of whitespace-separated numbers (gml:posList / gml:pos text) without
// splitting into strings. Gives exactly the same doubles as
// text.trim().split(/\s+/).map(parseFloat).
//
// Fast path: plain decimals with at most 15 digits and no exponent. Their
// digits form an integer below 2^53 and the divisor is a power of ten up to
// 10^15; both are exact doubles, so one IEEE division gives the correctly
// rounded value, which is what parseFloat returns. Everything else (exponents,
// longer numbers, malformed tokens) goes through parseFloat.
const POWERS_OF_TEN = [];
for (let i = 0; i <= 15; i++) {
POWERS_OF_TEN.push(10 ** i);
}
const MAX_FAST_DIGITS = 15;
// Same set as \s in JavaScript regular expressions.
function isWhitespace(c) {
if (c <= 32) {
return c === 32 || (c >= 9 && c <= 13);
}
return c === 0xA0 || c === 0x1680 || (c >= 0x2000 && c <= 0x200A) || c === 0x2028 || c === 0x2029 ||
c === 0x202F || c === 0x205F || c === 0x3000 || c === 0xFEFF;
}
// Returns the numbers in text, in order.
export function parseNumbers(text) {
const numbers = [];
const n = text.length;
let i = 0;
while (i < n) {
while (i < n && isWhitespace(text.charCodeAt(i))) {
i++;
}
if (i >= n) {
break;
}
const start = i;
let c = text.charCodeAt(i);
let negative = false;
if (c === 45 || c === 43) { // - +
negative = c === 45;
c = text.charCodeAt(++i);
}
let mantissa = 0;
let digits = 0;
let fractionDigits = 0;
while (c >= 48 && c <= 57) {
mantissa = mantissa * 10 + (c - 48);
digits++;
c = text.charCodeAt(++i);
}
if (c === 46) { // .
c = text.charCodeAt(++i);
while (c >= 48 && c <= 57) {
mantissa = mantissa * 10 + (c - 48);
digits++;
fractionDigits++;
c = text.charCodeAt(++i);
}
}
if (digits > 0 && digits <= MAX_FAST_DIGITS && (i >= n || isWhitespace(c))) {
const value = mantissa / POWERS_OF_TEN[fractionDigits];
numbers.push(negative ? -value : value);
continue;
}
// Slow path: the whole token through parseFloat.
while (i < n && !isWhitespace(text.charCodeAt(i))) {
i++;
}
numbers.push(parseFloat(text.slice(start, i)));
}
return numbers;
}
// Triangulation of planar polygons (exterior ring + interior rings, each a
// list of [x, y, z]). Triangles and quadrilaterals without holes, the bulk
// of LoD2 surfaces, are split directly; everything else goes through libtess.
import { libtess, glMatrix } from '../vendor.js';
const vec3 = glMatrix.vec3;
// Returns { vertices: [x, y, z, ...] (three per triangle), normal: vec3 }.
export function triangulate(polygon) {
const normal = calculateNormal(polygon[0]);
const vertices = polygon.length === 1 ? triangulateSimple(polygon[0], normal) : null;
return {
vertices: vertices !== null ? vertices : triangulateWithLibtess(polygon, normal),
normal: normal,
};
}
// Corners of a ring without the repeated closing point, or null if two
// consecutive corners are equal (left to libtess, which drops them).
function corners(ring) {
let n = ring.length;
if (n > 1 && samePoint(ring[0], ring[n - 1])) {
n--;
}
const points = ring.slice(0, n);
for (let i = 0; i < n; i++) {
if (samePoint(points[i], points[(i + 1) % n])) {
return null;
}
}
return points;
}
function samePoint(a, b) {
return a[0] === b[0] && a[1] === b[1] && a[2] === b[2];
}
// libtess projects the polygon to 2D by dropping the coordinate along the
// normal's largest component and merges points that then coincide. Such
// polygons are left to libtess so that the result stays the same.
function coincideInProjection(points, normal) {
const ax = Math.abs(normal[0]), ay = Math.abs(normal[1]), az = Math.abs(normal[2]);
const drop = ax > ay && ax > az ? 0 : ay > az ? 1 : 2;
const s = (drop + 1) % 3, t = (drop + 2) % 3;
for (let i = 0; i < points.length; i++) {
for (let j = i + 1; j < points.length; j++) {
if (points[i][s] === points[j][s] && points[i][t] === points[j][t]) {
return true;
}
}
}
return false;
}
// Turn at corner b (from a over b to c) relative to the polygon normal:
// positive for a left turn, negative for a right turn, near 0 if collinear.
function turn(a, b, c, normal) {
const e1x = b[0] - a[0], e1y = b[1] - a[1], e1z = b[2] - a[2];
const e2x = c[0] - b[0], e2y = c[1] - b[1], e2z = c[2] - b[2];
const cx = e1y * e2z - e1z * e2y, cy = e1z * e2x - e1x * e2z, cz = e1x * e2y - e1y * e2x;
const t = cx * normal[0] + cy * normal[1] + cz * normal[2];
const scale = Math.hypot(e1x, e1y, e1z) * Math.hypot(e2x, e2y, e2z);
return Math.abs(t) <= 1e-9 * scale ? 0 : t;
}
function pushTriangle(out, a, b, c) {
out.push(a[0], a[1], a[2], b[0], b[1], b[2], c[0], c[1], c[2]);
}
// Triangles and simple quadrilaterals, or null to use libtess.
function triangulateSimple(ring, normal) {
const p = corners(ring);
if (p === null || (p.length !== 3 && p.length !== 4) || coincideInProjection(p, normal)) {
return null;
}
const turns = p.map((b, i) => turn(p[(i + p.length - 1) % p.length], b, p[(i + 1) % p.length], normal));
if (turns.some(t => t === 0)) {
return null; // degenerate: collinear corners
}
const out = [];
if (p.length === 3) {
pushTriangle(out, p[0], p[1], p[2]);
return out;
}
// The normal points so that most turns are left turns. A convex
// quadrilateral has four of them; a concave one has exactly one right turn
// at its reflex corner, and only the diagonal from that corner lies inside.
// Two right turns mean the quadrilateral crosses itself.
const right = turns.filter(t => t < 0).length;
let r;
if (right === 0) {
r = 0;
} else if (right === 1) {
r = turns.findIndex(t => t < 0);
} else {
return null;
}
pushTriangle(out, p[r], p[(r + 1) % 4], p[(r + 2) % 4]);
pushTriangle(out, p[r], p[(r + 2) % 4], p[(r + 3) % 4]);
return out;
}
function triangulateWithLibtess(polygon, normal) {
const triangleVerts = [];
tessy.gluTessBeginPolygon(triangleVerts);
tessy.gluTessNormal(normal[0], normal[1], normal[2]);
for (const contour of polygon) {
tessy.gluTessBeginContour();
for (const point of contour) {
tessy.gluTessVertex(point, point);
}
tessy.gluTessEndContour();
}
tessy.gluTessEndPolygon();
return triangleVerts;
}
function calculateNormal(ring) {
const coords = [0, 0, 0];
for (let i = 0; i < ring.length - 1; i++) {
const current = ring[i];
const next = ring[i + 1];
coords[0] += (current[2] + next[2]) * (current[1] - next[1]);
coords[1] += (current[0] + next[0]) * (current[2] - next[2]);
coords[2] += (current[1] + next[1]) * (current[0] - next[0]);
}
if (coords[0] === 0 && coords[1] === 0 && coords[2] === 0) {
// no valid normal vector found
if (ring.length < 3) {
// no three points, return x-axis
return vec3.fromValues(1, 0, 0);
}
// vec3.create() takes no arguments; clone copies the points.
return calculateNormalWithCross(vec3.clone(ring[0]), vec3.clone(ring[1]), vec3.clone(ring[2]));
}
const v = vec3.fromValues(coords[0], coords[1], coords[2]);
vec3.normalize(v, v);
return v;
}
function calculateNormalWithCross(v1, v2, v3) {
const dir1 = vec3.create();
vec3.sub(dir1, v2, v1);
const dir2 = vec3.create();
vec3.sub(dir2, v3, v1);
const cross = vec3.create();
vec3.cross(cross, dir1, dir2);
vec3.normalize(cross, cross);
return cross;
}
// One tesselator per module instance (per thread, if the parser later runs
// in workers).
const tessy = (function initTesselator() {
// function called for each vertex of tesselator output
function vertexCallback(data, polyVertArray) {
polyVertArray[polyVertArray.length] = data[0];
polyVertArray[polyVertArray.length] = data[1];
polyVertArray[polyVertArray.length] = data[2];
}
function begincallback() {
}
function errorcallback() {
}
// callback for when segments intersect and must be split
function combinecallback(coords) {
return [coords[0], coords[1], coords[2]];
}
function edgeCallback() {
// don't really care about the flag, but need no-strip/no-fan behavior
}
const tessy = new libtess.GluTesselator();
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_VERTEX_DATA, vertexCallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_BEGIN, begincallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_ERROR, errorcallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_COMBINE, combinecallback);
tessy.gluTessCallback(libtess.gluEnum.GLU_TESS_EDGE_FLAG, edgeCallback);
return tessy;
})();
import test from 'node:test';
import assert from 'node:assert/strict';
import { parseNumbers } from '../public/src/loading/coordinates.js';
// The reference: what the parser did before (split on whitespace, parseFloat).
const reference = text => (text.trim() === '' ? [] : text.trim().split(/\s+/).map(parseFloat));
function assertSameAsParseFloat(text) {
const actual = parseNumbers(text);
const expected = reference(text);
assert.equal(actual.length, expected.length, `count for ${JSON.stringify(text)}`);
for (let i = 0; i < expected.length; i++) {
if (!Object.is(actual[i], expected[i])) {
assert.fail(`number ${i} of ${JSON.stringify(text.slice(0, 200))}: ${actual[i]} !== ${expected[i]}`);
}
}
}
test('coordinates separated by any whitespace are parsed in order', () => {
assert.deepEqual(parseNumbers(' 3513000.37\t5403456.71\n250 -0.5\r\n'), [3513000.37, 5403456.71, 250, -0.5]);
});
test('empty and whitespace-only text gives no numbers', () => {
assert.deepEqual(parseNumbers(''), []);
assert.deepEqual(parseNumbers(' \n\t '), []);
});
test('random decimal numbers give exactly the same doubles as parseFloat', () => {
// Seeded generator, so a failure can be reproduced.
let seed = 12345;
const random = () => (seed = (seed * 1103515245 + 12345) % 2147483648) / 2147483648;
const numbers = [];
for (let i = 0; i < 20000; i++) {
const magnitude = 10 ** Math.floor(random() * 9 - 2);
const value = (random() - 0.3) * magnitude;
numbers.push(value.toFixed(Math.floor(random() * 9)));
}
assertSameAsParseFloat(numbers.join(' '));
});
test('special notations give the same doubles as parseFloat', () => {
assertSameAsParseFloat(
'0 -0 +0 5. .5 -.5 +7.25 007.10 1e3 1.5E-2 -2.5e+4 ' +
'12345678901234567890 0.12345678901234567890 123456789012345.6 9007199254740993 ' +
'3513000.3700000001 1e400 -1e-400');
});
test('numbers with more digits than a double holds exactly give the same doubles as parseFloat', () => {
// 17 digits: accumulating the digits would round before the division.
assertSameAsParseFloat('259658909219030.06 70642.808107788870 904196606.41160686 6.7285402536930302');
});
test('malformed tokens behave like parseFloat', () => {
assertSameAsParseFloat('abc 1.5abc - . -. 1.2.3 1e 1e+ NaN Infinity -Infinity 0x10');
});
test('unicode whitespace separates numbers like split(/\\s+/)', () => {
assertSameAsParseFloat('1 2 3  45
 6');
});
import test from 'node:test';
import assert from 'node:assert/strict';
import { triangulate } from '../public/src/loading/triangulate.js';
// Sum of the areas of the triangles in a flat [x, y, z, ...] vertex list.
function area(vertices) {
let sum = 0;
for (let i = 0; i < vertices.length; i += 9) {
const ax = vertices[i + 3] - vertices[i], ay = vertices[i + 4] - vertices[i + 1], az = vertices[i + 5] - vertices[i + 2];
const bx = vertices[i + 6] - vertices[i], by = vertices[i + 7] - vertices[i + 1], bz = vertices[i + 8] - vertices[i + 2];
sum += Math.hypot(ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx) / 2;
}
return sum;
}
// Closed ring in the plane z = 0 from [x, y] pairs.
const ring = (...points) => [...points, points[0]].map(([x, y]) => [x, y, 0]);
const near = (actual, expected) => assert.ok(Math.abs(actual - expected) < 1e-9, `${actual} != ${expected}`);
test('a triangle gives one triangle with its corners', () => {
const { vertices } = triangulate([ring([0, 0], [2, 0], [0, 1])]);
assert.deepEqual(vertices, [0, 0, 0, 2, 0, 0, 0, 1, 0]);
});
test('a convex quadrilateral gives two triangles covering its area', () => {
const { vertices } = triangulate([ring([0, 0], [3, 0], [3, 2], [0, 2])]);
assert.equal(vertices.length, 18);
near(area(vertices), 6);
});
test('a vertical wall is triangulated in its own plane', () => {
const wall = [[0, 0, 0], [4, 0, 0], [4, 0, 3], [0, 0, 3], [0, 0, 0]];
const { vertices, normal } = triangulate([wall]);
assert.equal(vertices.length, 18);
near(area(vertices), 12);
near(Math.abs(normal[1]), 1);
});
test('a concave quadrilateral is split at its reflex corner so no triangle lies outside', () => {
// (2, 1) is the reflex corner; splitting along (0,0)-(4,0) would cover the notch.
const { vertices } = triangulate([ring([0, 0], [2, 1], [4, 0], [2, 4])]);
assert.equal(vertices.length, 18);
near(area(vertices), 6);
});
test('the reflex corner is found at any position in the ring', () => {
const points = [[0, 0], [2, 1], [4, 0], [2, 4]];
for (let shift = 0; shift < 4; shift++) {
const rotated = points.map((_, i) => points[(i + shift) % 4]);
near(area(triangulate([ring(...rotated)]).vertices), 6);
}
});
test('a ring without the repeated closing point is handled the same way', () => {
const open = [[0, 0, 0], [3, 0, 0], [3, 2, 0], [0, 2, 0]];
near(area(triangulate([open]).vertices), 6);
});
test('a self-intersecting quadrilateral is left to the general triangulation', () => {
// Bow-tie: two triangles of area 0.25 each.
const { vertices } = triangulate([ring([0, 0], [1, 1], [1, 0], [0, 1])]);
near(area(vertices), 0.5);
});
test('a quadrilateral with a hole is left to the general triangulation', () => {
const outer = ring([0, 0], [4, 0], [4, 4], [0, 4]);
const hole = ring([1, 1], [1, 3], [3, 3], [3, 1]);
near(area(triangulate([outer, hole]).vertices), 12);
});
test('repeated and collinear points are left to the general triangulation', () => {
near(area(triangulate([ring([0, 0], [2, 0], [2, 0], [0, 2])]).vertices), 2);
near(area(triangulate([ring([0, 0], [1, 0], [2, 0], [0, 2])]).vertices), 2);
});
test('corners that coincide when projected along the normal are merged like libtess does', () => {
// Real roof from Oberstadt4736.xml: two corners differ only by 1 mm in z,
// the main axis of the normal. libtess projects them onto the same point
// and gives one triangle; the fast path must not add a 1 mm sliver.
const roof = [[0, 0, 0], [0.645, -4.137, -0.811], [6.813, -2.88, -0.812], [6.813, -2.88, -0.811], [0, 0, 0]];
assert.equal(triangulate([roof]).vertices.length, 9);
});
test('larger polygons are triangulated as before', () => {
// L-shape with six corners, area 3.
const { vertices } = triangulate([ring([0, 0], [2, 0], [2, 1], [1, 1], [1, 2], [0, 2])]);
near(area(vertices), 3);
});
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment