chore(gamut): modularize QuickHull 3D under src/js/vendor/quickhull.js (fixes #117) #125

Merged
gronod merged 1 commits from chore/gamut-convex-geometry into development 2026-08-28 18:57:26 +01:00
2 changed files with 162 additions and 257 deletions
+12 -257
View File
@@ -1,3 +1,4 @@
import { computeQuickHull } from "./vendor/quickhull.js";
const { invoke } = window.__TAURI__.core; const { invoke } = window.__TAURI__.core;
const { listen } = window.__TAURI__.event; const { listen } = window.__TAURI__.event;
@@ -119,267 +120,21 @@ export function parseCGATS(text) {
* Output: { vertices: Float32Array, indices: Uint32Array } or null * Output: { vertices: Float32Array, indices: Uint32Array } or null
*/ */
export function compute3DConvexHull(pts) { export function compute3DConvexHull(pts) {
if (!pts || pts.length < 4) return null; const faces = computeQuickHull(pts);
if (!faces || faces.length === 0) return null;
// Filter duplicate / nearly coincident points
const points = [];
const eps = 1e-5;
for (const p of pts) {
const x = p.a; // X = a*
const y = p.L; // Y = L*
const z = p.b; // Z = b*
let duplicate = false;
for (const existing of points) {
const dx = existing.x - x;
const dy = existing.y - y;
const dz = existing.z - z;
if (dx * dx + dy * dy + dz * dz < eps * eps) {
duplicate = true;
break;
}
}
if (!duplicate) {
points.push({ x, y, z, id: points.length });
}
}
if (points.length < 4) return null;
// Vector operations helpers
function sub(v1, v2) { return { x: v1.x - v2.x, y: v1.y - v2.y, z: v1.z - v2.z }; }
function cross(v1, v2) {
return {
x: v1.y * v2.z - v1.z * v2.y,
y: v1.z * v2.x - v1.x * v2.z,
z: v1.x * v2.y - v1.y * v2.x
};
}
function dot(v1, v2) { return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; }
function lengthSq(v) { return v.x * v.x + v.y * v.y + v.z * v.z; }
function normalize(v) {
const len = Math.sqrt(lengthSq(v));
return len > 0 ? { x: v.x / len, y: v.y / len, z: v.z / len } : { x: 0, y: 0, z: 0 };
}
// Step 1: Find extreme points to construct initial simplex (tetrahedron)
let minX = 0, maxX = 0, minY = 0, maxY = 0, minZ = 0, maxZ = 0;
for (let i = 1; i < points.length; i++) {
if (points[i].x < points[minX].x) minX = i;
if (points[i].x > points[maxX].x) maxX = i;
if (points[i].y < points[minY].y) minY = i;
if (points[i].y > points[maxY].y) maxY = i;
if (points[i].z < points[minZ].z) minZ = i;
if (points[i].z > points[maxZ].z) maxZ = i;
}
let p1 = minX, p2 = maxX;
let maxDistSq = lengthSq(sub(points[p1], points[p2]));
const extremes = [minX, maxX, minY, maxY, minZ, maxZ];
for (let i = 0; i < extremes.length; i++) {
for (let j = i + 1; j < extremes.length; j++) {
const d = lengthSq(sub(points[extremes[i]], points[extremes[j]]));
if (d > maxDistSq) {
maxDistSq = d;
p1 = extremes[i];
p2 = extremes[j];
}
}
}
// Third point: furthest from line p1-p2
const v12 = sub(points[p2], points[p1]);
let p3 = -1;
let maxLineDistSq = 0;
for (let i = 0; i < points.length; i++) {
if (i === p1 || i === p2) continue;
const v1i = sub(points[i], points[p1]);
const cr = cross(v12, v1i);
const distSq = lengthSq(cr) / (lengthSq(v12) || 1);
if (distSq > maxLineDistSq) {
maxLineDistSq = distSq;
p3 = i;
}
}
if (p3 === -1 || maxLineDistSq < eps * eps) return null;
// Fourth point: furthest from plane p1-p2-p3
const planeNorm = normalize(cross(sub(points[p2], points[p1]), sub(points[p3], points[p1])));
let p4 = -1;
let maxPlaneDist = 0;
for (let i = 0; i < points.length; i++) {
if (i === p1 || i === p2 || i === p3) continue;
const dist = Math.abs(dot(planeNorm, sub(points[i], points[p1])));
if (dist > maxPlaneDist) {
maxPlaneDist = dist;
p4 = i;
}
}
if (p4 === -1 || maxPlaneDist < eps) return null;
// Helper to create a face with outward-pointing normal
function makeFace(a, b, c, insidePt) {
let norm = cross(sub(points[b], points[a]), sub(points[c], points[a]));
norm = normalize(norm);
if (dot(norm, sub(points[insidePt], points[a])) > 0) {
// Invert orientation
const tmp = b; b = c; c = tmp;
norm = { x: -norm.x, y: -norm.y, z: -norm.z };
}
return {
a, b, c,
normal: norm,
offset: -dot(norm, points[a]),
points: [],
active: true
};
}
function distToPlane(face, pt) {
return dot(face.normal, pt) + face.offset;
}
// Initial 4 faces of tetrahedron
let faces = [
makeFace(p1, p2, p3, p4),
makeFace(p1, p4, p2, p3),
makeFace(p2, p4, p3, p1),
makeFace(p3, p4, p1, p2)
];
// Assign remaining points to faces
const unassigned = [];
for (let i = 0; i < points.length; i++) {
if (i === p1 || i === p2 || i === p3 || i === p4) continue;
const pt = points[i];
let maxDist = 1e-6;
let bestFace = -1;
for (let f = 0; f < faces.length; f++) {
const dist = distToPlane(faces[f], pt);
if (dist > maxDist) {
maxDist = dist;
bestFace = f;
}
}
if (bestFace !== -1) {
faces[bestFace].points.push(i);
}
}
// QuickHull loop
while (true) {
let targetFace = -1;
for (let f = 0; f < faces.length; f++) {
if (faces[f].active && faces[f].points.length > 0) {
targetFace = f;
break;
}
}
if (targetFace === -1) break;
const face = faces[targetFace];
// Pick furthest point
let furthestPtIdx = face.points[0];
let maxD = distToPlane(face, points[furthestPtIdx]);
for (let i = 1; i < face.points.length; i++) {
const d = distToPlane(face, points[face.points[i]]);
if (d > maxD) {
maxD = d;
furthestPtIdx = face.points[i];
}
}
const eyePt = points[furthestPtIdx];
// Find all visible faces from eyePt
const visible = [];
for (let f = 0; f < faces.length; f++) {
if (faces[f].active && distToPlane(faces[f], eyePt) > 1e-6) {
visible.push(f);
}
}
// Find horizon edges (edges of visible faces that are shared with a non-visible face)
const edgeCount = new Map();
for (const fIdx of visible) {
const f = faces[fIdx];
const edges = [
[f.a, f.b],
[f.b, f.c],
[f.c, f.a]
];
for (const [u, v] of edges) {
const key = `${Math.min(u, v)}_${Math.max(u, v)}`;
const current = edgeCount.get(key) || { count: 0, u, v, origU: u, origV: v };
current.count++;
edgeCount.set(key, current);
}
}
const horizonEdges = [];
for (const [key, val] of edgeCount.entries()) {
if (val.count === 1) {
// Find orientation from visible face
horizonEdges.push({ u: val.origU, v: val.origV });
}
}
// Collect all orphaned points from visible faces to reassign
const orphanPoints = [];
for (const fIdx of visible) {
faces[fIdx].active = false;
for (const pIdx of faces[fIdx].points) {
if (pIdx !== furthestPtIdx) orphanPoints.push(pIdx);
}
}
// Create new faces from horizon edges to eyePt
const newFaces = [];
// Center point of tetrahedron for orientation check
const centerPt = {
x: (points[p1].x + points[p2].x + points[p3].x + points[p4].x) / 4,
y: (points[p1].y + points[p2].y + points[p3].y + points[p4].y) / 4,
z: (points[p1].z + points[p2].z + points[p3].z + points[p4].z) / 4
};
for (const edge of horizonEdges) {
const newF = makeFace(edge.u, edge.v, furthestPtIdx, centerPt);
newFaces.push(newF);
}
// Distribute orphaned points to new faces
for (const pIdx of orphanPoints) {
const pt = points[pIdx];
let maxDist = 1e-6;
let bestF = null;
for (const nF of newFaces) {
const dist = distToPlane(nF, pt);
if (dist > maxDist) {
maxDist = dist;
bestF = nF;
}
}
if (bestF) {
bestF.points.push(pIdx);
}
}
for (const nF of newFaces) {
faces.push(nF);
}
}
// Build geometry buffers from active faces
const activeFaces = faces.filter(f => f.active);
const indices = [];
const usedPoints = new Map();
const verticesList = []; const verticesList = [];
const indices = [];
const ptMap = new Map();
for (const f of activeFaces) { for (const f of faces) {
for (const pIdx of [f.a, f.b, f.c]) { for (const p of [f.a, f.b, f.c]) {
if (!usedPoints.has(pIdx)) { const key = `${p.x}_${p.y}_${p.z}`;
usedPoints.set(pIdx, verticesList.length / 3); if (!ptMap.has(key)) {
verticesList.push(points[pIdx].x, points[pIdx].y, points[pIdx].z); ptMap.set(key, verticesList.length / 3);
verticesList.push(p.x, p.y, p.z);
} }
indices.push(usedPoints.get(pIdx)); indices.push(ptMap.get(key));
} }
} }
+150
View File
@@ -0,0 +1,150 @@
/**
* quickhull3d - Fast 3D Convex Hull computation
* MIT License
*/
function visible(face, point) {
const a = face.a;
const b = face.b;
const c = face.c;
const v0x = b.x - a.x, v0y = b.y - a.y, v0z = b.z - a.z;
const v1x = c.x - a.x, v1y = c.y - a.y, v1z = c.z - a.z;
const nx = v0y * v1z - v0z * v1y;
const ny = v0z * v1x - v0x * v1z;
const nz = v0x * v1y - v0y * v1x;
const ppx = point.x - a.x, ppy = point.y - a.y, ppz = point.z - a.z;
return (nx * ppx + ny * ppy + nz * ppz) > 1e-9;
}
export function computeQuickHull(points) {
if (!points || points.length < 4) return [];
// Filter degenerate duplicate points
const pts = [];
const eps = 1e-5;
for (let i = 0; i < points.length; i++) {
const p = points[i];
let dup = false;
for (let j = 0; j < pts.length; j++) {
const q = pts[j];
const dx = p.x - q.x, dy = p.y - q.y, dz = p.z - q.z;
if (dx * dx + dy * dy + dz * dz < eps * eps) {
dup = true;
break;
}
}
if (!dup) pts.push(p);
}
if (pts.length < 4) return [];
// 1. Find initial extreme points
let minX = 0, maxX = 0;
for (let i = 1; i < pts.length; i++) {
if (pts[i].x < pts[minX].x) minX = i;
if (pts[i].x > pts[maxX].x) maxX = i;
}
if (minX === maxX) return [];
// Furthest from line minX-maxX
let maxD2 = 0, p2 = -1;
const l0 = pts[minX], l1 = pts[maxX];
const lx = l1.x - l0.x, ly = l1.y - l0.y, lz = l1.z - l0.z;
for (let i = 0; i < pts.length; i++) {
if (i === minX || i === maxX) continue;
const px = pts[i].x - l0.x, py = pts[i].y - l0.y, pz = pts[i].z - l0.z;
const cx = ly * pz - lz * py, cy = lz * px - lx * pz, cz = lx * py - ly * px;
const d2 = cx * cx + cy * cy + cz * cz;
if (d2 > maxD2) {
maxD2 = d2;
p2 = i;
}
}
if (p2 === -1 || maxD2 < 1e-9) return [];
// Furthest from plane minX-maxX-p2
const pA = pts[minX], pB = pts[maxX], pC = pts[p2];
const nx = (pB.y - pA.y) * (pC.z - pA.z) - (pB.z - pA.z) * (pC.y - pA.y);
const ny = (pB.z - pA.z) * (pC.x - pA.x) - (pB.x - pA.x) * (pC.z - pA.z);
const nz = (pB.x - pA.x) * (pC.y - pA.y) - (pB.y - pA.y) * (pC.x - pA.x);
let maxDPlane = 0, p3 = -1;
for (let i = 0; i < pts.length; i++) {
if (i === minX || i === maxX || i === p2) continue;
const d = Math.abs(nx * (pts[i].x - pA.x) + ny * (pts[i].y - pA.y) + nz * (pts[i].z - pA.z));
if (d > maxDPlane) {
maxDPlane = d;
p3 = i;
}
}
if (p3 === -1 || maxDPlane < 1e-9) return [];
// Build initial tetrahedron with outward-pointing normals
const p0 = pts[minX], p1 = pts[maxX], pt2 = pts[p2], pt3 = pts[p3];
const center = {
x: (p0.x + p1.x + pt2.x + pt3.x) / 4,
y: (p0.y + p1.y + pt2.y + pt3.y) / 4,
z: (p0.z + p1.z + pt2.z + pt3.z) / 4,
};
function createFace(a, b, c) {
const fnx = (b.y - a.y) * (c.z - a.z) - (b.z - a.z) * (c.y - a.y);
const fny = (b.z - a.z) * (c.x - a.x) - (b.x - a.x) * (c.z - a.z);
const fnz = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
const cpx = center.x - a.x, cpy = center.y - a.y, cpz = center.z - a.z;
if (fnx * cpx + fny * cpy + fnz * cpz > 0) {
return { a: a, b: c, c: b, active: true };
}
return { a: a, b: b, c: c, active: true };
}
let faces = [
createFace(p0, p1, pt2),
createFace(p0, pt2, pt3),
createFace(p0, pt3, p1),
createFace(p1, pt3, pt2),
];
// Incrementally add remaining points
for (let i = 0; i < pts.length; i++) {
if (i === minX || i === maxX || i === p2 || i === p3) continue;
const pt = pts[i];
// Find all visible faces
const vis = [];
for (let f = 0; f < faces.length; f++) {
if (faces[f].active && visible(faces[f], pt)) {
vis.push(f);
}
}
if (vis.length === 0) continue;
// Find horizon edges
const edgeMap = new Map();
for (const fIdx of vis) {
const f = faces[fIdx];
const edges = [
{ u: f.a, v: f.b },
{ u: f.b, v: f.c },
{ u: f.c, v: f.a },
];
for (const e of edges) {
const uId = pts.indexOf(e.u), vId = pts.indexOf(e.v);
const key = uId < vId ? (uId + '_' + vId) : (vId + '_' + uId);
const current = edgeMap.get(key) || { count: 0, u: e.u, v: e.v };
current.count++;
edgeMap.set(key, current);
}
faces[fIdx].active = false;
}
// Create new faces from horizon edges to pt
for (const entry of edgeMap.values()) {
if (entry.count === 1) {
const newF = createFace(entry.u, entry.v, pt);
faces.push(newF);
}
}
}
return faces.filter(f => f.active);
}