Classes from an old processing library project for modelling manifold n-gon meshes.

Edge.java

/**
 * 
 */
package manifold;

import processing.core.*;

/**
 * a class for edges
 */
public class Edge extends Topology {

	PApplet myParent;
	Vertex start;
	Vertex end;
	Face left;
	Face right;

	/**
	 * constructor for Edge class
	 * @param theParent
	 * 			the parent PApplet
	 * @param start
	 * 			vertex at which to start the edge
	 * @param end
	 * 			vertex at which to end the edge
	 * @param left
	 * 			face on the left of the edge
	 */
	Edge(PApplet theParent, Vertex start, Vertex end, Face left) {
		myParent = theParent;
		this.start = start;
		this.end = end;
		this.left = left;
		this.right = null;
	}

	/**
	 * draw edge (boundary edges drawn in thicker stroke)
	 */
	void draw() {
		myParent.stroke(0); // TODO: set colour here
		if (this.boundary() == false) {
			myParent.strokeWeight(2);
		}
		else {
			myParent.strokeWeight(5);
		}
		myParent.line(this.start.position.x, this.start.position.y, this.start.position.z,
				this.end.position.x, this.end.position.y, this.end.position.z);
	}

	/**
	 * returns the length of the edge
	 * @return float
	 */
	float length() {
		return this.start.position.dist(this.end.position);
	}

	/**
	 * returns a vector in direction of edge, with magnitude of edge's length
	 * @return PVector
	 */
	PVector vector() {
		return PVector.sub(this.end.position, this.start.position);
	}

	/**
	 * as vector() but starting at specified vertex
	 * @param v
	 * 			the vertex at which the vector should begin
	 * @return PVector
	 */
	PVector vectorFrom(Vertex v) { 
		if (this.start == v) return this.vector();
		else if (this.end == v) return PVector.mult(this.vector(), -1);
		else return null;
	}

	/**
	 * returns a normalised vector starting at a specified vertex and in the direction of the other vertex
	 * @param v
	 * 			the vertex at which the vector should begin
	 * @return PVector
	 */
	PVector unitVectorFrom(Vertex v) {
		return this.vectorFrom(v).normalize(null);
	}

	/**
	 * returns true if the edge is a boundary edge
	 * @return boolean
	 */
	boolean boundary() {
		return (this.right == null);
	}

	/**
	 * returns a PVector with the coordinates of the mid-point of the edge
	 * @return PVector
	 */
	PVector midPoint() {
		return PVector.div(PVector.add(this.start.position, this.end.position), 2);
	}

	/**
	 * reverses the direction of the edge
	 */
	void reverse() {
		Vertex tempVertex = this.start;
		this.start = this.end;
		this.end = tempVertex;
		Face tempFace = this.left;
		this.left = this.right;
		this.right = tempFace;
	}

	/**
	 * returns the vertex opposite the specified vertex (or null if the vertex isn't on the edge) 
	 * @param v
	 * 			the known vertex
	 * @return Vertex
	 */
	Vertex otherVertex(Vertex v) {
		if (this.start == v) return this.end;
		else if (this.end == v) return this.start;
		else return null;
	}
}

Face.java

/**
 * 
 */
package manifold;

import processing.core.*;

import java.util.ArrayList;
import java.util.List;
import java.math.*;

/**
 * @author Will
 *
 */
public class Face extends Topology implements PConstants{

	PApplet myParent;
	Vertex[] vertices;
	//Edge[] edges;

	Face(PApplet theParent, Vertex[] vertices) {
		this.myParent = theParent;
		this.vertices = vertices;
		//this.edges = new Edge[vertices.length];
	}

	private void draw(boolean smooth) {
		myParent.noStroke();
		myParent.fill(200, 200);
		if (this.vertices.length <= 3 || smooth) {
			// draw face normally for triangles
			// or with vertex normals for a smooth surface
			myParent.beginShape();
			for (Vertex v : this.vertices) {
				if (smooth) {
					PVector n = v.normal();
					myParent.normal(n.x, n.y, n.z);
				}
				//myParent.fill(v.color());
				myParent.vertex(v.position.x, v.position.y, v.position.z);
			}
			myParent.endShape(CLOSE);
		}
		else {
			// draw 'fan' around centroid of face
			PVector c = this.centroid();
			for (int i = 0; i < this.vertices.length; i++) {
				PVector a = this.vertices[i].position;
				PVector b = this.vertices[(i+1)%this.vertices.length].position;
				myParent.beginShape();
				//myParent.fill(this.vertices[i].color());
				myParent.vertex(a.x, a.y, a.z);
				//myParent.fill(this.vertices[(i+1)%this.vertices.length].color());
				myParent.vertex(b.x, b.y, b.z);
				myParent.vertex(c.x, c.y, c.z);
				myParent.endShape(CLOSE);
			}
		}
	}

	void drawFaceNormals() { 
		// draw normals
		PVector a = this.centroid();
		PVector b = PVector.add(this.centroid(), PVector.mult(this.normal(), 1));
		myParent.strokeWeight(1);
		myParent.stroke(0);
		myParent.line(a.x, a.y, a.z, b.x, b.y, b.z);
	}
	
	void draw(boolean normals, boolean smooth) {
		this.draw(smooth);
		if (normals) {
			this.drawFaceNormals();
		}
	}

	PVector centroid() {
		PVector c = new PVector(0, 0, 0);
		for (Vertex v : this.vertices) {
			c.add(v.position);
		}
		c.div(this.vertices.length);
		return c;
	}

	PVector normal() {
		PVector n = new PVector();
		for (int i=0; i<this.vertices.length; i++) {
			n.add(this.vertices[i].position.cross(this.vertices[(i+1)%this.vertices.length].position));
		}
		return n.normalize(null);
	}

	float area() {
		float area;
		if (this.vertices.length == 3) {
			// triangle
			PVector A = this.vertices[1].vectorFrom(this.vertices[0]);
			PVector B = this.vertices[1].vectorFrom(this.vertices[2]);
			area = 0.5f * A.cross(B).mag(); 
		}
		else {
			// planar polygon
			// TODO: use fan method instead
			float A = 0;
			for (int i = 0; i < this.vertices.length - 1; i++) {
				A += (this.vertices[i].position.x * this.vertices[i+1].position.y) -
						(this.vertices[i+1].position.x * this.vertices[i].position.y);
			}
			area = myParent.abs((float) (A * 0.5));
		}
		return area;
	}
	
