diff --git a/src/js/gamut_viewer.js b/src/js/gamut_viewer.js index bc2e0f5..ad45f5f 100644 --- a/src/js/gamut_viewer.js +++ b/src/js/gamut_viewer.js @@ -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)); } } diff --git a/src/js/vendor/quickhull.js b/src/js/vendor/quickhull.js new file mode 100644 index 0000000..bdd3cb2 --- /dev/null +++ b/src/js/vendor/quickhull.js @@ -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); +}