You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
202 lines
6.7 KiB
202 lines
6.7 KiB
|
|
function implicitSurfaceTriangleMesh(implicitFunction, n, args, id = 8) {
|
|
|
|
// HERE IS WHERE MOST OF THE WORK HAPPENS
|
|
|
|
let marchingTetrahedra = function(V, ni, nj) {
|
|
|
|
// CONVENIENCE FUNCTIONS TO COMPUTE (i,j,k) FROM VOLUME INDEX n
|
|
|
|
function n2i(n) { return n % ni; }
|
|
function n2j(n) { return (n / dj >>> 0) % nj; }
|
|
function n2k(n) { return n / dk >>> 0 ; }
|
|
|
|
// ADD A VERTEX, AND RETURN A UNIQUE ID FOR THAT VERTEX
|
|
|
|
function E(a, b) {
|
|
if (a > b) { let tmp = a; a = b; b = tmp; }
|
|
let ai = n2i(a), aj = n2j(a), ak = n2k(a),
|
|
bi = n2i(b), bj = n2j(b), bk = n2k(b);
|
|
let m = (n << 6) + (ai & bi ? 1 << 6 : ai | bi << 3)
|
|
+ (aj & bj ? dj << 6 : aj << 1 | bj << 4)
|
|
+ (ak & bk ? dk << 6 : ak << 2 | bk << 5);
|
|
|
|
// ADD TO VERTEX ARRAY ONLY THE FIRST TIME THE VERTEX IS ENCOUNTERED
|
|
|
|
if (vertexID[m] === undefined) {
|
|
vertexID[m] = P.length / 3;
|
|
let t = -V[n+a] / (V[n+b] - V[n+a]),
|
|
c = function(i,a,b) { return (i + (1-t)*a + t*b) / ni * 2 - 1; };
|
|
P.push( c(i,ai,bi), c(j,aj,bj), c(k,ak,bk) );
|
|
}
|
|
|
|
return vertexID[m];
|
|
}
|
|
|
|
// CASE WHERE WE ADD ONE TRIANGLE IN A TETRAHEDRON
|
|
|
|
function tri(a, b, c, d) {
|
|
T.push(E(a,b), E(a,c), E(a,d));
|
|
}
|
|
|
|
// CASE WHERE WE ADD TWO TRIANGLES IN A TETRAHEDRON
|
|
|
|
function quad(a, b, c, d) {
|
|
let ac = E(a,c), bc = E(b,c), ad = E(a,d), bd = E(b,d);
|
|
T.push(bc, ac, ad);
|
|
T.push(ad, bd, bc);
|
|
}
|
|
|
|
// DECLARE VARIABLES
|
|
|
|
let nk = V.length / (ni * nj), di = 1, dj = ni, dk = ni * nj;
|
|
let dij = di + dj, dik = di + dk, djk = dj + dk, dijk = di + dj + dk;
|
|
let P = [], T = [], vertexID = [], i, j, k, m = 0, n, S = [0,di,dij,dijk];
|
|
let lo = new Array(nj * nk),
|
|
hi = new Array(nj * nk);
|
|
|
|
// THE SIX POSSIBLE INTERMEDIATE PATHS THROUGH A TETRAHEDRON
|
|
|
|
let S1 = [di , dj , dk , di , dj , dk ];
|
|
let S2 = [dij, djk, dik, dik, dij, djk];
|
|
|
|
// THERE ARE 16 CASES TO CONSIDER
|
|
|
|
let cases = [ [0 ], [1, 0,1,2,3], [1, 1,2,0,3], [2, 0,1,2,3],
|
|
[1, 2,3,0,1], [2, 0,2,3,1], [2, 1,2,0,3], [1, 3,1,2,0],
|
|
[1, 3,0,2,1], [2, 0,3,1,2], [2, 1,3,2,0], [1, 2,1,0,3],
|
|
[2, 2,3,0,1], [1, 1,3,0,2], [1, 0,3,2,1], [0 ], ];
|
|
|
|
// FOR EACH (Y,Z), DON'T DO ANY WORK OUTSIDE OF X RANGE WHERE SURFACE MIGHT BE
|
|
|
|
for (k = 0 ; k < nk ; k++)
|
|
for (j = 0 ; j < nj ; j++, m++) {
|
|
let n0 = m * ni, n1 = n0 + ni - 1;
|
|
for (n = n0 ; n <= n1 && V[n] > 0 ; n++) ;
|
|
lo[m] = Math.max(0, n-1 - n0);
|
|
for (n = n1 ; n >= n0 && V[n] > 0 ; --n) ;
|
|
hi[m] = Math.min(ni-1, n+1 - n0);
|
|
}
|
|
|
|
// FOR ALL Y AND Z IN THE VOLUME
|
|
|
|
for (k = 0 ; k < nk - 1 ; k++) {
|
|
let i0, i1, m = k * nj, n1, s0, s1;
|
|
for (j = 0 ; j < nj - 1 ; j++, m++) {
|
|
i0 = Math.min(lo[m], lo[m+1], lo[m+ni], lo[m+1+ni]);
|
|
i1 = Math.max(hi[m], hi[m+1], hi[m+ni], hi[m+1+ni]);
|
|
|
|
// GO THROUGH RANGE OF X WHERE THE SURFACE MIGHT BE (IE: WITH ANY POSITIVE VALUES)
|
|
|
|
if (i0 <= i1) {
|
|
n = m * ni + i0;
|
|
n1 = m * ni + i1;
|
|
s0 = (V[n]>0) + (V[n+dj]>0) + (V[n+dk]>0) + (V[n+djk]>0);
|
|
for (i = i0 ; n <= n1 ; i++, n++, s0 = s1) {
|
|
|
|
// FOR EACH CUBE
|
|
|
|
s1 = (V[n+di]>0) + (V[n+dij]>0) + (V[n+dik]>0) + (V[n+dijk]>0);
|
|
if (s0 + s1 & 7) {
|
|
let C14 = (V[n] > 0) | (V[n+dijk] > 0) << 3;
|
|
|
|
// CYCLE THROUGH THE SIX TETRAHEDRA THAT TILE THE CUBE
|
|
|
|
for (let p = 0 ; p < 6 ; p++) {
|
|
let C = cases [ C14 | (V[n+S1[p]] > 0) << 1 | (V[n+S2[p]] > 0) << 2 ];
|
|
|
|
// FOR EACH TETRAHEDRON, OUTPUT EITHER ZERO, ONE OR TWO TRIANGLES
|
|
|
|
if (C[0]) { // C[0] == number of triangles to be created.
|
|
S[1] = S1[p]; // assign 2nd and 3rd corners of simplex.
|
|
S[2] = S2[p];
|
|
(C[0]==1 ? tri : quad)(S[C[1]], S[C[2]], S[C[3]], S[C[4]]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// MAKE SURE ALL TRIANGLE VERTICES ARE LISTED IN COUNTERCLOCKWISE ORDER
|
|
|
|
for (let m = 0 ; m < T.length ; m += 3) {
|
|
let a = 3 * T[m], b = 3 * T[m+1], c = 3 * T[m+2],
|
|
n = Math.floor(ni*(P[a ]+1)/2) +
|
|
Math.floor(ni*(P[a+1]+1)/2) * dj +
|
|
Math.floor(ni*(P[a+2]+1)/2) * dk,
|
|
u = cross([P[b] - P[a], P[b+1] - P[a+1], P[b+2] - P[a+2]],
|
|
[P[c] - P[b], P[c+1] - P[b+1], P[c+2] - P[b+2]]),
|
|
v = [ V[n+1] - V[n], V[n+dj] - V[n], V[n+dk] - V[n] ];
|
|
if (dot(u, v) < 0) { let tmp = T[m]; T[m] = T[m + 2]; T[m + 2] = tmp; }
|
|
}
|
|
|
|
// RETURN POINTS AND TRIANGLES
|
|
|
|
return [P, T];
|
|
}
|
|
|
|
// SAMPLE THE VOLUME
|
|
|
|
let F = i => (i - n/2) / (n/2);
|
|
let volume = [];
|
|
|
|
for (let k = 0 ; k < n ; k++)
|
|
for (let j = 0 ; j < n ; j++)
|
|
for (let i = 0 ; i < n ; i++)
|
|
volume.push(implicitFunction(F(i), F(j), F(k), args));
|
|
|
|
// FIND ALL VERTICES AND TRIANGLES IN THE VOLUME
|
|
|
|
let VT = marchingTetrahedra(volume, n, n);
|
|
let V = VT[0];
|
|
let T = VT[1];
|
|
|
|
// COMPUTE SURFACE NORMALS
|
|
|
|
let N = new Array(V.length);
|
|
for (let i = 0 ; i < V.length ; i += 3) {
|
|
let x = V[i], y = V[i+1], z = V[i+2], e = .001,
|
|
f0 = implicitFunction(x ,y ,z , args),
|
|
fx = implicitFunction(x+e,y ,z , args),
|
|
fy = implicitFunction(x ,y+e,z , args),
|
|
fz = implicitFunction(x ,y ,z+e, args),
|
|
normal = normalize([f0-fx,f0-fy,f0-fz]);
|
|
for (let j = 0 ; j < 3 ; j++)
|
|
N[i+j] = normal[j];
|
|
}
|
|
|
|
// CONSTRUCT AND RETURN THE TRIANGLES MESH
|
|
|
|
let mesh = [];
|
|
for (let i = 0; i < T.length; i += 3) {
|
|
let a = 3 * T[i ],
|
|
b = 3 * T[i + 1],
|
|
c = 3 * T[i + 2];
|
|
mesh.push( id, V[a],V[a+1],V[a+2] , N[a],N[a+1],N[a+2] ,
|
|
id, V[b],V[b+1],V[b+2] , N[b],N[b+1],N[b+2] ,
|
|
id, V[c],V[c+1],V[c+2] , N[c],N[c+1],N[c+2] );
|
|
}
|
|
return new Float32Array(mesh);
|
|
}
|
|
|
|
let blob = (center, radius, x, y, z) => {
|
|
x -= center[0];
|
|
y -= center[1];
|
|
z -= center[2];
|
|
return Math.max(0, 1 - .16 * (x*x + y*y + z*z) / (radius * radius));
|
|
}
|
|
|
|
var implicitFunction = (x,y,z,args) => {
|
|
let ret = -.5;
|
|
let x4 = _x => _x*_x*_x*_x;
|
|
args.forEach((paras, _)=>{
|
|
const center = paras[0],
|
|
radius = paras[1];
|
|
ret += x4(blob(center, radius, x, y, z));
|
|
});
|
|
return ret;
|
|
}
|
|
|