	public PVector circumcenter() {
		// assume triangle
		PVector P1, P2, P3, A, B, C;
		P1 = this.vertices[0].position;
		P2 = this.vertices[1].position;
		P3 = this.vertices[2].position;
		A = PVector.sub(P2, P1);
		B = PVector.sub(P3, P2);
		C = PVector.sub(P1, P3);
		float K = this.area();
		return PVector.add(PVector.mult(PVector.add(P1, P3), 0.5f), 
				PVector.mult((C.cross(A.cross(B))), (float) ((A.dot(B)) / (8 * Math.pow(K, 2)))));
	}
	
	public float circumradius() {
		// assume triangle
		return this.circumcenter().dist(this.vertices[0].position);
	}

	Edge[] edges() {
		// For each pair of vertices, get the edge!
		//logger(this, "DEBUG", "edges; Finding edges for face " + this);
		int n = this.vertices.length;
		Edge[] edges = new Edge[n];
		for (int i = 0; i < n; i++) {
			edges[i] = this.vertices[i].edgeWith(this.vertices[(i+1)%n]);
		} 
		return edges;
	}
}

Manifold.java

/**
 * Manifold
 * A collection of topological classes for constructing and manipulating meshes.
 * http://pearswj.co.uk/manifold
 *
 * Copyright © 2012 Will Pearson http://pearswj.co.uk
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 * 
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General
 * Public License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place, Suite 330,
 * Boston, MA  02111-1307  USA
 * 
 * @author      Will Pearson http://pearswj.co.uk
 * @modified    03/25/2013
 * @version     0.1.1 (1)
 */

package manifold;


import java.io.File;
import java.io.PrintWriter;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

import processing.core.PApplet;
import processing.core.PConstants;
import processing.core.PVector;
//import java.util.logging.*;
//import org.apache.log4j.*;

/**
 * a class for manifolds
 */

public class Manifold implements PConstants {

	// myParent is a reference to the parent sketch
	PApplet myParent;

	private List<Vertex> vertices;
	private List<Edge> edges;
	private List<Face> faces;

	public final static String VERSION = "0.1.1";
	
	//private static final Logger logger = Logger.getLogger(Manifold.class.getName());
	
	// manage the displaying of elements
	private boolean showVertices = false;
	private boolean showVertexNormals = false;
	private boolean showCurvatureTensors = false;
	private boolean showVertexLabels = false;
	private boolean showEdges = true;
	private boolean showFaces = false;
	private boolean showFaceNormals = false;
	private boolean showFaceLabels = false;
	private boolean smoothFaces = false;
	
	private float[] modelRotations = new float[]{0, 0, 0};

	/**
	 * a Constructor, usually called in the setup() method in your sketch to
	 * initialize and start the library.
	 * 
	 * @example Hello
	 * @param theParent
	 */
	public Manifold(PApplet theParent) {
		myParent = theParent;
		welcome();
		this.vertices = new ArrayList<Vertex>(); 
		this.edges = new ArrayList<Edge>();
		this.faces = new ArrayList<Face>();
		//theParent.registerDraw(this);
	}


	private void welcome() {
		System.out.println("Manifold 0.1.1 by Will Pearson");
	}


	/**
	 * return the version of the library.
	 * 
	 * @return String
	 */
	public static String version() {
		return VERSION;
	}

	//---------------------------------------------------------//
	//              Add vertices, faces and edges              //
	//---------------------------------------------------------//

	/**
	 * Add a new vertex
	 * 
	 * @param v
	 * 			the vertex being added to the manifold
	 * @return Vertex
	 */
	public Vertex addVertex(Vertex v) {
		this.vertices.add(v);
		return v;
	}

	/**
	 * add a new vertex
	 * 
	 * @param position
	 * 			the position of a new vertex
	 * @return Vertex
	 */
	public Vertex addVertex(PVector position) {
		return this.addVertex(new Vertex(myParent, position));
	}

	// Add a new face:

	/**
	 * add a new face by its vertices
	 * 
	 * @param vertices
	 * 			an array of the vertices that make the face (supplied in an anti-clockwise order)
	 * @return Face
	 */
	public Face addFace(Vertex[] vertices) {
		Face f = new Face(myParent, vertices);
		// Find/create edges.
		for (int i = 0; i < vertices.length; i++) {
			Vertex start = vertices[i];
			Vertex end = vertices[(i+1)%vertices.length];
			Edge e = end.edgeTo(start); // find an edge that exists in the opposite dir
			if (e != null) {
				e.right = f;
				//f.edges[i] = e; // Add edge to face.
			} 
			else if (start.edgeTo(end) == null) {
				// check if an edge already exists in the same direction and if not
				// create an edge.
				this.addEdge(start, end, f);
				//f.edges[i] = this.addEdge(start, end, f);
			} 
			else {
				// don't add the new face if there's an existing edge in the same dir.
				// (something's probably wrong in the vertex order...)
				return null;
			}
			// Note: should check that there isn't already an edge in this direction
			// (i.e. face created wrong...)
		}   
		this.faces.add(f);
		return f;
	}

	/**
	 * add a new face by the indices of its vertices
	 * 
	 * @param vIndex
	 * 			an array of the indices of the vertices that make the face (supplied in an anti-clockwise order)
	 * @return Face
	 */
	public Face addFace(int[] vIndex) {
		// Add a new face, described by vertex indices
		// TODO: catch 'index out of bounds' errors
		Vertex[] vertices = new Vertex[vIndex.length];
		for (int i = 0; i < vIndex.length; i++) {
			vertices[i] = this.vertices.get(vIndex[i]);
		}
		return this.addFace(vertices);
	}


	// Add a new edge:

	private Edge addEdge(Edge e) {
		this.edges.add(e);
		if (!e.start.edges.contains(e)) {
			e.start.edges.add(e);
		}
		if (!e.end.edges.contains(e)) {
			e.end.edges.add(e);
		}
		return e;
	}

	private Edge addEdge(Vertex start, Vertex end, Face left) {
		return this.addEdge(new Edge(myParent, start, end, left));
	}

	//---------------------------------------------------------//
	//            Remove vertices, faces and edges             //
	//---------------------------------------------------------//

	/**
	 * remove a face
	 * 
	 * @param f
	 * 			the face object being removed
	 * @return boolean
	 */
	public boolean removeFace(Face f) {
		// remove a face
		//logger(this, "INFO", "removeFace; removing face " + f + "...");
		for (Edge e : f.edges()) {
			//logger(this, "DEBUG", "removeFace; removing face " + f + " from edge " + e);
			if (e.left == f) {
				// f is on the left of the edge
				if (e.right != null) {
					// the edge has a right face, reverse the edge
					e.reverse();
					e.right = null;
				}
				else {
					// no face on right, just remove the edge
					//this.removeEdge(e);
					e.start.edges.remove(e);
					e.end.edges.remove(e);
					this.edges.remove(e);
				}
			}
			else {
				// f is on the left of the edge
				e.right = null;
			}
		}
		return this.faces.remove(f);
	}

