chore(gamut): modularize QuickHull 3D under src/js/vendor/quickhull.js (fixes #117) #125
+12
-257
@@ -1,3 +1,4 @@
|
||||
import { computeQuickHull } from "./vendor/quickhull.js";
|
||||
const { invoke } = window.__TAURI__.core;
|
||||
const { listen } = window.__TAURI__.event;
|
||||
|
||||
@@ -119,267 +120,21 @@ export function parseCGATS(text) {
|
||||
* Output: { vertices: Float32Array, indices: Uint32Array } or null
|
||||
*/
|
||||
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 indices = [];
|
||||
const ptMap = new Map();
|
||||
|
||||
for (const f of activeFaces) {
|
||||
for (const pIdx of [f.a, f.b, f.c]) {
|
||||
if (!usedPoints.has(pIdx)) {
|
||||
usedPoints.set(pIdx, verticesList.length / 3);
|
||||
verticesList.push(points[pIdx].x, points[pIdx].y, points[pIdx].z);
|
||||
for (const f of faces) {
|
||||
for (const p of [f.a, f.b, f.c]) {
|
||||
const key = `${p.x}_${p.y}_${p.z}`;
|
||||
if (!ptMap.has(key)) {
|
||||
ptMap.set(key, verticesList.length / 3);
|
||||
verticesList.push(p.x, p.y, p.z);
|
||||
}
|
||||
indices.push(usedPoints.get(pIdx));
|
||||
indices.push(ptMap.get(key));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Vendored
+150
@@ -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);
|
||||
}
|
||||
Reference in New Issue
Block a user