Resources Numerical accuracy Calculation policy
Finds which points of a drone or scanner cloud are ground, so a surface can be built from them.
The progressive morphological filter of Zhang et al. (2003), which the file's header cites, run on a grid that holds the lowest point of each cell in a Float64Array. The result is a selection: ground points are returned with their original coordinates. Which points are selected depends on the cell size and on the slope and height thresholds.
5 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests bin/test-civilcraft-ground.mjs
The heights the filter uses when the caller gives none.
h₀ = 0.3 × cell; h_max = 6 × cell; low noise = slope × cell × 3 + 3 × h₀
Plain 64-bit products of the cell size. Each can be replaced by a value from the caller.
const h0 = opts.initialHeight ?? cell * 0.3, hMax = opts.maxHeight ?? cell * 6, lowNoise = opts.lowNoise ?? slope * cell * 3 + 3 * h0;Sorts the points into grid cells and keeps the lowest of each.
column = ⌊(x − minₓ) ÷ cell⌋; row = ⌊(y − min_y) ÷ cell⌋
Cell indexes are floored from 64-bit quotients measured from the data's own minimum, so the quotients stay small. Each cell keeps the lowest elevation and the index of that point.
const nx = Math.floor((maxX - minX) / cell) + 1, ny = Math.floor((maxY - minY) / cell) + 1;
if (nx * ny > maxCells) throw Object.assign(new Error('classifyGround: too many cells'), { code: 'cells', cells: nx * ny, max: maxCells }); // i18n-ignore: the UI words it from the code
const N = nx * ny, low = new Float64Array(N).fill(Infinity), pick = new Int32Array(N).fill(-1);
points.forEach((p, i) => {
const c = Math.floor((p.y - minY) / cell) * nx + Math.floor((p.x - minX) / cell);
if (p.z < low[c]) { low[c] = p.z; pick[c] = i; }
});
const has = new Uint8Array(N); // 1 = a real lowest point; 0 = empty
for (let c = 0; c < N; c++) if (pick[c] >= 0) { has[c] = 1; counts.cells++; }Drops single points that lie far below their neighbors.
stray when median(neighbors) − z > low-noise limit
A cell with at least 4 of its 8 neighbors filled is dropped when it lies more than the limit under their median. For an even count the median is the mean of the two middle values.
// stray low points: well under the middle of the 8 neighbours that have points
const stray = new Uint8Array(N), nb = [];
for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) {
const c = j * nx + i;
if (!has[c]) continue;
nb.length = 0;
for (let dj = -1; dj <= 1; dj++) for (let di = -1; di <= 1; di++) {
if (!di && !dj) continue;
const a = i + di, b = j + dj;
if (a >= 0 && a < nx && b >= 0 && b < ny && has[b * nx + a]) nb.push(low[b * nx + a]);
}
if (nb.length >= 4) { nb.sort((p, q) => p - q); const med = (nb[(nb.length - 1) >> 1] + nb[nb.length >> 1]) / 2; if (med - low[c] > lowNoise) stray[c] = 1; }
}Marks cells that stand above the opened surface as objects, round by round.
wₖ = 2ᵏ + 1; dhₖ = min(h_max, slope × (wₖ − wₖ₋₁) × cell + h₀); object when z − opened > dhₖ
Each round opens the previous round's surface and marks cells that rise more than the threshold above it. The loop ends when the window covers the largest object or is larger than the grid.
const object = new Uint8Array(N), tmp = new Float64Array(N);
let cur = z, prevW = 1;
const windows = [];
for (let k = 1; ; k++) {
const w = 2 ** k + 1;
windows.push(w);
const r = (w - 1) >> 1;
const opened = open(cur, nx, ny, r, tmp);
const dh = Math.min(hMax, slope * (w - prevW) * cell + h0);
for (let c = 0; c < N; c++) if (has[c] && !stray[c] && cur[c] - opened[c] > dh) object[c] = 1;
cur = opened; prevW = w;
if (w * cell >= maxObject || w > Math.max(nx, ny)) break;
}The minimum, then the maximum, over a square window around each cell.
Each of the two passes is done along x and then along y with the edges clamped (the sweep function below it). Only comparisons of 64-bit elevations, so no rounding is introduced.
/** Opening (erosion then dilation) over a (2r+1)-cell square, edges clamped. Returns a new array (valid until the next call). */
function open(src, nx, ny, r, tmp) {
const a = new Float64Array(src.length);
sweep(src, a, nx, ny, r, Math.min, Infinity, tmp);
const b = new Float64Array(src.length);
sweep(a, b, nx, ny, r, Math.max, -Infinity, tmp);
return b;
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.