	private boolean removeEdge(Edge e) {
		// remove an edge
		// remove the edge from the vertices that list it (i.e. start and end)
		e.start.edges.remove(e);
		e.end.edges.remove(e);
		// remove the faces that depended on this edge
		this.removeFace(e.left);
		this.removeFace(e.right);
		return this.edges.remove(e);
	} 

	//---------------------------------------------------------//
	//                      Draw Methods                       //
	//---------------------------------------------------------//
	
	public void draw() {
		if (showVertices) {
		    this.drawVertices(showVertexNormals, showCurvatureTensors, showVertexLabels);
		  }
		  if (showEdges) {
		    this.drawEdges();
		  }
		  if (showFaces) {
		    this.drawFaces(smoothFaces, showFaceNormals, showFaceLabels);
		  }
	}

	private void drawVertices(boolean normals, boolean tensors, boolean labels) {
		for (Vertex v : this.vertices) {
			v.draw(normals, tensors);
			if (labels) this.drawLabel(""+this.vertices.indexOf(v), PVector.add(v.position, v.normal()));
		}
	}

	private void drawVertices() {
		// quick draw
		for (Vertex v : this.vertices) {
			v.draw();
		}
	}

	private void drawEdges() {
		for (Edge e : this.edges) {
			e.draw();
		}
	}

	private void drawFaces(boolean smooth, boolean normals, boolean labels) {
		for (Face f : this.faces) {
			f.draw(normals, smooth);
			if (labels) {
				String label = ""+this.faces.indexOf(f);
                PVector textPosition = PVector.add(f.centroid(), f.normal());
                this.drawLabel(label, textPosition);
			}
		}
	}

	private void drawFaces() {
		// quick draw
		for (Face f : this.faces) {
			f.draw(false, false);
		}
	}
	
	private void drawLabel(String label, PVector textPosition) {
        myParent.textSize(0.1f);
        myParent.fill(0);
        //float[] rota = cam.getRotations();
        myParent.pushMatrix();
        myParent.translate(textPosition.x, textPosition.y, textPosition.z);
        myParent.rotateX(this.modelRotations[0]);
        myParent.rotateY(this.modelRotations[1]);
        myParent.rotateZ(this.modelRotations[2]);
        myParent.text(label, 0, 0, 0);
        myParent.popMatrix();
	}
	
	public void showVertices(boolean showVertices) {
		this.showVertices = showVertices;
	}


	public void showVertexNormals(boolean showVertexNormals) {
		this.showVertexNormals = showVertexNormals;
	}
	
	
	public void showCurvatureTensors(boolean showCurvatureTensors) {
		this.showCurvatureTensors = showCurvatureTensors;
	}


	public void showVertexLabels(boolean showVertexLabels) {
		this.showVertexLabels = showVertexLabels;
	}


	public void showEdges(boolean showEdges) {
		this.showEdges = showEdges;
	}


	public void showFaces(boolean showFaces) {
		this.showFaces = showFaces;
	}


	public void showFaceNormals(boolean showFaceNormals) {
		this.showFaceNormals = showFaceNormals;
	}


	public void showFaceLabels(boolean showFaceLabels) {
		this.showFaceLabels = showFaceLabels;
	}
	
	public void smoothFaces(boolean smoothFaces) {
		this.smoothFaces = smoothFaces;
	}
	
	public void setModelRotations(float rotateX, float rotateY, float rotateZ) {
        this.modelRotations = new float[]{rotateX, rotateY, rotateZ};
	}

	//---------------------------------------------------------//
	//                    Accessor Methods                     //
	//---------------------------------------------------------//

	// return arrays of vertices, edges and faces

	public Vertex[] vertices() {
		return this.vertices.toArray(new Vertex[this.vertices.size()]);
	}

	public Edge[] edges() {
		return this.edges.toArray(new Edge[this.edges.size()]);
	}

	public Face[] faces() {
		return this.faces.toArray(new Face[this.faces.size()]);
	}

	//---------------------------------------------------------//
	//                    Conway Operations                    //
	//---------------------------------------------------------//

	// Dual

	public Manifold dual() {
		// Find the dual of the current manifold.
		Manifold d = new Manifold(myParent);
		// vertices from faces
		for (Face f : this.faces) {
			d.addVertex(f.centroid());
		}
		// vertices from boundary edges
		Vertex[] boundaryPoints = new Vertex[this.edges.size()];
		//for (Edge e : this.edges) {
		for (int i = 0; i < this.edges.size(); i++) {
			if (this.edges.get(i).boundary() == true) {
				// add boundaryPoint to dual (also store in boundaryPoints for indexing)
				boundaryPoints[i] = d.addVertex(this.edges.get(i).midPoint());
			}
		}
		// faces from vertices (and boundary points)
		for (Vertex v : this.vertices) {
			if (v.edges.size() > 0) { // Check that vertex isn't an orphan 
				//Vertex[] fv = new Vertex[v.edges.size()];
				ArrayList<Vertex> fv = new ArrayList<Vertex>();
				Face[] vf = v.faces();
				//logger(this, "DEBUG", "dual; " + vf.length + " faces on vertex " + v);
				if (v.edges.get(0).boundary() == true) {
					//logger(this, "DEBUG", "dual; adding edge-point for boundary edge");// " + e);
					fv.add(boundaryPoints[this.edges.indexOf(v.edges.get(0))]);
				}
				//for (int i = 0; i < vf.length; i++) {
				for (Face f : vf) {
					//fv[i] = d.vertices.get(this.faces.indexOf(vf[i]));
					//logger(this, "DEBUG", "dual; adding face-point for face");// " + f);
					fv.add(d.vertices.get(this.faces.indexOf(f)));
				}
				if (v.edges.get(v.edges.size()-1).boundary() == true) {
					//logger(this, "DEBUG", "dual; adding edge-point for boundary edge");// " + e);
					fv.add(boundaryPoints[this.edges.indexOf(v.edges.get(v.edges.size()-1))]);
				}
				d.addFace(fv.toArray(new Vertex[0]));
			}
		}
		this.set(d);
		return d;
	}

	//---------------------------------------------------------//
	//                       Subdivision                       //
	//---------------------------------------------------------//

	// Catmull-Clark

