Surface meshing in the browser


precision mediump float;
varying vec3 vTransformedNormal;
varying vec4 vPosition;
uniform vec3 uMaterialDiffuseColor;
void main(void) {
vec3 materialDiffuseColor = uMaterialDiffuseColor;
vec3 ambientLightWeighting = vec3(.2,.2,.2);
vec3 directionalColor = vec3(.8,.8,.8);
vec3 normal = normalize(vTransformedNormal);
vec3 uLightingDirection = vec3(-0.57,-0.57,-0.57);
float diffuseLightBrightness = abs(dot(normal,uLightingDirection));
vec3 diffuseLightWeighting = diffuseLightBrightness*directionalColor;
vec3 r = -reflect(uLightingDirection, normal);
r = normalize(r);
vec3 v = -vPosition.xyz;
v = normalize(v);
vec4 spec = vec4(.5,.5,.5,1.0)*pow(max(0.0,dot(r, v)), 16.0);
gl_FragColor = vec4(ambientLightWeighting * materialDiffuseColor
+ materialDiffuseColor * diffuseLightWeighting,
1.0)+spec;
// gl_FragColor = vec4(0,0,0,1.0);
}

attribute vec3 aVertexPosition;
attribute vec3 aVertexNormal;
uniform mat4 uMVMatrix;
uniform mat4 uPMatrix;
uniform mat3 uNMatrix;
uniform vec3 uAmbientLightingColor;
uniform vec3 uDirectionalDiffuseColor;
varying vec3 vTransformedNormal;
varying vec4 vPosition;
void main(void) {
gl_PointSize = 3.0;
gl_Position = uPMatrix * uMVMatrix * vec4(aVertexPosition, 1.0);
vPosition = uMVMatrix * vec4(aVertexPosition, 1.0);
vTransformedNormal =  aVertexNormal;
}

