utils/TerrainWorker.worker.jsjavascript
"use strict";
//
//Terrain worker
//
import { getMedian, getMeshHeights } from "../terrain/terrainUtils";
import { doubleToTwoFloatsV3 } from "../math/doubleToTwoFloats";
function binarySearchFast(arr, x) {
let start = 0,
end = arr.length - 1;
while (start <= end) {
let k = Math.floor((start + end) * 0.5);
if (Math.abs(arr[k] - x) < 1e-3) return k;
else if (arr[k] < x) start = k + 1;
else end = k - 1;
}
return -1;
}
function checkNoDataValue(noDataValues, value) {
return binarySearchFast(noDataValues, value) !== -1;
}
var Vec3 = function (x, y, z) {
this.x = x;
this.y = y;
this.z = z;
};
var blerp = function (x, y, fQ11, fQ21, fQ12, fQ22) {
return fQ11 * (1.0 - x) * (1.0 - y) + fQ21 * x * (1.0 - y) + fQ12 * (1.0 - x) * y + fQ22 * x * y;
};
var _tempVec = new Vec3(0.0, 0.0, 0.0);
var _rtcTempVec = new Vec3(0.0, 0.0, 0.0);
var _tempHigh = new Vec3(0.0, 0.0, 0.0),
_tempLow = new Vec3(0.0, 0.0, 0.0);
var _meshHeights = new Float32Array(0);
// Squared vertex radii, their median is the LOD radius
var _r2 = new Float64Array(0);
// Geocentric file grid vertices, the quad normals are made from them
var _world = new Float64Array(0);
self.onmessage = function (e) {
var elevations = e.data.elevations,
this_plainVertices = e.data.this_plainVertices,
this_plainNormals = e.data.this_plainNormals,
this_normalMapVertices = e.data.this_normalMapVertices,
this_normalMapNormals = e.data.this_normalMapNormals,
heightFactor = e.data.heightFactor,
gridSize = e.data.gridSize,
noDataValues = e.data.noDataValues,
id = e.data.id;
let rtc_x = e.data.relativeCenter[0],
rtc_y = e.data.relativeCenter[1],
rtc_z = e.data.relativeCenter[2];
var xmin = 549755748352.0,
xmax = -549755748352.0,
ymin = 549755748352.0,
ymax = -549755748352.0,
zmin = 549755748352.0,
zmax = -549755748352.0;
// Squared geocentric distance of the same vertices. The axis aligned bounds above cannot tell
// relief from segment size, because the box of a tilted patch is far larger than the patch, so
// the radius range is accumulated here instead. Kept squared to avoid a root per vertex.
var r2min = Infinity,
r2max = -Infinity;
const fileGridSize = Math.sqrt(elevations.length) - 1;
const fileGridSize_one = fileGridSize + 1;
const fileGridSize_one_x2 = fileGridSize_one * fileGridSize_one;
const tgs = gridSize;
const dg = fileGridSize / tgs;
const gs = tgs + 1;
const hf = heightFactor;
var vInd = 0;
var gsgs3 = gs * gs * 3;
if (_r2.length !== gs * gs) {
_r2 = new Float64Array(gs * gs);
}
var terrainVertices = new Float64Array(gsgs3),
terrainVerticesHigh = new Float32Array(gsgs3),
terrainVerticesLow = new Float32Array(gsgs3);
var normalMapNormals, normalMapVertices;
var nv = this_normalMapVertices,
nn = this_normalMapNormals;
if (fileGridSize >= tgs) {
normalMapNormals = new Float32Array(fileGridSize_one_x2 * 3);
normalMapVertices = new Float64Array(fileGridSize_one_x2 * 3);
// Mesh heights keep the elevations peaks and pits, moved to the nearest mesh vertices
let meshHeights = null;
if (dg > 1) {
if (_meshHeights.length !== gs * gs) {
_meshHeights = new Float32Array(gs * gs);
}
meshHeights = getMeshHeights(
elevations,
fileGridSize,
tgs,
noDataValues.length ? noDataValues[0] : NaN,
_meshHeights
);
}
if (_world.length < fileGridSize_one_x2 * 3) {
_world = new Float64Array(fileGridSize_one_x2 * 3);
}
//
// Vertices, each one once
//
for (let k = 0; k < fileGridSize_one_x2; k++) {
var j = k % fileGridSize_one,
i = ~~(k / fileGridSize_one);
var vInd0 = k * 3;
var currElv = elevations[k];
if (checkNoDataValue(noDataValues, currElv)) {
currElv = 0.0;
}
var h0 = hf * currElv;
var v0x = rtc_x + nv[vInd0] + h0 * nn[vInd0],
v0y = rtc_y + nv[vInd0 + 1] + h0 * nn[vInd0 + 1],
v0z = rtc_z + nv[vInd0 + 2] + h0 * nn[vInd0 + 2];
_world[vInd0] = v0x;
_world[vInd0 + 1] = v0y;
_world[vInd0 + 2] = v0z;
normalMapVertices[vInd0] = v0x - rtc_x;
normalMapVertices[vInd0 + 1] = v0y - rtc_y;
normalMapVertices[vInd0 + 2] = v0z - rtc_z;
//
// The vertex goes into screen buffer
if (i % dg === 0 && j % dg === 0) {
let mx = v0x,
my = v0y,
mz = v0z;
if (meshHeights) {
let meshElv = meshHeights[(i / dg) * gs + j / dg];
// Only the vertices a peak or pit is moved to
if (meshElv !== elevations[k] && !checkNoDataValue(noDataValues, meshElv)) {
let mh = hf * meshElv;
mx = rtc_x + nv[vInd0] + mh * nn[vInd0];
my = rtc_y + nv[vInd0 + 1] + mh * nn[vInd0 + 1];
mz = rtc_z + nv[vInd0 + 2] + mh * nn[vInd0 + 2];
}
}
if (mx < xmin) xmin = mx;
if (mx > xmax) xmax = mx;
if (my < ymin) ymin = my;
if (my > ymax) ymax = my;
if (mz < zmin) zmin = mz;
if (mz > zmax) zmax = mz;
let r2_v0 = mx * mx + my * my + mz * mz;
if (r2_v0 < r2min) r2min = r2_v0;
if (r2_v0 > r2max) r2max = r2_v0;
_r2[vInd / 3] = r2_v0;
_rtcTempVec.x = mx - rtc_x;
_rtcTempVec.y = my - rtc_y;
_rtcTempVec.z = mz - rtc_z;
doubleToTwoFloatsV3(_rtcTempVec, _tempHigh, _tempLow);
terrainVerticesHigh[vInd] = _tempHigh.x;
terrainVerticesLow[vInd] = _tempLow.x;
terrainVertices[vInd++] = _rtcTempVec.x;
terrainVerticesHigh[vInd] = _tempHigh.y;
terrainVerticesLow[vInd] = _tempLow.y;
terrainVertices[vInd++] = _rtcTempVec.y;
terrainVerticesHigh[vInd] = _tempHigh.z;
terrainVerticesLow[vInd] = _tempLow.z;
terrainVertices[vInd++] = _rtcTempVec.z;
}
}
//
// Normals of the quads V0 V1 / V2 V3, added to their vertices in the quads order
//
for (let i = 0; i < fileGridSize; i++) {
for (let j = 0; j < fileGridSize; j++) {
let vInd0 = (i * fileGridSize_one + j) * 3,
vInd1 = vInd0 + 3,
vInd2 = vInd0 + fileGridSize_one * 3,
vInd3 = vInd2 + 3;
let x0 = _world[vInd0],
y0 = _world[vInd0 + 1],
z0 = _world[vInd0 + 2];
let e10x = _world[vInd1] - x0,
e10y = _world[vInd1 + 1] - y0,
e10z = _world[vInd1 + 2] - z0;
let e20x = _world[vInd2] - x0,
e20y = _world[vInd2 + 1] - y0,
e20z = _world[vInd2 + 2] - z0;
let e30x = _world[vInd3] - x0,
e30y = _world[vInd3 + 1] - y0,
e30z = _world[vInd3 + 2] - z0;
// sw = e20 x e30, ne = e30 x e10, n0 = ne + sw, normalized
let swx = e20y * e30z - e20z * e30y,
swy = e20z * e30x - e20x * e30z,
swz = e20x * e30y - e20y * e30x;
let l = 1.0 / Math.sqrt(swx * swx + swy * swy + swz * swz);
swx *= l;
swy *= l;
swz *= l;
let nex = e30y * e10z - e30z * e10y,
ney = e30z * e10x - e30x * e10z,
nez = e30x * e10y - e30y * e10x;
l = 1.0 / Math.sqrt(nex * nex + ney * ney + nez * nez);
nex *= l;
ney *= l;
nez *= l;
let n0x = nex + swx,
n0y = ney + swy,
n0z = nez + swz;
l = 1.0 / Math.sqrt(n0x * n0x + n0y * n0y + n0z * n0z);
n0x *= l;
n0y *= l;
n0z *= l;
normalMapNormals[vInd0] += n0x;
normalMapNormals[vInd0 + 1] += n0y;
normalMapNormals[vInd0 + 2] += n0z;
normalMapNormals[vInd1] += nex;
normalMapNormals[vInd1 + 1] += ney;
normalMapNormals[vInd1 + 2] += nez;
normalMapNormals[vInd2] += swx;
normalMapNormals[vInd2 + 1] += swy;
normalMapNormals[vInd2 + 2] += swz;
normalMapNormals[vInd3] += n0x;
normalMapNormals[vInd3 + 1] += n0y;
normalMapNormals[vInd3 + 2] += n0z;
}
}
} else {
normalMapNormals = new Float32Array(gsgs3);
normalMapVertices = new Float64Array(gsgs3);
var oneSize = tgs / fileGridSize;
var gsgs = gsgs3 / 3;
var fgsOne = fileGridSize + 1;
for (let i = 0; i < gsgs; i++) {
let ii = Math.floor(i / gs),
ij = i % gs;
let qii = ii % oneSize,
qij = ij % oneSize;
let hlt_ind = Math.floor(ii / oneSize) * fgsOne + Math.floor(ij / oneSize);
if (ij === tgs) {
hlt_ind -= 1;
qij = oneSize;
}
if (ii === tgs) {
hlt_ind -= fgsOne;
qii = oneSize;
}
let hrt_ind = hlt_ind + 1,
hlb_ind = hlt_ind + fgsOne,
hrb_ind = hlb_ind + 1;
let h_lt = elevations[hlt_ind],
h_rt = elevations[hrt_ind],
h_lb = elevations[hlb_ind],
h_rb = elevations[hrb_ind];
if (checkNoDataValue(noDataValues, h_lt)) {
h_lt = 0.0;
}
if (checkNoDataValue(noDataValues, h_rt)) {
h_rt = 0.0;
}
if (checkNoDataValue(noDataValues, h_lb)) {
h_lb = 0.0;
}
if (checkNoDataValue(noDataValues, h_rb)) {
h_rb = 0.0;
}
let hi = blerp(qij / oneSize, qii / oneSize, h_lt, h_rt, h_lb, h_rb);
let i3 = i * 3;
_rtcTempVec.x = this_plainVertices[i3] + hi * this_plainNormals[i3];
_rtcTempVec.y = this_plainVertices[i3 + 1] + hi * this_plainNormals[i3 + 1];
_rtcTempVec.z = this_plainVertices[i3 + 2] + hi * this_plainNormals[i3 + 2];
doubleToTwoFloatsV3(_rtcTempVec, _tempHigh, _tempLow);
terrainVertices[i3] = _rtcTempVec.x;
terrainVertices[i3 + 1] = _rtcTempVec.y;
terrainVertices[i3 + 2] = _rtcTempVec.z;
terrainVerticesHigh[i3] = _tempHigh.x;
terrainVerticesHigh[i3 + 1] = _tempHigh.y;
terrainVerticesHigh[i3 + 2] = _tempHigh.z;
terrainVerticesLow[i3] = _tempLow.x;
terrainVerticesLow[i3 + 1] = _tempLow.y;
terrainVerticesLow[i3 + 2] = _tempLow.z;
_tempVec.x = _rtcTempVec.x + rtc_x;
_tempVec.y = _rtcTempVec.y + rtc_y;
_tempVec.z = _rtcTempVec.z + rtc_z;
if (_tempVec.x < xmin) xmin = _tempVec.x;
if (_tempVec.x > xmax) xmax = _tempVec.x;
if (_tempVec.y < ymin) ymin = _tempVec.y;
if (_tempVec.y > ymax) ymax = _tempVec.y;
if (_tempVec.z < zmin) zmin = _tempVec.z;
if (_tempVec.z > zmax) zmax = _tempVec.z;
let r2_tmp = _tempVec.x * _tempVec.x + _tempVec.y * _tempVec.y + _tempVec.z * _tempVec.z;
if (r2_tmp < r2min) r2min = r2_tmp;
if (r2_tmp > r2max) r2max = r2_tmp;
_r2[i] = r2_tmp;
}
normalMapVertices.set(terrainVertices);
for (let k = 0; k < gsgs; k++) {
let j = k % gs,
i = ~~(k / gs);
if (i !== tgs && j !== tgs) {
let v0ind = k * 3,
v1ind = v0ind + 3,
v2ind = v0ind + gs * 3,
v3ind = v2ind + 3;
let x0 = terrainVertices[v0ind] + rtc_x,
y0 = terrainVertices[v0ind + 1] + rtc_y,
z0 = terrainVertices[v0ind + 2] + rtc_z;
// Edges are normalized here before the cross products
let e10x = terrainVertices[v1ind] + rtc_x - x0,
e10y = terrainVertices[v1ind + 1] + rtc_y - y0,
e10z = terrainVertices[v1ind + 2] + rtc_z - z0;
let l = 1.0 / Math.sqrt(e10x * e10x + e10y * e10y + e10z * e10z);
e10x *= l;
e10y *= l;
e10z *= l;
let e20x = terrainVertices[v2ind] + rtc_x - x0,
e20y = terrainVertices[v2ind + 1] + rtc_y - y0,
e20z = terrainVertices[v2ind + 2] + rtc_z - z0;
l = 1.0 / Math.sqrt(e20x * e20x + e20y * e20y + e20z * e20z);
e20x *= l;
e20y *= l;
e20z *= l;
let e30x = terrainVertices[v3ind] + rtc_x - x0,
e30y = terrainVertices[v3ind + 1] + rtc_y - y0,
e30z = terrainVertices[v3ind + 2] + rtc_z - z0;
l = 1.0 / Math.sqrt(e30x * e30x + e30y * e30y + e30z * e30z);
e30x *= l;
e30y *= l;
e30z *= l;
// sw = e20 x e30, ne = e30 x e10, n0 = ne + sw, normalized
let swx = e20y * e30z - e20z * e30y,
swy = e20z * e30x - e20x * e30z,
swz = e20x * e30y - e20y * e30x;
l = 1.0 / Math.sqrt(swx * swx + swy * swy + swz * swz);
swx *= l;
swy *= l;
swz *= l;
let nex = e30y * e10z - e30z * e10y,
ney = e30z * e10x - e30x * e10z,
nez = e30x * e10y - e30y * e10x;
l = 1.0 / Math.sqrt(nex * nex + ney * ney + nez * nez);
nex *= l;
ney *= l;
nez *= l;
let n0x = nex + swx,
n0y = ney + swy,
n0z = nez + swz;
l = 1.0 / Math.sqrt(n0x * n0x + n0y * n0y + n0z * n0z);
n0x *= l;
n0y *= l;
n0z *= l;
normalMapNormals[v0ind] += n0x;
normalMapNormals[v0ind + 1] += n0y;
normalMapNormals[v0ind + 2] += n0z;
normalMapNormals[v1ind] += nex;
normalMapNormals[v1ind + 1] += ney;
normalMapNormals[v1ind + 2] += nez;
normalMapNormals[v2ind] += swx;
normalMapNormals[v2ind + 1] += swy;
normalMapNormals[v2ind + 2] += swz;
normalMapNormals[v3ind] += n0x;
normalMapNormals[v3ind + 1] += n0y;
normalMapNormals[v3ind + 2] += n0z;
}
}
}
self.postMessage(
{
id: id,
normalMapNormals: normalMapNormals,
normalMapVertices: normalMapVertices,
terrainVertices: terrainVertices,
terrainVerticesHigh: terrainVerticesHigh,
terrainVerticesLow: terrainVerticesLow,
bounds: [xmin, ymin, zmin, xmax, ymax, zmax],
radiusRange: r2max >= r2min ? [Math.sqrt(r2min), Math.sqrt(r2max)] : null,
lodRadius: Math.sqrt(getMedian(_r2, gs * gs)),
relativeCenter: e.data.relativeCenter
},
[
normalMapNormals.buffer,
normalMapVertices.buffer,
terrainVertices.buffer,
terrainVerticesHigh.buffer,
terrainVerticesLow.buffer
]
);
};