Deforming a random Delaunay triangulation into its Tutte embedding

tutte.js

		 var scene = new THREE.Scene();
		 var camera = new THREE.PerspectiveCamera(45, window.innerWidth / window.innerHeight, 0.1, 1000);
		 camera.position.z = 300;
		 var renderer = new THREE.WebGLRenderer();
		 renderer.setSize(window.innerWidth, window.innerHeight);
		 renderer.setClearColor( 0xffffff, 1 );
		 document.body.appendChild(renderer.domElement);

		 var material = new THREE.MeshBasicMaterial({
		     color: "white",
		     wireframe: true
		 });

		 var size = { x: 200, z: 200 };
		 var pointsCount = 30;
		 var points3d = [];
		 var vertices = new Float32Array(3 * pointsCount);
		 for (let i = 0; i < pointsCount; i++) {
		     let x = THREE.Math.randFloatSpread(size.x);
		     let y = THREE.Math.randFloatSpread(size.z);
		     
		     points3d.push(new THREE.Vector3(x, y, 0));
		     vertices[3 * i] = x;
		     vertices[3 * i + 1] = y;
		     vertices[3 * i + 2] =  0;
		 }
		 
		 var indexDelaunay = Delaunator.from(
		     points3d.map(v => {
			 return [v.x, v.y];
		     })
		 );

		 var indices = new Uint16Array(3 * indexDelaunay.triangles.length);
		 for (let i = 0; i < indexDelaunay.triangles.length; i++) indices[i] = indexDelaunay.triangles[i];

		 var v = pointsCount;
		 var _L = math.zeros(2 * v, 2 * v, "sparse");

		 var edges = {};
		 for (var i = 0; i < indexDelaunay.triangles.length / 3; i++) {
		     var a = indices[3 * i];
		     var b = indices[3 * i + 1];
		     if (a > b) {a += b; b = a - b; a -= b;}
		     if (isNaN(edges[[a, b]])) edges[[a, b]]= 0;
		     edges[[a, b]] += 1;

		     a = indices[3 * i + 1];
		     b = indices[3 * i + 2];
		     if (a > b) {a += b; b = a - b; a -= b;}
		     if (isNaN(edges[[a, b]])) edges[[a, b]]= 0;
		     edges[[a, b]] += 1;

		     a = indices[3 * i + 2];
		     b = indices[3 * i];
		     if (a > b) {a += b; b = a - b; a -= b;}
		     if (isNaN(edges[[a, b]])) edges[[a, b]]= 0;
		     edges[[a, b]] += 1;
		 }

		 var interior = new Array(v);
		 interior.fill(true);
		 for (var ed in edges) {
		     var str = ed.toString().split(",");
		     var a = parseInt(str[0]);
		     var b = parseInt(str[1]);
		     if (edges[ed] == 1) {interior[a] = false; interior[b] = false;}
		 }
		 
		 for (var j = 0; j < 2 * v; j++) _L.set([j, j], 0);
		 for (var j = 0; j < indexDelaunay.triangles.length / 3; j++)
		 {
		     var a = indexDelaunay.triangles[3 * j];
		     var b = indexDelaunay.triangles[3 * j + 1];
		     var c = indexDelaunay.triangles[3 * j + 2];

		     if (interior[a])
		     {
			 _L.set([2 * a, 2 * b], -1);
			 _L.set([2 * a, 2 * c], -1);
			 _L.set([2 * a, 2 * a], _L.get([2 * a, 2 * a]) + 1);
			 
			 _L.set([2 * a + 1, 2 * b + 1], -1);
			 _L.set([2 * a + 1, 2 * c + 1], -1);
			 _L.set([2 * a + 1, 2 * a + 1], _L.get([2 * a + 1, 2 * a + 1]) + 1);
		     }
		     else
		     {
			 _L.set([2 * a, 2 * a], 1);
			 _L.set([2 * a + 1, 2 * a + 1], 1);
		     }
		     
		     if (interior[b])
		     {
			 
			 _L.set([2 * b, 2 * a], -1);
			 _L.set([2 * b, 2 * c], -1);
			 _L.set([2 * b, 2 * b], _L.get([2 * b, 2 * b]) + 1);
			 
			 _L.set([2 * b + 1, 2 * a + 1], -1);
			 _L.set([2 * b + 1, 2 * c + 1], -1);
			 _L.set([2 * b + 1, 2 * b + 1], _L.get([2 * b + 1, 2 * b + 1]) + 1);
		     }
		     else
		     {
			 _L.set([2 * b, 2 * b], 1);
			 _L.set([2 * b + 1, 2 * b + 1], 1);
		     }
		     
		     if (interior[c])
		     {
			 _L.set([2 * c, 2 * a], -1);
			 _L.set([2 * c, 2 * b], -1);
			 _L.set([2 * c, 2 * c], _L.get([2 * c, 2 * c]) + 1);
			 
			 _L.set([2 * c + 1, 2 * a + 1], -1);
			 _L.set([2 * c + 1, 2 * b + 1], -1);
			 _L.set([2 * c + 1, 2 * c + 1], _L.get([2 * c + 1, 2 * c + 1]) + 1);		    
		     }
		     else
		     {
			 _L.set([2 * c, 2 * c], 1);
			 _L.set([2 * c + 1, 2 * c + 1], 1);
		     }
		 }

		 _a = new Array(2 * v);
		 _a.fill(0);
		 for (var i = 0; i < v; i++)
		 {
		     if (!interior[i])
		     {
			 _a[2 * i] =  vertices[3 * i];
			 _a[2 * i + 1] = vertices[3 * i + 1];
		     }
		 }
		 
		 _x = math.lusolve(_L, _a);
		 var t = 0;
		 var render = function() {
		     while(scene.children.length > 0){ 
			 scene.remove(scene.children[0]); 
		     }

		     t += .11;
		     var p = (Math.sin(t) + 1) /2;
		     
		     var verticesDef = new Float32Array(3 * pointsCount);
		     
		     for (var j  = 0; j < v; j++) {
			 verticesDef[3 * j] = p * vertices[3 * j] + (1 - p ) * _x.get([2 * j, 0]);
			 verticesDef[3 * j + 1] = p * vertices[3 * j + 1] + (1 - p) * _x.get([2 * j + 1, 0]);
		     }

		     var geomDef = new THREE.BufferGeometry();
		     geomDef.addAttribute('position', new THREE.BufferAttribute( verticesDef, 3 ) );
		     geomDef.addAttribute( 'index', new THREE.BufferAttribute(indices, 3 ) );
		     var cloudDef = new THREE.Points(
			 geomDef,
			 new THREE.PointsMaterial({ color: "red", size: 10 })
		     );
		     scene.add(cloudDef);
		     
		     var plDef = new THREE.Mesh(
			 geomDef,
			 new THREE.MeshLambertMaterial({color: "red", wireframe: true})
		     );
		     scene.add(plDef);
		     
		     requestAnimationFrame(render);
		     renderer.render(scene, camera);
		 };
		 render();

requires math.js, three.js, and delaunator.