var canvas;
var mouseX,mouseY,pmouseX,pmouseY,isMouseDown,mouseDragging;
var ctx;
var paused = false;
function init() {
canvas = document.getElementById("test");
depth = canvas.width/8.0;
//ctx = canvas.getContext("2d");
initGL(canvas);
initBuffers();
NScamera.distance = 500;
vec3.set(NScamera.center, canvas.width/2, canvas.height/2,20);
initPts();
initialize();
step();
gl.enable(gl.DEPTH_TEST);
}
var ptGrid;
var density = 10;
var depth = 0;
function initPts() {
pts = new Array();
var weightsX = [Math.random()*80+20,Math.random()*70+15,Math.random()*65+5,Math.random()*60+5,Math.random()*45+5,Math.random()*25+5];
var weightsY = [Math.random()*80+20,Math.random()*70+15,Math.random()*65+5,Math.random()*60+5,Math.random()*45+5,Math.random()*25+5];
var totalWeight = 0;
var len = weightsX.length;
for(var i=0;i< len;++i) {
totalWeight += weightsX[i] + weightsY[i];
}
vec3.set(NScamera.center, canvas.width/2, canvas.height/2,totalWeight/2.0);
var gw = Math.ceil(canvas.width/density);
var gh = Math.ceil(canvas.height/density);
var gd = Math.ceil(totalWeight/density);
var gwh = gw*gh;
ptGrid = new Array(gw*gh*gd);
for(var i=0;i< ptGrid.length ;++i) ptGrid[i] = new Array();
var tries = 0;
var x = 0, y = 0;
var p;
var dx = vec3.create();
var dy = vec3.create();
while(tries < 2000) {
x = Math.random()*canvas.width;
y = Math.random()*canvas.width;
if((x-canvas.width/2.0)*(x-canvas.width/2.0)+(y-canvas.height/2.0)*(y-canvas.height/2.0) < 250*250) {
z = 0;
for(var i=0;i< len ;++i) {
z += (Math.cos(x/canvas.width*Math.PI*(2+i*1.5))+1)*0.5*weightsX[i];
z += (Math.sin(y/canvas.height*Math.PI*(2+i*1.5))+1)*0.5*weightsY[i];
}
var xi=Math.floor(x/density);
var yi=Math.floor(y/density);
var zi=Math.floor(z/density);
var hit = false;
for(var i=Math.max(0,xi-1);i<=Math.min(gw-1,xi+1);++i) {
for(var j=Math.max(0,yi-1);j<=Math.min(gh-1,yi+1);++j) {
for(var k=Math.max(0,zi-1);k<=Math.min(gd-1,zi+1);++k) {
for(var q=0;q < ptGrid[k*gwh+j*gw+i].length;++q) {
p = ptGrid[k*gwh+j*gw+i][q];
if((p.pt[0]-x)*(p.pt[0]-x)+(p.pt[1]-y)*(p.pt[1]-y) < density*density) hit = true;
}
}
}
}
if(!hit) {
tries = 0;
var pt = new Vertex();
vec3.set(pt.pt,x,y,z);//(x*x+y*y)/10000.0
ptPos[pts.length*3] = pt.pt[0];
ptPos[pts.length*3+1] = pt.pt[1];
ptPos[pts.length*3+2] = pt.pt[2];
ptPosBuff.numItems++;
pts.push(pt);
ptGrid[zi*gwh+yi*gw+xi].push(pt);
}
tries++;
}
}
gl.bindBuffer(gl.ARRAY_BUFFER, ptPosBuff);
gl.bufferData(gl.ARRAY_BUFFER, ptPos, gl.STATIC_DRAW);
}
var pivotEdge = null;
function keyPress() {
reset();
}
function reset() {
ptPosBuff.numItems = 0;
initPts();
initialize();
vertexPosBuff.numItems = 0;
}
var normalMatrix = mat3.create();
function draw() {
NScamera.step();
NScamera.velocityY = (.005+NScamera.velocityY*3)/4.0;
NScamera.velocityX = (.003+NScamera.velocityX*3)/4.0;
gl.viewport(0, 0, gl.viewportWidth, gl.viewportHeight);
gl.clear(gl.COLOR_BUFFER_BIT | gl.DEPTH_BUFFER_BIT);
mat4.perspective(pMatrix,45, gl.viewportWidth / gl.viewportHeight, 0.1, 10000.0);
mat4.identity(mvMatrix);
NScamera.feed(mvMatrix);
mat4.toInverseMat3(mvMatrix,normalMatrix);
mat3.transpose(normalMatrix,normalMatrix);
gl.bindBuffer(gl.ARRAY_BUFFER, ptPosBuff);
gl.vertexAttribPointer(shaderProgram.vertexPositionAttribute, ptPosBuff.itemSize, gl.FLOAT, false, 0, 0);
gl.bindBuffer(gl.ARRAY_BUFFER, vertexNormBuff);
gl.vertexAttribPointer(shaderProgram.vertexNormalAttribute, vertexNormBuff.itemSize, gl.FLOAT, false, 0, 0);
gl.drawArrays(gl.POINTS,0,ptPosBuff.numItems);
var norm = vec3.create();
for(var i=vertexPosBuff.numItems/3;i < triangles.length;++i) {
planeNormal(norm,triangles[i][0].pt,triangles[i][1].pt,triangles[i][2].pt);
vec3.normalize(norm,norm);
for(var j=0;j < 3;++j) {
vertexPos[vertexPosBuff.numItems*3] = triangles[i][j].pt[0];
vertexPos[vertexPosBuff.numItems*3+1] = triangles[i][j].pt[1];
vertexPos[vertexPosBuff.numItems*3+2] = triangles[i][j].pt[2];
vertexNorm[vertexPosBuff.numItems*3] = norm[0];
vertexNorm[vertexPosBuff.numItems*3+1] = norm[1];
vertexNorm[vertexPosBuff.numItems*3+2] = norm[2];
vertexPosBuff.numItems++;
}
}
gl.bindBuffer(gl.ARRAY_BUFFER, vertexPosBuff);
gl.bufferData(gl.ARRAY_BUFFER, vertexPos, gl.DYNAMIC_DRAW);
gl.vertexAttribPointer(shaderProgram.vertexPositionAttribute, vertexPosBuff.itemSize, gl.FLOAT, false, 0, 0);
gl.bindBuffer(gl.ARRAY_BUFFER, vertexNormBuff);
gl.bufferData(gl.ARRAY_BUFFER, vertexNorm, gl.DYNAMIC_DRAW);
gl.vertexAttribPointer(shaderProgram.vertexNormalAttribute, vertexNormBuff.itemSize, gl.FLOAT, false, 0, 0);
gl.drawArrays(gl.TRIANGLES,0,vertexPosBuff.numItems);
}
(function() {
var requestAnimationFrame = window.requestAnimationFrame || window.mozRequestAnimationFrame || window.webkitRequestAnimationFrame;
window.requestAnimationFrame = requestAnimationFrame;
}) ();
function step() {
window.requestAnimationFrame(step);
//just one boundary for now
if(!paused) {
for(var i=0;i < boundaries.length;++i) {
var bound = boundaries[i].e;
var curr = bound.next;
var dumb = curr.active == 0;
if(curr == bound) dumb = false;
while(dumb ) {
curr = curr.next;
dumb = curr.active == 0;
if(curr == bound) dumb = false;
}
if(curr.active == 1) {
pivot(curr);
pivotEdge = curr;
}
}
}
draw();
}
var vertexPos;
var vertexPosLen = 0;
var vertexNorm;
var ptPos;
var ptPosBuff;
var vertexPosBuff;
var vertexNormBuff;
function initBuffers() {
ptPos = new Float32Array(3*10000);
ptPosBuff = gl.createBuffer();
vertexPos = new Float32Array(3*3*10000);
vertexPosBuff = gl.createBuffer();
vertexNorm = new Float32Array(3*3*10000);
vertexNormBuff = gl.createBuffer();
gl.bindBuffer(gl.ARRAY_BUFFER, vertexPosBuff);
gl.bufferData(gl.ARRAY_BUFFER, vertexPos, gl.DYNAMIC_DRAW);
vertexPosBuff.itemSize = 3;
vertexPosBuff.numItems = 0;
gl.bindBuffer(gl.ARRAY_BUFFER, vertexNormBuff);
gl.bufferData(gl.ARRAY_BUFFER, vertexNorm, gl.DYNAMIC_DRAW);
vertexNormBuff.itemSize = 3;
vertexNormBuff.numItems = 0;
gl.bindBuffer(gl.ARRAY_BUFFER, ptPosBuff);
gl.bufferData(gl.ARRAY_BUFFER, ptPos, gl.DYNAMIC_DRAW);
ptPosBuff.itemSize = 3;
ptPosBuff.numItems = 0;
}
}
}
function mouseDraggedMan(event) {
NScamera.mouseDragged(mouseX-pmouseX,-(mouseY-pmouseY), mouseX, mouseY,mouseButton);
}
function mouseUpMan(event) {
isMouseDown = false;
}
function mouseMoveMan(event) {
pmouseX = mouseX;
pmouseY = mouseY;
mouseX = event.clientX;
mouseY = event.clientY;
if(isMouseDown) {
mouseDragging = true;
mouseDraggedMan(event);
}
}
function mouseDownMan(event) {
isMouseDown = true;
mouseDragging = false;
mouseButton = event.which;
}
mat4.toInverseMat3=function(a,b){var c=a[0],d=a[1],e=a[2],g=a[4],f=a[5],h=a[6],i=a[8],j=a[9],k=a[10],l=k*f-h*j,o=-k*g+h*i,m=j*g-f*i,n=c*l+d*o+e*m;if(!n)return null;n=1/n;b||(b=mat3.create());b[0]=l*n;b[1]=(-k*d+e*j)*n;b[2]=(h*d-e*f)*n;b[3]=o*n;b[4]=(k*c-e*i)*n;b[5]=(-h*c+e*g)*n;b[6]=m*n;b[7]=(-j*c+d*i)*n;b[8]=(f*c-d*g)*n;return b};