	public Manifold catmullClark() {
		// subdivide the manifold following the Catmull-Clark algorithm.
		// see http://rosettacode.org/wiki/Catmull-Clark_subdivision_surface (1)
		// and https://graphics.stanford.edu/wikis/cs148-09-summer/Assignment3Description (2)
		Manifold cc = new Manifold(myParent);

		// for each face, a FACE POINT is created which is the average
		// of all the points of the face.
		Vertex[] facePoints = new Vertex[this.faces.size()];
		for (Face f : this.faces) {
			facePoints[this.faces.indexOf(f)] = cc.addVertex(f.centroid());
		}

		// for each edge, an EDGE POINT is created which is the average
		// between the center of the edge and the center of the segment
		// made with the face points of the two adjacent faces.
		Vertex[] edgePoints = new Vertex[this.edges.size()];
		for (Edge e : this.edges) {
			PVector edgePoint = new PVector();
			if (!e.boundary()) {
				edgePoint.add(e.left.centroid());
				edgePoint.add(e.right.centroid());
				edgePoint.add(e.start.position);
				edgePoint.add(e.end.position);
				edgePoint.div(4);
			}
			else {
				edgePoint = e.midPoint();
			}
			edgePoints[this.edges.indexOf(e)] = cc.addVertex(edgePoint);
		}

		// for each ORIGINAL POINT, update location based on (2)
		Vertex[] origPoints = new Vertex[this.vertices.size()];
		for (Vertex v : this.vertices) {
			PVector newCoords;
			// old coordinates
			PVector oldCoords = v.position.get(); // S
			if (!v.boundary()) { // (-Q + 4E + (n-3)*S)/n
				// average of the face points of the faces the point belongs to
				PVector avgFacePoints = new PVector(); // Q
				for (Face f : v.faces()) {
					avgFacePoints.add(facePoints[this.faces.indexOf(f)].position);
				}
				avgFacePoints.div(v.faces().length);
				// average of the centers of edges the point belongs to
				PVector avgEdgePoints = new PVector(); // E
				for (Edge e : v.edges) {
					avgEdgePoints.add(edgePoints[this.edges.indexOf(e)].position);
				}
				avgEdgePoints.div(v.edges.size());
				// calculate new coordinates
				float n = v.faces().length; // number of faces a point belongs to
				newCoords = PVector.sub(PVector.mult(avgEdgePoints, 4), avgFacePoints);
				newCoords.add(PVector.mult(oldCoords, (n-3)));
				newCoords.div(n);
			}
			else { // S/2 + M/4
				PVector avgBoundaryEdgePoints = new PVector(); // M
				for (Edge e : v.edges()) {
					if (e.boundary()) {
						avgBoundaryEdgePoints.add(e.midPoint());
					}
				}
				newCoords = PVector.div(oldCoords, 2);
				newCoords.add(PVector.div(avgBoundaryEdgePoints, 4));
			}
			origPoints[this.vertices.indexOf(v)] = cc.addVertex(newCoords); // update
		}

		// add faces by linking up each original (moved) point with a face point and
		// the two corresponding edge points
		for (Face f : this.faces) {
			for (int i = 0; i < f.vertices.length; i++) {
				Vertex prev = f.vertices[i];
				Vertex curr = f.vertices[(i+1)%f.vertices.length];
				Vertex next = f.vertices[(i+2)%f.vertices.length];
				Edge nextEdge = curr.edgeWith(next);
				Edge prevEdge = curr.edgeWith(prev);
				Vertex[] subFace = new Vertex[4];
				subFace[0] = origPoints[this.vertices.indexOf(curr)];
				subFace[1] = edgePoints[this.edges.indexOf(nextEdge)];
				subFace[2] = facePoints[this.faces.indexOf(f)];
				subFace[3] = edgePoints[this.edges.indexOf(prevEdge)];
				cc.addFace(subFace);
			}
		}
		this.set(cc);
		return cc;
	}

	// Loop

	public Manifold loop() {
		// http://www.cs.cmu.edu/afs/cs/academic/class/15462-s12/www/lec_slides/lec07.pdf (1)
		// http://graphics.stanford.edu/~mdfisher/subdivision.html (2)
		// Note: triangulates non triangular faces first. 
		Manifold l = new Manifold(myParent);
		this.triangulate();
		// for each edge, an EDGE POINT is created
		Vertex[] edgePoints = new Vertex[this.edges.size()];
		for (Edge e : this.edges) {
			PVector edgePoint = new PVector();
			if (!e.boundary()) {
				for (Vertex v : e.left.vertices) {
					if (v != e.start && v != e.end) {
						edgePoint.add(PVector.mult(v.position, (float) 0.125));
						break;
					}
				}
				for (Vertex v : e.right.vertices) {
					if (v != e.start && v != e.end) {
						edgePoint.add(PVector.mult(v.position, (float) 0.125));
						break;
					}
				}
				edgePoint.add(PVector.mult(e.start.position, (float) 0.375));
				edgePoint.add(PVector.mult(e.end.position, (float) 0.375));
			}
			else {
				edgePoint = PVector.mult(e.start.position, (float) 0.5);
				edgePoint.add(PVector.mult(e.end.position, (float) 0.5));
			}
			edgePoints[this.edges.indexOf(e)] = l.addVertex(edgePoint);
		}

		// for each ORIGINAL POINT, create a new point
		Vertex[] origPoints = new Vertex[this.vertices.size()];
				for (Vertex v : this.vertices) {
					PVector origPoint;
					if (!v.boundary()) {
						float n = v.faces().length;
						float beta; // (2)
						if ( n > 3) {
							beta = 3 / (8 * n);
						}
						else {
							beta = (float) 0.1875;
						}
						//float beta = (1/n) * (0.625 - sq(0.375 + 0.25 * cos(2 * PI / n))); // Loop's original algorithm (1)
						//println(beta);
						origPoint = PVector.mult(v.position, 1 - n * beta);
						for (Face f : v.faces()) {
							origPoint.add(PVector.mult(f.vertices[(Arrays.asList(f.vertices).indexOf(v) + 1) % f.vertices.length].position, beta));
						}
					}
					else {
						origPoint = PVector.mult(v.position, (float) 0.75);
						for (Edge e : v.edges()) {
							if (e.boundary()) {
								if (v == e.start) {
									origPoint.add(PVector.mult(e.end.position, (float) 0.125));
								}
								else {
									origPoint.add(PVector.mult(e.start.position, (float) 0.125));
								}
							}
						}
					}
					origPoints[this.vertices.indexOf(v)] = l.addVertex(origPoint);
				}

				// draw faces
				for (Face f : this.faces) {
					for (Vertex v : f.vertices) { // this part ONLY works with triangular original meshes
						Vertex[] newf = new Vertex[3];
						newf[0] = origPoints[this.vertices.indexOf(v)];
						newf[1] = edgePoints[this.edges.indexOf(v.edgeWith(f.vertices[(Arrays.asList(f.vertices).indexOf(v) + 1) % f.vertices.length]))];
						newf[2] = edgePoints[this.edges.indexOf(v.edgeWith(f.vertices[(Arrays.asList(f.vertices).indexOf(v) + 2) % f.vertices.length]))];
						l.addFace(newf);
					}
					Vertex[] newf = new Vertex[f.vertices.length];
					for (int i = 0; i < f.vertices.length; i++) {
						newf[i] = edgePoints[this.edges.indexOf(f.vertices[i].edgeWith(f.vertices[(i+1)%f.vertices.length]))];
					}
					l.addFace(newf);
				}
				this.set(l);
				return l;
	}


