diff --git a/index.html b/index.html index cd92e74..69cb379 100644 --- a/index.html +++ b/index.html @@ -10,11 +10,28 @@
+ + + + + -
@@ -24,27 +41,46 @@

- - + + +

+

+ +

- + +

- + +

- - + +

+ +
+

+ + +

+ + +
+

diff --git a/main.js b/main.js index c6cc431..5ad8808 100644 --- a/main.js +++ b/main.js @@ -1,119 +1,214 @@ -import { Niivue } from '@niivue/niivue'; -import { createMZ3, downloadMesh } from './marching-cubes.js'; +import { Niivue, NVMeshUtilities } from '@niivue/niivue' +// Dynamically load the worker modules +async function createVoxelWorker() { + if (window.Worker) { + const { default: VoxelWorker } = await import('./voxelWorker?worker') + return new VoxelWorker() + } + throw new Error('Web Workers are not supported in this environment.') +} -// Dynamically load the worker module -async function createWorker() { +async function createMeshWorker() { if (window.Worker) { - const { default: MeshWorker } = await import('./meshWorker?worker'); - return new MeshWorker(); + const { default: MeshWorker } = await import('./meshWorker?worker') + return new MeshWorker() } - throw new Error('Web Workers are not supported in this environment.'); + throw new Error('Web Workers are not supported in this environment.') } +function formatNumber(value) { + if (Math.abs(value) >= 1) { + // For numbers >= 1, use up to 1 decimal place + return value.toFixed(1) + } else { + // For numbers < 1, use up to 3 significant digits + return value.toPrecision(3) + } +} async function main() { - const worker = await createWorker(); - const loadingCircle = document.getElementById('loadingCircle'); - - worker.onmessage = async function (e) { - const { vertices, triangles } = e.data; - - const verticesArray = new Float32Array(vertices); - const trianglesArray = new Uint32Array(triangles); - - const meshBuffer = createMZ3(verticesArray, trianglesArray, false); - await nv1.loadFromArrayBuffer(meshBuffer, 'test.mz3'); - loadingCircle.classList.add('hidden'); - }; - + const MeshWorker = await createMeshWorker() + const VoxelWorker = await createVoxelWorker() + const loadingCircle = document.getElementById('loadingCircle') + function meshStatus() { + const str = `Mesh has ${nv1.meshes[0].pts.length / 3} vertices and ${nv1.meshes[0].tris.length / 3} triangles` + document.getElementById('location').innerHTML = str + shaderSelect.onchange() + nv1.setMeshProperty(nv1.meshes[0].id, 'visible', visibleCheck.checked) + } + async function loadMesh(vertices, triangles) { + if (nv1.meshes.length > 0) { + nv1.removeMesh(nv1.meshes[0]) + } + const verticesArray = new Float32Array(vertices) + const trianglesArray = new Uint32Array(triangles) + const meshBuffer = NVMeshUtilities.createMZ3(verticesArray, trianglesArray, false) + await nv1.loadFromArrayBuffer(meshBuffer, 'test.mz3') + loadingCircle.classList.add('hidden') + meshStatus() + } + MeshWorker.onmessage = async function (e) { + const { vertices, triangles } = e.data + await loadMesh(vertices, triangles) + } + VoxelWorker.onmessage = async function (e) { + const { vertices, triangles } = e.data + await loadMesh(vertices, triangles) + } saveBtn.onclick = function () { if (nv1.meshes.length < 1) { - window.alert("No mesh open for saving. Use 'Create Mesh'."); + window.alert("No mesh open for saving. Use 'Create Mesh'.") } else { - saveDialog.show(); + saveDialog.show() } - }; - + } + volumeSelect.onchange = function () { + const selectedOption = volumeSelect.options[volumeSelect.selectedIndex] + const txt = selectedOption.text + let fnm = './' + txt + if (volumeSelect.selectedIndex > 3) { + fnm = 'https://niivue.github.io/niivue/images/' + txt + } else if (volumeSelect.selectedIndex > 0) { + fnm = 'https://niivue.github.io/niivue-demo-images/' + txt + } + if (nv1.meshes.length > 0) { + nv1.removeMesh(nv1.meshes[0]) + } + if (nv1.volumes.length > 0) { + nv1.removeVolumeByIndex(0) + } + if (volumeSelect.selectedIndex > 4) { + nv1.loadMeshes([{ url: fnm }]) + } else { + if (!fnm.endsWith('.mgz')) { + fnm += '.nii.gz' + } + nv1.loadVolumes([{ url: fnm }]) + } + } applySaveBtn.onclick = function () { if (nv1.meshes.length < 1) { - return; - } - if (formatSelect.selectedIndex == 0) { - downloadMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, 'simplified_mesh.mz3', true); + return } - if (formatSelect.selectedIndex == 1) { - downloadMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, 'simplified_mesh.obj'); + let format = 'obj' + if (formatSelect.selectedIndex === 0) { + format = 'mz3' } - if (formatSelect.selectedIndex == 2) { - downloadMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, 'simplified_mesh.stl'); + if (formatSelect.selectedIndex === 2) { + format = 'stl' } - }; - + NVMeshUtilities.saveMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, `mesh.${format}`, true) + } remeshBtn.onclick = function () { - remeshDialog.show(); - }; - - applyBtn.onclick = async function () { - if (nv1.meshes.length > 0) { - nv1.removeMesh(nv1.meshes[0]); + if (nv1.volumes.length < 1) { + window.alert('No voxel-based image open for meshing. Drag and drop an image.') + } else { + remeshDialog.show() } - const img = new Uint8ClampedArray(nv1.volumes[0].img); - const dims = [ - nv1.volumes[0].hdr.dims[1], - nv1.volumes[0].hdr.dims[2], - nv1.volumes[0].hdr.dims[3] - ]; - const isoValue = isoSlide.value; - const largestCheckValue = largestCheck.checked; - const bubbleCheckValue = bubbleCheck.checked; - const affine = nv1.volumes[0].hdr.affine; - const shrinkValue = shrinkSlide.value / shrinkSlide.max; - - loadingCircle.classList.remove('hidden'); - - worker.postMessage({ - img: img.buffer, - dims, - isoValue, - largestCheck: largestCheckValue, - bubbleCheck: bubbleCheckValue, - affine, - shrinkValue - }, [img.buffer]); - }; - - aboutBtn.onclick = function () { + } + simplifyBtn.onclick = function () { if (nv1.meshes.length < 1) { - window.alert("No mesh open for saving. Use 'Create Mesh'."); + window.alert('No mesh open to simplify. Drag and drop a mesh or create a mesh from a voxel based image.') } else { - window.alert( - `Mesh has ${nv1.meshes[0].pts.length / 3} vertices and ${nv1.meshes[0].tris.length / 3} triangles` - ); + simplifyDialog.show() } - }; - + } + applySimpleBtn.onclick = async function () { + if (nv1.meshes.length < 1) { + console.log('No mesh open to simplify.') + return + } + const verts = nv1.meshes[0].pts.slice() + const tris = nv1.meshes[0].tris.slice() + const shrinkValue = Math.min(Math.max(Number(shrinkSimplePct.value) / 100, 0.01), 1) + loadingCircle.classList.remove('hidden') + MeshWorker.postMessage({ + verts, + tris, + shrinkValue + }) + } + applyBtn.onclick = async function () { + if (nv1.volumes.length < 1) { + console.log('No volume open to meshify.') + return + } + const isoValue = Number(isoNumber.value) + let img = new Float32Array(nv1.volumes[0].img) + const scl_slope = nv1.volumes[0].hdr.scl_slope + const scl_inter = nv1.volumes[0].hdr.scl_inter + if (scl_slope !== 1.0 || scl_inter !== 0) { + img = new Float32Array(img) + for (let i = 0; i < img.length; i++) { + img[i] = img[i] * scl_slope + scl_inter + } + } + const dims = [nv1.volumes[0].hdr.dims[1], nv1.volumes[0].hdr.dims[2], nv1.volumes[0].hdr.dims[3]] + const largestCheckValue = largestCheck.checked + const bubbleCheckValue = bubbleCheck.checked + const affine = nv1.volumes[0].hdr.affine + const shrinkValue = Math.min(Math.max(Number(shrinkPct.value) / 100, 0.01), 1) + const smoothValue = smoothSlide.value + loadingCircle.classList.remove('hidden') + VoxelWorker.postMessage( + { + img, + dims, + isoValue, + largestCheck: largestCheckValue, + bubbleCheck: bubbleCheckValue, + smoothValue, + affine, + shrinkValue + }, + [img.buffer] + ) + } visibleCheck.onchange = function () { - nv1.setMeshProperty(nv1.meshes[0].id, 'visible', this.checked); - }; - + nv1.setMeshProperty(nv1.meshes[0].id, 'visible', this.checked) + } function handleLocationChange(data) { - document.getElementById('location').innerHTML = '  ' + data.string; + document.getElementById('location').innerHTML = '  ' + data.string + } + shaderSelect.onchange = function () { + nv1.setMeshShader(nv1.meshes[0].id, this.value) + } + function handleMeshLoaded() { + meshStatus() } - const defaults = { + onMeshLoaded: handleMeshLoaded, onLocationChange: handleLocationChange, + backColor: [0.2, 0.2, 0.3, 1], show3Dcrosshair: true - }; - const nv1 = new Niivue(defaults); - nv1.attachToCanvas(gl1); - nv1.setClipPlane([0.2, 0, 120]); - nv1.opts.dragMode = nv1.dragModes.pan; - nv1.setRenderAzimuthElevation(245, 15); - nv1.opts.multiplanarForceRender = true; - nv1.opts.yoke3Dto2DZoom = true; - nv1.opts.crosshairGap = 11; - nv1.setInterpolation(true); - await nv1.loadVolumes([{ url: './bet.nii.gz' }]); - applyBtn.onclick(); + } + const nv1 = new Niivue(defaults) + nv1.attachToCanvas(gl1) + nv1.isAlphaClipDark = true + function imageStatus() { + const otsu = nv1.findOtsu(3) + isoLabel.textContent = + 'Isosurface Threshold (' + + formatNumber(nv1.volumes[0].cal_min) + + '...' + + formatNumber(nv1.volumes[0].cal_max) + + ')' + isoNumber.value = formatNumber(otsu[1]) + const str = `Image has ${nv1.volumes[0].dims[1]}×${nv1.volumes[0].dims[2]}×${nv1.volumes[0].dims[3]} voxels` + document.getElementById('location').innerHTML = str + nv1.setSliceType(nv1.sliceTypeMultiplanar) + } + nv1.onImageLoaded = () => { + imageStatus() + } + nv1.setClipPlane([0.2, 0, 120]) + nv1.opts.dragMode = nv1.dragModes.pan + nv1.setRenderAzimuthElevation(245, 15) + nv1.opts.multiplanarForceRender = true + nv1.opts.yoke3Dto2DZoom = true + nv1.setInterpolation(true) + await nv1.loadVolumes([{ url: './bet.nii.gz' }]) + imageStatus() + applyBtn.onclick() } -main(); +main() diff --git a/main_no_web_worker.js b/main_no_web_worker.js deleted file mode 100644 index 8f432ec..0000000 --- a/main_no_web_worker.js +++ /dev/null @@ -1,82 +0,0 @@ -import { Niivue } from '@niivue/niivue' -import { voxels2mesh, createMZ3, downloadMesh } from './marching-cubes.js' -import { simplifyJS } from './simplify.js' - -async function main() { - saveBtn.onclick = function () { - if (nv1.meshes.length < 1) { - window.alert("No mesh open for saving. Use 'Create Mesh'.") - } else { - // downloadArrayBuffer(meshBuffer, 'simplified_mesh.mz3') - saveDialog.show() - } - } - applySaveBtn.onclick = function () { - if (nv1.meshes.length < 1) { - return - } - if (formatSelect.selectedIndex == 0) { - downloadMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, 'simplified_mesh.mz3', true) - } - if (formatSelect.selectedIndex == 1) { - downloadMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, 'simplified_mesh.obj') - } - if (formatSelect.selectedIndex == 2) { - downloadMesh(nv1.meshes[0].pts, nv1.meshes[0].tris, 'simplified_mesh.stl') - } - } - remeshBtn.onclick = function () { - remeshDialog.show() - } - applyBtn.onclick = async function () { - if (nv1.meshes.length > 0) { - nv1.removeMesh(nv1.meshes[0]) - } - console.log(nv1.volumes[0]) - const img = new Uint8ClampedArray(nv1.volumes[0].img) - const dims = [nv1.volumes[0].hdr.dims[1], nv1.volumes[0].hdr.dims[2], nv1.volumes[0].hdr.dims[3]] - let mesh = voxels2mesh( - img, - dims, - isoSlide.value, - largestCheck.checked, - bubbleCheck.checked, - nv1.volumes[0].hdr.affine - ) - mesh = simplifyJS(mesh.vertices, mesh.triangles, shrinkSlide.value / shrinkSlide.max) - // let meshBuffer = createMZ3(mesh.vertices, mesh.triangles, false) - const meshBuffer = downloadMesh(mesh.vertices, mesh.triangles, '.mz3') - await nv1.loadFromArrayBuffer(meshBuffer, 'test.mz3') - // nv1.setMeshShader(nv1.meshes[0].id, "Edge") - } - aboutBtn.onclick = function () { - if (nv1.meshes.length < 1) { - window.alert("No mesh open for saving. Use 'Create Mesh'.") - } else { - window.alert(`Mesh has ${nv1.meshes[0].pts.length / 3} vertices and ${nv1.meshes[0].tris.length / 3} triangles`) - } - } - visibleCheck.onchange = function () { - nv1.setMeshProperty(nv1.meshes[0].id, 'visible', this.checked) - } - function handleLocationChange(data) { - document.getElementById('location').innerHTML = '  ' + data.string - } - const defaults = { - onLocationChange: handleLocationChange, - show3Dcrosshair: true - } - const nv1 = new Niivue(defaults) - nv1.attachToCanvas(gl1) - nv1.setClipPlane([0.2, 0, 120]) - nv1.opts.dragMode = nv1.dragModes.pan - nv1.setRenderAzimuthElevation(245, 15) - nv1.opts.multiplanarForceRender = true - nv1.opts.yoke3Dto2DZoom = true - nv1.opts.crosshairGap = 11 - nv1.setInterpolation(true) - await nv1.loadVolumes([{ url: './bet.nii.gz' }]) - applyBtn.onclick() -} - -main() diff --git a/marching-cubes.js b/marching-cubes.js index cb686f0..096ac86 100644 --- a/marching-cubes.js +++ b/marching-cubes.js @@ -1,4 +1,3 @@ -import { gzipSync } from 'fflate' import { BWLabeler } from './bwlabels.js' // https://raw.githubusercontent.com/Twinklebear/webgl-marching-cubes/master/js/marching-cubes.js @@ -448,7 +447,7 @@ function sortVerticesByDistance(vertices) { return distances } -function unifyVertices(vertices, triangles, verbose = false) { +function unifyVertices(vertices, triangles, isVerbose = false) { // Extract the first vertex const nVtx = vertices.length / 3 const xyz0 = [vertices[0], vertices[1], vertices[2]] @@ -494,7 +493,7 @@ function unifyVertices(vertices, triangles, verbose = false) { nVtxUnique++ } if (nVtx === nVtxUnique) { - if (verbose) { + if (isVerbose) { console.log('Unify vertices found no shared vertices') } return { vertices, triangles } @@ -512,13 +511,35 @@ function unifyVertices(vertices, triangles, verbose = false) { for (let i = 0; i < triangles.length; i++) { triangles[i] = triRemaps[triangles[i]] } - if (verbose) { + if (isVerbose) { console.log(`Vertex welding ${nVtx} -> ${nVtxUnique}`) } return { vertices: newVertices, triangles } } -function verticesVox2mm(vertices, affine) { +function isPositiveDeterminant(matrix) { + // Calculate the determinant of the 3x3 matrix + const det = + matrix[0][0] * (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) - + matrix[0][1] * (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]) + + matrix[0][2] * (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) + + // Return true if the determinant is positive, otherwise false + return det > 0 +} + +function reverseMeshFaces(tris) { + for (let i = 0; i < tris.length; i += 3) { + // Swap the second and third elements of each triangle + const temp = tris[i + 1] + tris[i + 1] = tris[i + 2] + tris[i + 2] = temp + } + return tris +} + +function verticesVox2mm(mesh, affine) { + const vertices = mesh.vertices for (let i = 0; i < vertices.length; i += 3) { const xyz = [vertices[i], vertices[i + 1], vertices[i + 2]] // translate 0.5 voxels: NIfTI SForm from voxel centers @@ -529,10 +550,14 @@ function verticesVox2mm(vertices, affine) { vertices[i + 1] = xyz[0] * affine[1][0] + xyz[1] * affine[1][1] + xyz[2] * affine[1][2] + affine[1][3] vertices[i + 2] = xyz[0] * affine[2][0] + xyz[1] * affine[2][1] + xyz[2] * affine[2][2] + affine[2][3] } - return vertices + const isPos = isPositiveDeterminant(affine) + if (!isPos) { + mesh.triangles = reverseMeshFaces(mesh.triangles) + } + return mesh } -function cullDegenerateTriangles(triangles, verbose = true) { +function cullDegenerateTriangles(triangles, isVerbose = true) { let tris = triangles const nTri = tris.length / 3 let nOK = 0 @@ -548,7 +573,7 @@ function cullDegenerateTriangles(triangles, verbose = true) { if (nOK === 0) { throw new Error('invalid mesh: no valid triangles') } - if (verbose) { + if (isVerbose) { console.log(`${nTri - nOK} of the ${nTri} triangles are degenerate`) } tris = triangles.slice() @@ -657,7 +682,7 @@ function fillh(imgBin, dim, is26 = true) { return imgBin } -function fillBubbles(img, dims, isovalue = 1) { +function fillBubbles(img, dims, isovalue = 1, isVerbose = false) { // TODO check with binary input: isovalue might be off by 1 let img01 = img.map((value) => (value >= isovalue ? 1 : 0)) const countNonZero = img01.reduce((count, value) => count + (value !== 0 ? 1 : 0), 0) @@ -666,7 +691,9 @@ function fillBubbles(img, dims, isovalue = 1) { if (countNonZero === countFilledNonZero) { return img } - console.log(`${countFilledNonZero - countNonZero} voxels were bubbles`) + if (isVerbose) { + console.log(`${countFilledNonZero - countNonZero} voxels were bubbles`) + } const isovalue1 = isovalue + 1 for (let i = 0; i < img.length; i++) { if (img01[i] === 0) { @@ -727,30 +754,158 @@ function zeroVolumeEdges(img, dims) { return img } +function laplacian_smooth(verts, tris) { + const nvert = verts.length / 3 + const ntri = tris.length / 3 + const sum = new Float64Array(nvert * 3).fill(0) + const num = new Float64Array(nvert).fill(0) + + for (let i = 0; i < ntri; i++) { + // Each point of a triangle has two neighbors: + const p0 = tris[i * 3] + const p1 = tris[i * 3 + 1] + const p2 = tris[i * 3 + 2] + + num[p0] += 2 + sum[p0 * 3] += verts[p1 * 3] + verts[p2 * 3] + sum[p0 * 3 + 1] += verts[p1 * 3 + 1] + verts[p2 * 3 + 1] + sum[p0 * 3 + 2] += verts[p1 * 3 + 2] + verts[p2 * 3 + 2] + + num[p1] += 2 + sum[p1 * 3] += verts[p0 * 3] + verts[p2 * 3] + sum[p1 * 3 + 1] += verts[p0 * 3 + 1] + verts[p2 * 3 + 1] + sum[p1 * 3 + 2] += verts[p0 * 3 + 2] + verts[p2 * 3 + 2] + + num[p2] += 2 + sum[p2 * 3] += verts[p0 * 3] + verts[p1 * 3] + sum[p2 * 3 + 1] += verts[p0 * 3 + 1] + verts[p1 * 3 + 1] + sum[p2 * 3 + 2] += verts[p0 * 3 + 2] + verts[p1 * 3 + 2] + } + + for (let i = 0; i < nvert; i++) { + if (num[i] <= 0) { + continue + } + verts[i * 3] = sum[i * 3] / num[i] + verts[i * 3 + 1] = sum[i * 3 + 1] / num[i] + verts[i * 3 + 2] = sum[i * 3 + 2] / num[i] + } +} + +function laplacian_smoothHC(verts, tris, iter = 100, alpha = 0.1, beta = 0.5, lockEdges = true) { + if (iter < 1) { + return + } + const nvert = verts.length / 3 + const ntri = tris.length / 3 + const alpha1 = 1.0 - alpha + const beta1 = 1.0 - beta + + // let p = new Array(nvert * 3).fill(0); + const q = new Float32Array(nvert * 3).fill(0) + const b = new Array(nvert * 3).fill(0) + // Copy verts to p + // p.set(verts); + const p = new Float32Array(verts) + + for (let j = 0; j < iter; j++) { + // Copy p to q + q.set(p) + + // Laplacian smooth on p + laplacian_smooth(p, tris) + + for (let i = 0; i < nvert; i++) { + const pi = i * 3 + const vi = i * 3 + + b[pi] = p[pi] - (alpha * verts[vi] + alpha1 * q[pi]) + b[pi + 1] = p[pi + 1] - (alpha * verts[vi + 1] + alpha1 * q[pi + 1]) + b[pi + 2] = p[pi + 2] - (alpha * verts[vi + 2] + alpha1 * q[pi + 2]) + } + + // Copy b to q + q.set(b) + + // Laplacian smooth on q + laplacian_smooth(q, tris) + + for (let i = 0; i < nvert; i++) { + const pi = i * 3 + + p[pi] = p[pi] - (beta * b[pi] + beta1 * q[pi]) + p[pi + 1] = p[pi + 1] - (beta * b[pi + 1] + beta1 * q[pi + 1]) + p[pi + 2] = p[pi + 2] - (beta * b[pi + 2] + beta1 * q[pi + 2]) + } + } + + if (!lockEdges) { + verts.set(p) + return + } + + // Find border edges and only update non-border vertices + const vertices = new Array(nvert).fill(null).map(() => ({ border: false, p: [0, 0, 0] })) + for (let i = 0; i < nvert; i++) { + vertices[i].p = [verts[i * 3], verts[i * 3 + 1], verts[i * 3 + 2]] + } + + const triangles = new Array(ntri).fill(null).map(() => ({ deleted: 0, v: [0, 0, 0] })) + for (let i = 0; i < ntri; i++) { + triangles[i].v[0] = tris[i * 3] + triangles[i].v[1] = tris[i * 3 + 1] + triangles[i].v[2] = tris[i * 3 + 2] + } + + // Placeholder functions for mesh update logic + function update_mesh() { + // Implement mesh update logic + } + + update_mesh(0, triangles, vertices, [], 0, ntri, nvert) + + for (let i = 0; i < nvert; i++) { + if (vertices[i].border) { + continue + } + verts[i * 3] = p[i * 3] + verts[i * 3 + 1] = p[i * 3 + 1] + verts[i * 3 + 2] = p[i * 3 + 2] + } +} + export function voxels2mesh( img, dims, isovalue, isLargestClusterOnly = false, isFillBubbles = false, + smoothIterations = 5, affine = [], - verbose = true + isVerbose = true ) { - if (!(img instanceof Uint8Array) && !(img instanceof Uint8ClampedArray)) { - throw new Error('img must be a Uint8Array') - } - if (isovalue < 0 || isovalue > 255) { - throw new Error('isovalue must be in the range 0..255') - } if (dims[0] < 3 || dims[1] < 3 || dims[2] < 3) { throw new Error('volume too small for meshification') } + if (isovalue < 0) { + // marching cubes implementation assumes positive values + let mn = img[0] + for (let i = 0; i < img.length; i++) { + mn = Math.min(img[i], mn) + } + if (mn < 0) { + for (let i = 0; i < img.length; i++) { + img[i] = img[i] - mn + } + isovalue = isovalue - mn + } + } img = zeroVolumeEdges(img, dims) if (isLargestClusterOnly) { img = largestClusterOnly(img, dims, isovalue) } if (isFillBubbles) { - img = fillBubbles(img, dims, isovalue) + img = fillBubbles(img, dims, isovalue, isVerbose) } let mesh = [] mesh.vertices = marchingCubesJS(img, dims, isovalue) @@ -762,148 +917,16 @@ export function voxels2mesh( for (let i = 0; i < mesh.triangles.length; i++) { mesh.triangles[i] = i } - mesh = unifyVertices(mesh.vertices, mesh.triangles, verbose) - if (verbose) { + mesh = unifyVertices(mesh.vertices, mesh.triangles, isVerbose) + if (isVerbose) { console.log(`Vertices ${mesh.vertices.length / 3} Triangles ${mesh.triangles.length / 3}`) } // these algoritms should never yield degenerate triangles, but lets check - mesh.triangles = cullDegenerateTriangles(mesh.triangles, verbose) + mesh.triangles = cullDegenerateTriangles(mesh.triangles, isVerbose) + laplacian_smoothHC(mesh.vertices, mesh.triangles, smoothIterations) // optional affine is a 4x4 spatial transformation matrix if (affine.length === 4) { - mesh.vertices = verticesVox2mm(mesh.vertices, affine) + mesh = verticesVox2mm(mesh, affine) } return { vertices: mesh.vertices, triangles: mesh.triangles } } - -export function createMZ3(vertices, indices, compress = false) { - // generate binary MZ3 format mesh - // n.b. small, precise and small but support is not widespread - // n.b. result can be compressed with gzip - // https://github.com/neurolabusc/surf-ice/tree/master/mz3 - const magic = 23117 - const attr = 3 - const nface = indices.length / 3 - const nvert = vertices.length / 3 - const nskip = 0 - // Calculate buffer size - const headerSize = 16 - const indexSize = nface * 3 * 4 // Uint32Array - const vertexSize = nvert * 3 * 4 // Float32Array - const totalSize = headerSize + indexSize + vertexSize - const buffer = new ArrayBuffer(totalSize) - const writer = new DataView(buffer) - // Write header - writer.setUint16(0, magic, true) - writer.setUint16(2, attr, true) - writer.setUint32(4, nface, true) - writer.setUint32(8, nvert, true) - writer.setUint32(12, nskip, true) - // Write indices - let offset = headerSize - new Uint32Array(buffer, offset, indices.length).set(indices) - offset += indexSize - // Write vertices - new Float32Array(buffer, offset, vertices.length).set(vertices) - if (compress) { - return gzipSync(new Uint8Array(buffer)) - } - return buffer -} - -function createOBJ(vertices, indices) { - // generate binary OBJ format mesh - // n.b. widespread support, but large and slow due to ASCII - // https://en.wikipedia.org/wiki/Wavefront_.obj_file - let objContent = '' - // Add vertices to OBJ content - for (let i = 0; i < vertices.length; i += 3) { - objContent += `v ${vertices[i]} ${vertices[i + 1]} ${vertices[i + 2]}\n` - } - // Add faces to OBJ content (OBJ indices start at 1, not 0) - for (let i = 0; i < indices.length; i += 3) { - objContent += `f ${indices[i] + 1} ${indices[i + 1] + 1} ${indices[i + 2] + 1}\n` - } - // Encode the OBJ content as an ArrayBuffer - const encoder = new TextEncoder() - const arrayBuffer = encoder.encode(objContent).buffer - return arrayBuffer -} - -function createSTL(vertices, indices) { - // generate binary STL format mesh - // n.b. inefficient and slow as vertices are not reused - // https://en.wikipedia.org/wiki/STL_(file_format)#Binary - const numTriangles = indices.length / 3 - const bufferLength = 84 + numTriangles * 50 - const arrayBuffer = new ArrayBuffer(bufferLength) - const dataView = new DataView(arrayBuffer) - // Write header (80 bytes) - for (let i = 0; i < 80; i++) { - dataView.setUint8(i, 0) - } - // Write number of triangles (4 bytes) - dataView.setUint32(80, numTriangles, true) - let offset = 84 - for (let i = 0; i < indices.length; i += 3) { - const i0 = indices[i] * 3 - const i1 = indices[i + 1] * 3 - const i2 = indices[i + 2] * 3 - // Normal vector (12 bytes, set to zero) - dataView.setFloat32(offset, 0, true) // Normal X - dataView.setFloat32(offset + 4, 0, true) // Normal Y - dataView.setFloat32(offset + 8, 0, true) // Normal Z - offset += 12 - // Vertex 1 (12 bytes) - dataView.setFloat32(offset, vertices[i0], true) // Vertex 1 X - dataView.setFloat32(offset + 4, vertices[i0 + 1], true) // Vertex 1 Y - dataView.setFloat32(offset + 8, vertices[i0 + 2], true) // Vertex 1 Z - offset += 12 - // Vertex 2 (12 bytes) - dataView.setFloat32(offset, vertices[i1], true) // Vertex 2 X - dataView.setFloat32(offset + 4, vertices[i1 + 1], true) // Vertex 2 Y - dataView.setFloat32(offset + 8, vertices[i1 + 2], true) // Vertex 2 Z - offset += 12 - // Vertex 3 (12 bytes) - dataView.setFloat32(offset, vertices[i2], true) // Vertex 3 X - dataView.setFloat32(offset + 4, vertices[i2 + 1], true) // Vertex 3 Y - dataView.setFloat32(offset + 8, vertices[i2 + 2], true) // Vertex 3 Z - offset += 12 - // Attribute byte count (2 bytes, set to zero) - dataView.setUint16(offset, 0, true) - offset += 2 - } - return arrayBuffer -} - -function downloadArrayBuffer(buffer, filename) { - const blob = new Blob([buffer], { type: 'application/octet-stream' }) - const url = URL.createObjectURL(blob) - const a = document.createElement('a') - a.href = url - a.download = filename - document.body.appendChild(a) - a.style.display = 'none' - a.click() - setTimeout(() => { - document.body.removeChild(a) - URL.revokeObjectURL(url) - }, 0) -} - -export function downloadMesh(vertices, indices, filename, compress = false) { - let buff = [] - if (/\.obj$/i.test(filename)) { - buff = createOBJ(vertices, indices) - } else if (/\.stl$/i.test(filename)) { - buff = createSTL(vertices, indices) - } else { - if (!/\.mz3$/i.test(filename)) { - filename += '.mz3' - } - buff = createMZ3(vertices, indices, compress) - } - if (filename.length > 4) { - downloadArrayBuffer(buff, filename) - } - return buff -} diff --git a/meshWorker.js b/meshWorker.js index 8f49b6d..d2abe5a 100644 --- a/meshWorker.js +++ b/meshWorker.js @@ -1,15 +1,14 @@ -import { voxels2mesh, createMZ3, downloadMesh } from './marching-cubes.js'; -import { simplifyJS } from './simplify.js'; +import { simplifyJS } from './simplify.js' self.onmessage = function (e) { - const { img, dims, isoValue, largestCheck, bubbleCheck, affine, shrinkValue } = e.data; - const imgArray = new Uint8ClampedArray(img); - - let mesh = voxels2mesh(imgArray, dims, isoValue, largestCheck, bubbleCheck, affine); - mesh = simplifyJS(mesh.vertices, mesh.triangles, shrinkValue); - + const { verts, tris, shrinkValue, verbose = true } = e.data + const startTime = new Date() + const mesh = simplifyJS(verts, tris, shrinkValue) + if (verbose) { + console.log(new Date() - startTime + 'ms elapsed') + } postMessage({ vertices: mesh.vertices, triangles: mesh.triangles - }); -}; + }) +} diff --git a/package-lock.json b/package-lock.json index 82c34e0..06c90a4 100644 --- a/package-lock.json +++ b/package-lock.json @@ -1,14 +1,14 @@ { - "name": "niivue-brainchop", + "name": "niivue-mesh", "version": "0.1.0", "lockfileVersion": 3, "requires": true, "packages": { "": { - "name": "niivue-brainchop", + "name": "niivue-mesh", "version": "0.1.0", "dependencies": { - "@niivue/niivue": "^0.43.6", + "@niivue/niivue": "^0.44.1", "jslint": "^0.12.1" }, "devDependencies": { @@ -403,9 +403,9 @@ } }, "node_modules/@niivue/niivue": { - "version": "0.43.6", - "resolved": "https://registry.npmjs.org/@niivue/niivue/-/niivue-0.43.6.tgz", - "integrity": "sha512-mN+R+CcHSB20MkjPMwK9rvX+f5B8gJaDqcnyZzO0AhYJYjVxq44HYFvsKEqNANfuVtdUVZOHor2wWuhBjQ+FFg==", + "version": "0.44.1", + "resolved": "https://registry.npmjs.org/@niivue/niivue/-/niivue-0.44.1.tgz", + "integrity": "sha512-cJ3x6h9p3q7glMruit7M6KqFt0DPDI6+w4WeGh59Ly+jbW2EcBpe6JfRAKR8CcsxB21ePakGYRydPBRtGVEAug==", "license": "BSD-2-Clause", "dependencies": { "@lukeed/uuid": "^2.0.1", diff --git a/package.json b/package.json index 6e4e68c..0bd609f 100644 --- a/package.json +++ b/package.json @@ -1,5 +1,5 @@ { - "name": "niivue-brainchop", + "name": "niivue-mesh", "private": true, "version": "0.1.0", "type": "module", @@ -9,7 +9,7 @@ "preview": "vite preview" }, "dependencies": { - "@niivue/niivue": "^0.43.6", + "@niivue/niivue": "^0.44.1", "jslint": "^0.12.1" }, "devDependencies": { diff --git a/simplify.js b/simplify.js index 0b8c0d2..56e191d 100644 --- a/simplify.js +++ b/simplify.js @@ -13,7 +13,7 @@ class QuadricSimplifyMesh { constructor(vs, ts, targetCount, aggressiveness = 7, finishLossless = false, verbose = true) { this.vertices = [] this.triangles = [] - this.refs = [] + this.refs = new Array(ts.length) // Pre-allocate memory this.targetCount = targetCount this.aggressiveness = aggressiveness this.finishLossless = finishLossless @@ -27,7 +27,7 @@ class QuadricSimplifyMesh { p: { x: vs[i], y: vs[i + 1], z: vs[i + 2] }, tstart: 0, tcount: 0, - q: Array(10).fill(0), + q: new Float32Array(10).fill(0), border: 0 }) } @@ -35,7 +35,7 @@ class QuadricSimplifyMesh { for (let i = 0; i < ts.length; i += 3) { this.triangles.push({ v: [ts[i], ts[i + 1], ts[i + 2]], - err: Array(4).fill(0), + err: new Float32Array(4).fill(0), dirty: false, deleted: false, n: { x: 0, y: 0, z: 0 } @@ -44,9 +44,7 @@ class QuadricSimplifyMesh { } symMat1(ret, c) { - for (let i = 0; i < 10; i++) { - ret[i] = c - } + ret.fill(c) } symMat4(ret, a, b, c, d) { @@ -76,19 +74,9 @@ class QuadricSimplifyMesh { } symMatAdd(ret, n, m) { - this.symMat10( - ret, - n[0] + m[0], - n[1] + m[1], - n[2] + m[2], - n[3] + m[3], - n[4] + m[4], - n[5] + m[5], - n[6] + m[6], - n[7] + m[7], - n[8] + m[8], - n[9] + m[9] - ) + for (let i = 0; i < 10; i++) { + ret[i] = n[i] + m[i] + } } symMatDet(m, a11, a12, a13, a21, a22, a23, a31, a32, a33) { @@ -102,26 +90,22 @@ class QuadricSimplifyMesh { ) } - ptf(x, y, z) { - return { x, y, z } - } - vCross(v1, v2) { - return this.ptf(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x) + 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 } } vSum(a, b) { - return this.ptf(a.x + b.x, a.y + b.y, a.z + b.z) + return { x: a.x + b.x, y: a.y + b.y, z: a.z + b.z } } vSubtract(a, b) { - return this.ptf(a.x - b.x, a.y - b.y, a.z - b.z) + return { x: a.x - b.x, y: a.y - b.y, z: a.z - b.z } } vNormalize(v) { - let len = Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z) + const len = Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z) if (len <= 0) { - len = 0.001 + return } v.x /= len v.y /= len @@ -133,7 +117,7 @@ class QuadricSimplifyMesh { } vMult(a, v) { - return this.ptf(a.x * v, a.y * v, a.z * v) + return { x: a.x * v, y: a.y * v, z: a.z * v } } vertexError(q, x, y, z) { @@ -151,20 +135,42 @@ class QuadricSimplifyMesh { ) } + calculateErrorFast(id_v1, id_v2) { + const q = new Float32Array(10) + this.symMatAdd(q, this.vertices[id_v1].q, this.vertices[id_v2].q) + const border = this.vertices[id_v1].border + this.vertices[id_v2].border + const det = this.symMatDet(q, 0, 1, 2, 1, 4, 5, 2, 5, 7) + + if (det !== 0.0 && border === 0) { + const x = (-1.0 / det) * this.symMatDet(q, 1, 2, 3, 4, 5, 6, 5, 7, 8) + const y = (1.0 / det) * this.symMatDet(q, 0, 2, 3, 1, 5, 6, 2, 7, 8) + const z = (-1.0 / det) * this.symMatDet(q, 0, 1, 3, 1, 4, 6, 2, 5, 8) + return this.vertexError(q, x, y, z) + } + + const p1 = this.vertices[id_v1].p + const p2 = this.vertices[id_v2].p + const p3 = this.vMult(this.vSum(p1, p2), 0.5) + const error1 = this.vertexError(q, p1.x, p1.y, p1.z) + const error2 = this.vertexError(q, p2.x, p2.y, p2.z) + const error3 = this.vertexError(q, p3.x, p3.y, p3.z) + const error = Math.min(error1, Math.min(error2, error3)) + return error + } + calculateError(id_v1, id_v2, p_result) { - const q = Array(10).fill(0) + const q = new Float32Array(10) this.symMatAdd(q, this.vertices[id_v1].q, this.vertices[id_v2].q) const border = this.vertices[id_v1].border + this.vertices[id_v2].border const det = this.symMatDet(q, 0, 1, 2, 1, 4, 5, 2, 5, 7) if (det !== 0.0 && border === 0) { - p_result.x = (-1.0 / det) * this.symMatDet(q, 1, 2, 3, 4, 5, 6, 5, 7, 8) // vx = A41/det(q_delta) - p_result.y = (1.0 / det) * this.symMatDet(q, 0, 2, 3, 1, 5, 6, 2, 7, 8) // vy = A42/det(q_delta) - p_result.z = (-1.0 / det) * this.symMatDet(q, 0, 1, 3, 1, 4, 6, 2, 5, 8) // vz = A43/det(q_delta) + p_result.x = (-1.0 / det) * this.symMatDet(q, 1, 2, 3, 4, 5, 6, 5, 7, 8) + p_result.y = (1.0 / det) * this.symMatDet(q, 0, 2, 3, 1, 5, 6, 2, 7, 8) + p_result.z = (-1.0 / det) * this.symMatDet(q, 0, 1, 3, 1, 4, 6, 2, 5, 8) return this.vertexError(q, p_result.x, p_result.y, p_result.z) } - // det = 0 -> try to find best result const p1 = this.vertices[id_v1].p const p2 = this.vertices[id_v2].p const p3 = this.vMult(this.vSum(p1, p2), 0.5) @@ -173,13 +179,13 @@ class QuadricSimplifyMesh { const error3 = this.vertexError(q, p3.x, p3.y, p3.z) const error = Math.min(error1, Math.min(error2, error3)) if (error1 === error) { - Object.assign(p_result, p1) + p_result = p1 } if (error2 === error) { - Object.assign(p_result, p2) + p_result = p2 } if (error3 === error) { - Object.assign(p_result, p3) + p_result = p3 } return error } @@ -189,8 +195,7 @@ class QuadricSimplifyMesh { let dst = 0 for (let i = 0; i < this.triangles.length; i++) { if (!this.triangles[i].deleted) { - this.triangles[dst] = this.triangles[i] - dst++ + this.triangles[dst++] = this.triangles[i] } } this.triangles.length = dst @@ -214,16 +219,6 @@ class QuadricSimplifyMesh { vertex.tcount = 0 } - /* let maxRef = 0; - for (let i = 0; i < this.triangles.length; i++) { - let t = this.triangles[i]; - for (let j = 0; j < 3; j++) { - let v = this.vertices[t.v[j]]; - maxRef = Math.max(v.tstart + v.tcount, maxRef); - } - } -console.log("mork >>", maxRef, this.triangles.length * 3) */ - // Write References this.refs.length = this.triangles.length * 3 for (let i = 0; i < this.triangles.length; i++) { const t = this.triangles[i] @@ -233,6 +228,7 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ v.tcount++ } } + if (iteration !== 0) { return } @@ -241,8 +237,8 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ vertex.border = 0 } - const vids = new Array(this.vertices.length).fill(0) - const vcount = new Array(this.vertices.length).fill(0) + const vids = new Uint32Array(this.vertices.length) + const vcount = new Uint32Array(this.vertices.length) for (let i = 0; i < this.vertices.length; i++) { let nvcount = 0 @@ -289,7 +285,7 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ this.vNormalize(n) t.n = n for (let j = 0; j < 3; j++) { - const q = Array(10).fill(0) + const q = new Float32Array(10) this.symMat4(q, n.x, n.y, n.z, -this.vDot(n, p[0])) this.symMatAdd(this.vertices[t.v[j]].q, this.vertices[t.v[j]].q, q) } @@ -297,9 +293,8 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ for (let i = 0; i < this.triangles.length; i++) { const t = this.triangles[i] - const p = {} for (let j = 0; j < 3; j++) { - t.err[j] = this.calculateError(t.v[j], t.v[(j + 1) % 3], p) + t.err[j] = this.calculateErrorFast(t.v[j], t.v[(j + 1) % 3]) } t.err[3] = Math.min(t.err[0], Math.min(t.err[1], t.err[2])) } @@ -341,7 +336,6 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ } updateTriangles(i0, v, deleted, deletedTriangles) { - const p = {} for (let k = 0; k < v.tcount; k++) { const r = this.refs[v.tstart + k] const t = this.triangles[r.tid] @@ -355,9 +349,9 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ } t.v[r.tvertex] = i0 t.dirty = true - t.err[0] = this.calculateError(t.v[0], t.v[1], p) - t.err[1] = this.calculateError(t.v[1], t.v[2], p) - t.err[2] = this.calculateError(t.v[2], t.v[0], p) + t.err[0] = this.calculateErrorFast(t.v[0], t.v[1]) + t.err[1] = this.calculateErrorFast(t.v[1], t.v[2]) + t.err[2] = this.calculateErrorFast(t.v[2], t.v[0]) t.err[3] = Math.min(t.err[0], Math.min(t.err[1], t.err[2])) this.refs.push(r) } @@ -396,7 +390,6 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ simplify(verbose = true) { let deletedTriangles = 0 - // const vertexCount = this.vertices.length const triangleCount = this.triangles.length const deleted0 = new Array(triangleCount * 3).fill(false) // overprovision const deleted1 = new Array(triangleCount * 3).fill(false) // overprovision @@ -466,7 +459,7 @@ console.log("mork >>", maxRef, this.triangles.length * 3) */ if (v0.border !== v1.border) { continue } - const p = {} + const p = { x: 0, y: 0, z: 0 } this.calculateError(i0, i1, p) if (this.flipped(p, i0, i1, v0, v1, deleted0)) { continue diff --git a/voxelWorker.js b/voxelWorker.js new file mode 100644 index 0000000..33126f2 --- /dev/null +++ b/voxelWorker.js @@ -0,0 +1,19 @@ +import { voxels2mesh } from './marching-cubes.js'; +import { simplifyJS } from './simplify.js'; + +self.onmessage = function (e) { + const { img, dims, isoValue, largestCheck, bubbleCheck, smoothValue, affine, shrinkValue, verbose = true} = e.data; + const imgArray = new Float32Array(img); + let startTime = new Date(); + let mesh = voxels2mesh(imgArray, dims, isoValue, largestCheck, bubbleCheck, smoothValue, affine, verbose); + if (shrinkValue < 1.0) { + mesh = simplifyJS(mesh.vertices, mesh.triangles, shrinkValue); + } + if (verbose) { + console.log( new Date() - startTime + "ms elapsed"); + } + postMessage({ + vertices: mesh.vertices, + triangles: mesh.triangles + }); +};