init();



A common problem we come across is needing to make meshes of complex surfaces. This can be done with many different approaches, but regardless of the method, it is often quite difficult. Especially when doing simulations on surfaces, you need a well-behaved triangle mesh, that is one with a controllable density of elements and nicely shaped triangles. This can be done by generating a set of points on that surface and then triangulating them. But triangulating surfaces embedded in 3-space is no simple task. Recently, I came across an interesting approach to this problem which works especially well with generated point clouds called the Ball-Pivot Algorithm. I implemented it in javascript (you should see it running above).

Let’s start with the basic idea of triangulation in 2D. A triangulation of a set of points is just a set of triangles that covers all those points and don’t overlap.

Now not all triangulations are created equal. We want to make triangulations whose triangles are nicely shaped, as close to equilateral as possible. This leads to better numerical stability and convergence in simulations.

This leads to something called the Delaunay triangulation. This is a really neat triangulation defined such that circumcircle of any triangle in the triangulation does not contain any other points. The circumcircle is the smallest circle enclosing the triangle. One thing that’s not entirely obvious is that this triangulation always exists. And the thing that is great about this triangulation is that it maximizes the minimum angle of the triangles. This means it gives us the best shaped triangles we can get with the points we have.

The problem comes when we make the jump from 2D to 3D. We can generalize the idea of Delaunay triangulations to 3D, giving us tetrahedron instead of triangles and circumspheres instead of circumcircles. However, what we want is surface of triangles embedded in 3-space, not a 3 dimensional solid. One way to get around this is called alpha shapes. This is a subset of the 3D Delaunay triangulation but only including faces of a certain size. Alpha shapes can create good surface triangulations; however, it is necessary to compute the full 3D Delaunay, which is complex in computation and implementation as well as numerically unstable.

So in comes the Ball-Pivot algorithm. This algorithm is nice because it works in linear time for points of a constant density, and it is exceptionally cute conceptually. The idea is you have a ball of a certain radius sitting in a triangle. This ball rolls over the edge of the triangle and keeps rolling until it hits another point. Then, that edge and the point it hits make a new triangle which is added to the triangulation, and the ball rolls over a new edge. If you start with a good triangle, such that the ball contains no other points, then the new triangle the ball rolls onto won’t contain any other points either because any point would have been hit before hand. This maintains the Delaunay property! In fact, the ball pivot algorithm produces a subset of an alpha shape.

There are a few caveats. If the curvature and density of the surface points is too high, a larger ball won’t be able to fit in the necessary crevices, missing some points. If the ball is too small, it can fall through large triangles and leave holes. While these issues can worked around, the ball pivot algorithm works best for surfaces with a consistent density of points.

Feel free to look at the code inside, but it is not optimized and not pretty. This was thrown together in 2-3 days, so it definitely lacks polish.