	//---------------------------------------------------------//
	//                         Other...                        //
	//---------------------------------------------------------//

	public void toSphere() {
		// Project vertices onto a unit sphere.
		// Note: assumes that the centroid of the manifold is at the origin.
		for (Vertex v : this.vertices) {
			v.position.normalize();
		}
	}

	public void triangulate() {
		// Triangulate any faces with more than 3 sides
		// Draw fan around centroid
		Face[] originalFaces = this.faces.toArray(new Face[0]); // Clone
		for (Face f : originalFaces) {
			int n = f.vertices.length;
			if (n > 3) {
				//logger(this, "DEBUG", "triangulate; triangulating face " + f);
				Vertex centroid = this.addVertex(f.centroid());
				this.removeFace(f);
				for (int i = 0; i < n; i++) {
					//logger(this, "DEBUG", "triangulate; add new face");
					this.addFace(new Vertex[] {
							f.vertices[i], f.vertices[(i+1)%n], centroid
					}
							);
				}
			}
		}
	}
	
	/**
	 * flip edges where neighbouring triangles do not meet the Delaunay condition
	 */
	public int flipEdges() {
		// see http://en.wikipedia.org/wiki/Delaunay_triangulation#Visual_Delaunay_definition:_Flipping
		// TODO: use a recursive algorithm to ensure saturation?
		int tempSize = this.edges.size(); // avoid checking flipped edges
		for (int i = 0; i < tempSize; i++) {
			Edge e = this.edges.get(i);
			if (e.boundary()) continue;
			if (e.left.vertices.length == 3 && e.right.vertices.length == 3) {
				// get vertices A-D
				Vertex A = null, B = e.start, C = null, D = e.end;
				for (Vertex v : e.left.vertices) {
					if (v != B && v != D) {
						A = v;
						break;
					}
				}
				for (Vertex v : e.right.vertices) {
					if (v != B && v != D) {
						C = v;
						break;
					}
				}
				// get alpha and gamma
				// TODO: use circles (vectors) instead!!
				//float alpha = A.angleBetween(B, D);
				//float gamma = C.angleBetween(D, B);
				// check against Delaunay condition
				if (A.position.dist(e.right.circumcenter()) < e.right.circumradius()) {
					// flip edge
					Face ABD = e.left, BCD = e.right;
					this.removeFace(ABD);
					this.removeFace(BCD);
					this.addFace(new Vertex[]{B, C, A});
					this.addFace(new Vertex[]{D, A, C});
					// note: edge removed from the list and added to the end
					// must decrement the iterator (next edge now in current slot)
					i--;
					tempSize--;
				}
			}
		}
		return this.edges.size()-tempSize;
	}
	
	/**
	 * See {@link Vertex#laplacianSmoothing()}.
	 * @return The new "smooth" manifold.
	 */
	public Manifold laplacianSmoothing() {
		Manifold l = new Manifold(myParent);
		for (Vertex v : this.vertices) {
			l.addVertex(v.laplacianSmoothing());
		}
		l.edges = this.edges;
		l.faces = this.faces;
		this.set(l);
		return l;
	}

	public void set(Manifold m) {
		this.vertices = m.vertices;
		this.edges = m.edges;
		this.faces = m.faces;
	}

	// Debug...

	public void debug() {
		this.debug(true);
	}

	public void debug(boolean detail) {
		//logger(this, "INFO", "debug; printing topological/geometrical information...");

		myParent.println();

		myParent.println(" * Vertices: " + this.vertices.size());
		if (detail) {
			for (int i = 0; i < this.vertices.size(); i++) {
				Vertex v = this.vertices()[i];
				myParent.println("   " + i + ". " + v);
				myParent.println("     ~ " + v.position);
				myParent.print("     ~ " + v.edges.size() + " edge(s):");
				for (Edge e : v.edges) {
					myParent.print(" " + e);
				}
				myParent.println();
			}
			//println();
		}

		myParent.println(" * Edges: " + this.edges.size());
		if (detail) {
			for (Edge e : this.edges) {
				myParent.println("   - " + e);
				myParent.println("     ~ (l) " + e.left + ", (r) " + e.right);
				myParent.println("     ~ (s) " + e.start + ", (e) " + e.end);
			}
			//println();
		}

		myParent.println(" * Faces: " + this.faces.size());
		if (detail) {
			for (Face f : this.faces) {
				myParent.print("   - " + f + "");
				//for (Edge fe : f.edges) {
				//  print("\t" + fe);
				//}
				myParent.println();
			}
		}
		myParent.println();
	}

	//---------------------------------------------------------//
	//                         Export                          //
	//---------------------------------------------------------//

	// OBJ (http://en.wikipedia.org/wiki/Wavefront_.obj_file)

	public void exportOBJ() {
		//logger(this, "INFO", "exportOBJ; exporting OBJ file");
		PrintWriter obj = myParent.createWriter("export.obj");
		// vertices: "v x y z w"
		for (Vertex v : this.vertices) {
			obj.println("v " + v.position.x + " " + v.position.y + " " + v.position.z + " 1.0");
		}
		// vertex normals: "vn x y z w"
		for (Vertex v : this.vertices) {
			obj.println("vn " + v.normal().x + " " + v.normal().y + " " + v.normal().z);
		}
		// faces: "f v1 v2 v3 ..."
		for (Face f : this.faces) {
			obj.print("f");
			for (Vertex fv : f.vertices) {
				obj.print(" " + (this.vertices.indexOf(fv)+1) + "//" + (this.vertices.indexOf(fv)+1)); // vertices numbered from 1 (not 0)
			}
			obj.println();
		}
		obj.flush();
		obj.close();
	}

	// VRML (Indexed Face Sets: http://cs.iupui.edu/~aharris/mm/vrml4/vrml4.html)

	public void exportVRML() {
		//logger(this, "INFO", "exportVRML; exporting VRML file");
		PrintWriter vrml = myParent.createWriter("export.wrl");
		vrml.println("#VRML V2.0 utf8");
		vrml.println("Shape {");
		vrml.println("  geometry IndexedFaceSet {\n    coord Coordinate {");
		vrml.println("      point [");
		for (Vertex v : this.vertices) {
			vrml.println("        " + v.position.x + " " + v.position.y + " " + v.position.z + ",");
		}
		vrml.println("    }\n    coordIndex [");
		for (Face f : this.faces) {
			vrml.print("      ");
			for (int i = 0; i < f.vertices.length; i++) {
				vrml.print(this.vertices.indexOf(f.vertices[i]) + ",");
			}
			vrml.println("-1,"); // end face
		}
		vrml.println("    ]\n  }\n}");
		vrml.flush();
		vrml.close();
	}

	//---------------------------------------------------------//
	//                         Import                          //
	//---------------------------------------------------------//

	/**
	 * import a new manifold from an OBJ file
	 * @param file
	 * 			the file to import
	 * @return Manifold
	 */
	public Manifold importOBJ(String file) {
		//logger.info("loading manifold from file: " + file);
		Manifold m = new Manifold(myParent);
		//String file = myParent.selectInput("Select an OBJ file to import:");
		//myParent.println("You have selected: " + file);
		File f = new File(myParent.dataPath(file));
		if (!f.exists()) {
		  myParent.println("File does not exist");
		}
		else {
			for (String s : myParent.loadStrings(file)) {
				//logger(this, "DEBUG", "importOBJ; parsing line: \"" + s + "\"");
				String[] l = s.split(" ");
				if (l.length > 0) {
					// vertices
					if (l[0].equals("v")) {
						//logger.debug("vertex found");
						m.addVertex(new PVector(Float.valueOf(l[1]), Float.valueOf(l[2]), Float.valueOf(l[3])));
					}
					// faces
					else if (l[0].equals("f")) {
						//logger.debug("face found");
						int[] vertices = new int[l.length -1];
						for (int i = 1; i < l.length; i++) {
							vertices[i-1] = Integer.valueOf(l[i].split("/")[0])-1; // ignore texture-coordinates and normals
						}
						m.addFace(vertices);
					}
				}
			}
			this.set(m);
		}
		return m;
	}
	
	/**
	 * Create a prism with n-sided top/bottom, side length s and centred at the origin.
	 * @param theParent
	 * @param n
	 * @param s
	 * @return Manifold
	 */
	public static Manifold prism(PApplet theParent, int n, float s) {
		Manifold prism = new Manifold(theParent);

	    float a = 2 * PI / n; // angle
	    boolean ignoreSideLengths = (s == 0);
	    float R;
	    if (ignoreSideLengths) {
	      s = 1;
	      R = 0.7f;
	    }
	    else {
	      R = (0.5f * s) / theParent.sin(0.5f * a); // radius
	    }
	    
	    // Add vertices for top face (anticlockwise, looking down). 
	    for (int i = 0; i < n; i++) {
	      prism.addVertex(new PVector(R * theParent.cos(i * a), -0.5f * s, R * theParent.sin(i * a)));
	    }
	    // Add vertices for bottom face (anticlockwise, looking down).
	    for (int i = 0; i < n; i++) {
	      prism.addVertex(new PVector(R * theParent.cos(i * a), 0.5f * s, R * theParent.sin(i * a)));
	    }

	    // Add top/bottom faces.
	    int[] top = new int[n];
	    int[] bottom = new int[n];
	    for (int i = 0; i < n; i++) {
	      top[i] = i;
	      bottom[i] = 2 * n - (i + 1);
	    }
	    prism.addFace(top);
	    prism.addFace(bottom);

	    // Add side faces.
	    for (int i = 0; i < n; i++) {
	      prism.addFace(new int[] {
	        (i+1)%n, i, i+n, (i+1)%n+n
	      }
	      );
	    }
	    //logger.info("new prism created with " + n + " sides");
		return prism;
	}
    
    /**
     * Create a torus with 
     * @param theParent
     * @param R
     * 			the distance from the center of the tube to the center of the torus
     * @param r
     * 			the radius of the tube
     * @param toroidal
     * 			the number of divisions in the toroidal direction
     * @param poloidal
     * 			the number of divisions in the poloidal direction
     * @return the new manifold
     */
    public static Manifold torus(PApplet theParent, float R, float r, int toroidal, int poloidal) {
            Manifold t = new Manifold(theParent);
            for (int i = 0; i < poloidal; i++) {
                    for (int j = 0; j < toroidal; j++) {
                            float theta = i * PI*2/poloidal;
                            float phi = j * PI*2/toroidal;
                            float x = (R + r * theParent.cos(phi)) * theParent.cos(theta);
                            float y = (R + r * theParent.cos(phi)) * theParent.sin(theta);
                            float z = r * theParent.sin(phi);
                            t.addVertex(new PVector(x, y, z));
                    }
            }
            int n = toroidal * poloidal;
            for (int i = 0; i < poloidal; i++) {
                    for (int j = 0; j < toroidal; j++) {
                            int[] vertices = new int[]{
                                            ((i+1) * toroidal)%n + j,
                                            ((i+1) * toroidal)%n + (j+1)%toroidal,
                                            i*toroidal + (j+1)%toroidal,
                                            i*toroidal + j
                                            };
                            t.addFace(vertices);
                    }
            }
            return t;
    }
}

Topology.java

/**
 * 
 */
package manifold;

/**
 * abstract class for topological elements
 */
public abstract class Topology {
	
	public String toString() {
		return getClass().getSimpleName() + '@' + Integer.toHexString(hashCode());
	}
	
}

Vertex.java

/**
 * 
 */
package manifold;

import java.util.ArrayList;
import java.util.List;

import processing.core.*;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;

/**
 * a class for vertices
 */
public class Vertex extends Topology implements PConstants {

	PApplet myParent;
	public PVector position;
	public List<Edge> edges; // When you create a vertex, the number of edges is unknown.
	private int color;
	
	//private static final Logger logger = Logger.getLogger(Vertex.class.getName());

	public Vertex(PApplet theParent, PVector position) {
		this.myParent = theParent;
		this.position = position.get();
		this.edges = new ArrayList<Edge>();
		this.color = myParent.color(200, 200);
	}

	public void draw() {
		myParent.noStroke();
		myParent.fill(0);
		myParent.pushMatrix();
		myParent.translate(this.position.x, this.position.y, this.position.z);
		myParent.sphere((float) 0.02);
		myParent.popMatrix();
	}

	public void drawNormal() {
		myParent.pushMatrix();
		myParent.translate(this.position.x, this.position.y, this.position.z);
		PVector n = PVector.mult(this.normal(), 1);
		myParent.strokeWeight(1);
		myParent.stroke(0);
		myParent.line(0, 0, 0, n.x, n.y, n.z);
		myParent.popMatrix();
	}
	
	public void drawPrincipalCurvature() {
		myParent.pushMatrix();
		myParent.translate(this.position.x, this.position.y, this.position.z);
		PVector[] ct = this.curvatureTensors();
		ct[0].mult(4f);
		ct[1].mult(4f);
		myParent.stroke(200, 0, 0); // maximum
		myParent.line(-ct[0].x, -ct[0].y, -ct[0].z, ct[0].x, ct[0].y, ct[0].z);
		myParent.stroke(0, 0, 200); // minimum
		myParent.line(-ct[1].x, -ct[1].y, -ct[1].z, ct[1].x, ct[1].y, ct[1].z);
		myParent.popMatrix();
	}
	
	/**
	 * @param normal
	 * @param tensors
	 */
	public void draw(boolean normal, boolean tensors) {
		this.draw();
		if (normal) this.drawNormal();
		if (tensors) this.drawPrincipalCurvature();
	}

	public int color() {
		return color;
	}

	public void setColor(int color) {
		this.color = color;
	}

	/**
	 * Sorts the edges around a vertex in an anticlockwise direction
	 */
	public void sortEdges() {
		if (!this.edges.isEmpty()) {
			List<Edge> sorted = new ArrayList<Edge>();
			Edge e = this.edges.remove(0);
			sorted.add(e);
			int n = this.edges.size(); // Store this somewhere safe...
			boolean rs = false;
			for (int i = 0; i < n; i++) {
				Face f;
				if (rs == false) { // FORWARD SORT (default)
					//Edge e = sorted.get(i); // Last edge in sorted list.
					//logger(this, "DEBUG", "sortEdges; searching from " + e);
					if (this == e.start) {
						f = e.left;
					} 
					else {
						f = e.right;
					}
					if (f == null) { // Go back to first edge and reverse sort
						rs = true;
						e = sorted.get(0);
						i--;
						//logger(this, "DEBUG", "sortEdges; " + e + " is a boundary edge, returning to start and reversing");
						continue;
					}
					// Find the next edge (the one with face f on one side).
					for (Edge en: this.edges) {
						if (f == en.left || f == en.right) {
							sorted.add(en);
							this.edges.remove(en);
							//logger(this, "DEBUG", "sortEdges; edge found: " + en);
							e = en;
							break;
						}
					}
				}
				else { // REVERSE SORT (if boundary found) 
					//Edge e = sorted.get(0); // First edge in sorted list.
					//logger(this, "DEBUG", "sortEdges; searching from " + e);
					if (this == e.start) {
						f = e.right;
					} 
					else {
						f = e.left;
					}
					// Find the next edge (the one with face f on one side).
					for (Edge en: this.edges) {
						if (f == en.left || f == en.right) {
							sorted.add(0, en);
							this.edges.remove(en);
							e = en;
							//logger(this, "DEBUG", "sortEdges; edge found: " + en);
							break;
						}
					}
				}
			}
			// warn if there are edges left over (not a manifold)
			//if (this.edges.size() > 0) {
				//logger(this, "WARNING", "sortEdges; sort error - not a manifold! (" + this.edges.size() + " edge(s) left over.)");
			//}
			this.edges = sorted;
		}
	}

	/**
	 * returns array of edges connected to this vertex
	 * @return Edge[]
	 */
	public Edge[] edges() {
		return this.edges.toArray(new Edge[0]);
	}

	/**
	 * returns array of faces that include this vertex
	 * @return Face[]
	 */
	public Face[] faces() {
		// returns an array of faces to which this vertex belongs (ordered anticlockwise)
		List<Face> faces = new ArrayList<Face>();
		this.sortEdges();
		for (Edge e : this.edges) {
			if (this == e.start && e.right != null) {
				faces.add(e.right);
			} 
			else if (this == e.end && e.left != null) {
				faces.add(e.left);
			}
		}
		return faces.toArray(new Face[faces.size()]);
	}
	
	/**
	 * @return 1-ring neighbourhood
	 */
	public Vertex[] neighbourhood() {
		this.sortEdges();
		int n = this.edges.size();
		Vertex[] neighbours = new Vertex[n];
		for (int i = 0; i < n; i++) {
			neighbours[i] = this.edges.get(i).otherVertex(this);
		}
		return neighbours;
	}

	/**
	 * Find the edge that begins at this vertex and ends at b.
	 * @param b
	 * @return Edge
	 */
	public Edge edgeTo(Vertex b) {
		for (Edge e : this.edges) {
			if (e.start == this && e.end == b) {
				return e;
			}
		}
		return null;
	}

	/**
	 * Find the edge between this vertex and b, no direction assumed (if one exists).
	 * @param b
	 * @return Edge
	 */
	public Edge edgeWith(Vertex b) {
		Edge e = this.edgeTo(b);
		if (e != null) {
			return e;
		} else {
			return b.edgeTo(this);
		}
	}

	/**
	 * return the vertex normal
	 * @return PVector
	 */
	public PVector normal() {
		// sort edges --> .cross adjacent (normalised) pairs --> sum and normalise
		// average of face normals in 1-ring (weighted by area of face)
		// for quads and above a triangular face is approximated using the 1-ring
		this.sortEdges();
		PVector norm = new PVector();
		int n = this.edges().length;
		for (int i = 0; i < n; i++) {
		  norm.add(this.edges()[i].unitVectorFrom(this).cross(this.edges()[(i+1)%n].unitVectorFrom(this)));
		}
		return norm.normalize(null);
//		return this.meanCurvatureNormal();
	}

	/**
	 * returns true if this vertex lies on a boundary
	 * @return boolean
	 */
	public boolean boundary() {
		return (this.faces().length != this.edges().length);
	}

	/**
	 * returns a vector between this vertex and a specified vertex
	 * @param v
	 * 			the vertex at which the vector should end
	 * @return PVector
	 */
	PVector vectorFrom(Vertex v) {
		return PVector.sub(this.position, v.position);
	}
	/**
	 * returns a normalised vector between this vertex and a specified vertex
	 * @param v
	 * 			the vertex defining the vector's direction
	 * @return PVector
	 */
	PVector unitVectorFrom(Vertex v) {
		return this.vectorFrom(v).normalize(null);
	}

	/**
	 * returns the angle at this vertex between two edges
	 * @param a
	 * 			first edge
	 * @param b
	 * 			second edge
	 * @return float
	 */
	float angleBetween(Edge a, Edge b) {
		return PVector.angleBetween(a.unitVectorFrom(this), b.unitVectorFrom(this));
	}
	
	/**
	 * returns the angle formed by this and two other vertices
	 * @param a
	 * 			first vertex
	 * @param b
	 * 			last vertex
	 * @return float
	 */
	float angleBetween(Vertex a, Vertex b) {
		return PVector.angleBetween(this.vectorFrom(a), this.vectorFrom(b));
	}

	/**
	 * returns the cosine of the angle formed by this and two other vertices (quicker than calculating the angle itself)
	 * @param a
	 * 			first vertex
	 * @param b
	 * 			last vertex
	 * @return
	 */
	float cosAngleBetween(Vertex a, Vertex b) {
		return unitVectorFrom(a).dot(unitVectorFrom(b));
	}

	/**
	 * Compute the mean curvature at this vertex.
	 * @return \( H \) where \( H = {1 \over 2} (\kappa_1 + \kappa_2) \)
	 */
	public float meanCurvature() {
		PVector[] principals = this.curvatureTensors();
		return 0.5f * (principals[0].mag() + principals[1].mag());
	}

	/**
	 * Compute the gaussian curvature at this vertex.
	 * @return \( K \) where \( K = \kappa_1 \kappa_2 \)
	 */
	public float gaussianCurvature() {
		PVector[] principals = this.curvatureTensors();
		return principals[0].mag() * principals[1].mag();
	}
	
	/**
	 * Compute the principal directions and principal curvatures at this vector
	 * @return an array containing two PVectors, the first being the maximum principal direction and the second the minimum.
	 * The magnitude of each vector is equal to the corresponding principal curvature.
	 */
	public PVector[] curvatureTensors() {
		// Taubin tensors
		// see section 4 of http://pdf.aminer.org/000/234/737/curvature_approximation_for_triangulated_surfaces.pdf
		PVector[] ct = new PVector[2];
		PVector norm = this.normal();
		RealMatrix Nvi = new Array2DRowRealMatrix(new double[]{norm.x, norm.y, norm.z});
		RealMatrix I = MatrixUtils.createRealIdentityMatrix(3);
		
		this.sortEdges(); // TODO: dirty bit!
		
		// weighted areas, wij
		int n = this.edges.size();
		double[] wij = new double[n];
		double wijSum = 0;
		for (int i = 0; i < n; i++) {
			Edge e = this.edges()[i]; // current edge
			Edge prev = this.edges()[(i+(n-1))%n];
			Edge next = this.edges()[(i+1)%n];
			wij[i] = 0;
			if ((e.start == this && e.right != null) || (e.end == this && e.left != null)) { // face before edge
				//areas[i] += areas[(n+i-1)%n]/sumArea;
				wij[i] += 0.5 * prev.vectorFrom(this).cross(e.vectorFrom(this)).mag();
			}
			if ((e.start == this && e.left != null) || (e.end == this && e.right != null)) { // face after edge
				wij[i] += 0.5 * e.vectorFrom(this).cross(next.vectorFrom(this)).mag();
			}
			wijSum += wij[i];
		}
		
		// find matrix, Mvi
		RealMatrix MviEst = new Array2DRowRealMatrix(3, 3);
		for (Edge e : this.edges) {
			PVector edgeAsVector = this.vectorFrom(e.otherVertex(this));
			RealMatrix edgeAsMatrix = new Array2DRowRealMatrix(new double[]{edgeAsVector.x, edgeAsVector.y, edgeAsVector.z});
			RealMatrix Tij = (I.subtract(Nvi.multiply(Nvi.transpose()))).multiply(edgeAsMatrix);
			Tij = Tij.scalarMultiply(1/Tij.getFrobeniusNorm()); // TODO: double check
			double kij = Nvi.transpose().scalarMultiply(2).multiply(edgeAsMatrix).getEntry(0, 0) / Math.pow(e.length(), 2);
			RealMatrix AddToMviEst = Tij.multiply(Tij.transpose()).scalarMultiply((wij[this.edges.indexOf(e)]/wijSum)*kij);
			MviEst = MviEst.add(AddToMviEst);
		}
		
		// eigenvalues of Mvi (orthogonal)
		EigenDecomposition MviEiegVec = new EigenDecomposition(MviEst);
		double[][] v = new double[][]{
				MviEiegVec.getEigenvector(0).toArray(),
				MviEiegVec.getEigenvector(1).toArray(),
				MviEiegVec.getEigenvector(2).toArray()};
		double[] ev = MviEiegVec.getRealEigenvalues();
		
//		for (double vi : ev) {
//			System.out.print(vi + " ");
//			if (Math.abs(vi) < 0.000000001) System.out.print("(norm) ");
//		}
//		System.out.println();
		
		int[] index = new int[3]; // norm, max, min
		int[] mm = new int[2]; // separate min/max from normal
		
		// which one is the normal?!
		// TODO: check all three eigenvalues:
		// * min abs value is normal (closest to zero)
		// * if all three ~= 0, principal curvatures are zero vectors
		for (int i : new int[]{0, 1, 2}) {
			if (Math.abs(ev[i]) < 0.000000001) {
				mm = new int[]{(i+3-1)%3, (i+1)%3}; // TODO: send straight to 'index'?
				index[0] = i;
				break;
			}
		}
		
		// find min and max --> abs(max) > abs(min)
		if (Math.abs(ev[mm[0]]) > Math.abs(ev[mm[1]])) {
			index[1] = mm[0];
			index[2] = mm[1];
		}
		else {
			index[1] = mm[1];
			index[2] = mm[0];
		}
		
		// convert to PVectors and scale by curvature
		ct[0] = new PVector((float)v[index[1]][0], (float)v[index[1]][1], (float)v[index[1]][2]);
		//System.out.println(ct[0].mag());
		ct[0].mult((float)ev[index[1]]);
		ct[1] = new PVector((float)v[index[2]][0], (float)v[index[2]][1], (float)v[index[2]][2]);
		ct[1].mult((float)ev[index[2]]);
		
		return ct;
	}
	
	/**
	 * Laplacian smoothing.
	 * @return \( \displaystyle \bar{x}_{i}= \frac{1}{N} \sum_{j=1}^{N}x_j \)<br />
	 * Where \( N \) is the number of adjacent vertices to node \( i \) and \( \bar{x}_i \) is the new position for node \( x_i \).
	 */
	public PVector laplacianSmoothing() {
		PVector s = new PVector();
		Vertex[] n = this.neighbourhood();
		for (Vertex v : n) {
			s.add(v.position);
		}
		s.div(n.length);
		return s;
	}
}