Merge pull request #8 from CptWesley/bowyer-watson#6

Bowyer-Watson Triangulation

Author: CptWesley

src/ClippingPlane.cs

@@ -4,19 +4,20 @@ namespace RadiantMapToWavefrontObj
 {
     public class ClippingPlane
     {
-        public readonly double A, B, C, D;
+        public readonly double D;
         public readonly Vector Normal;
 
+        public double A => Normal.X;
+        public double B => Normal.Y;
+        public double C => Normal.Z;
+
         // Constructor for a clipping plane.
         public ClippingPlane(Vertex v1, Vertex v2, Vertex v3)
         {
             Vector vector1 = new Vector(v2.X - v1.X, v2.Y - v1.Y, v2.Z - v1.Z).Unit();
             Vector vector2 = new Vector(v3.X - v1.X, v3.Y - v1.Y, v3.Z - v1.Z).Unit();
 
             Normal = Vector.CrossProduct(vector1, vector2).Unit();
-            A = Normal.X;
-            B = Normal.Y;
-            C = Normal.Z;                                                       // NB: A,B,C = Normal.X,Y,Z
             D = -(A * v2.X + B * v2.Y + C * v2.Z);
         }
 

src/Edge.cs

@@ -0,0 +1,78 @@
+
+using System;
+
+namespace RadiantMapToWavefrontObj
+{
+    public struct Edge
+    {
+        public Vertex A { get; set; }
+        public Vertex B { get; set; }
+
+        // Constructor for an edge between two vertices.
+        public Edge(Vertex a, Vertex b)
+        {
+            A = a;
+            B = b;
+        }
+
+        // Returns the vector given by this edge.
+        public Vector GetVector()
+        {
+            return (Vector) B - (Vector) A;
+        }
+
+        // Returns the length of this edge.
+        public double Length()
+        {
+            return GetVector().Length();
+        }
+
+        // Returns the inverse edge of this edge.
+        public Edge GetInverse()
+        {
+            return new Edge(B, A);
+        }
+
+        // Returns a stringified version of this object.
+        public override string ToString()
+        {
+            return "<" + A + ", " + B + ">";
+        }
+
+        // Checks if two edges are equal.
+        public override bool Equals(object obj)
+        {
+            if (obj is Edge)
+            {
+                Edge that = (Edge)obj;
+                if (A == that.A && B == that.B)
+                    return true;
+            }
+            return false;
+        }
+
+        // Returns a hascode for the object.
+        public override int GetHashCode()
+        {
+            return A.GetHashCode() + 2 * B.GetHashCode();
+        }
+
+        // Override == operator.
+        public static bool operator ==(Edge a, Edge b)
+        {
+            if (ReferenceEquals(a, null))
+            {
+                if (ReferenceEquals(b, null))
+                    return true;
+                return false;
+            }
+            return a.Equals(b);
+        }
+
+        // Override != operator.
+        public static bool operator !=(Edge a, Edge b)
+        {
+            return !(a == b);
+        }
+    }
+}

src/Face.cs

@@ -1,4 +1,5 @@
-
+using System;
+
 namespace RadiantMapToWavefrontObj
 {
     public class Face
@@ -11,12 +12,21 @@ public Face()
             _vertices = new Vertex[3];
         }
 
-        // Empty contructor that creates a new face with 3 vertices.
+        // Contructor that creates a new face with 3 vertices.
         public Face(Vertex[] vertices)
         {
             _vertices = vertices;
         }
 
+        // Contructor that creates a new face with 3 vertices.
+        public Face(Vertex a, Vertex b, Vertex c)
+        {
+            _vertices = new Vertex[3];
+            _vertices[0] = a;
+            _vertices[1] = b;
+            _vertices[2] = c;
+        }
+
         // Returns the vertices of the face.
         public Vertex[] GetVertices()
         {
@@ -66,6 +76,34 @@ public Vector GetNormal()
             return Vector.CrossProduct(v1, v2) * -1;
         }
 
+        // Finds the center and radius of a circumcircle of this triangle.
+        public Tuple<Vertex, double> GetCircumsphere()
+        {
+            // Find center.
+            Vector v0 = (Vector)_vertices[1] - (Vector)_vertices[0];
+            Vector v1 = (Vector)_vertices[2] - (Vector)_vertices[0];
+
+            Vector vx = Vector.CrossProduct(v0, v1);
+
+            Vector centerVector = (Vector.CrossProduct(vx, v0) * v1.SquareLength() + Vector.CrossProduct(v1, vx) * v0.SquareLength()) / (2 * vx.SquareLength());
+            Vertex center = _vertices[0] + centerVector;
+
+            // Find radius.
+            double radius = centerVector.Length();
+
+            return new Tuple<Vertex, double>(center, radius);
+        }
+
+        // Finds and returns the edges of this triangle.
+        public Edge[] GetEdges()
+        {
+            Edge[] edges = new Edge[3];
+            edges[0] = new Edge(_vertices[0], _vertices[1]);
+            edges[1] = new Edge(_vertices[1], _vertices[2]);
+            edges[2] = new Edge(_vertices[2], _vertices[0]);
+            return edges;
+        }
+
         // Returns a stringified version of the object.
         public override string ToString()
         {

src/ObjObject.cs

@@ -1,6 +1,5 @@
 using System;
 using System.Collections.Generic;
-using System.Linq;
 
 namespace RadiantMapToWavefrontObj
 {
@@ -142,50 +141,142 @@ private static Vertex[] FindIntersections(ClippingPlane[] planes)
         private static Face[] CreateFaces(Vertex[] vertices, ClippingPlane[] planes)
         {
             List<Face> faces = new List<Face>();
-            double centerOffset = 1e-4;
 
             for (int i = 0; i < planes.Length; ++i)
             {
-                // Create a center point, we need a tiny offset to make sure we don't mess up for perfectly symetrical shapes.
-                double centerX = centerOffset;
-                double centerY = centerOffset;
-                double centerZ = centerOffset;
-
-                List<Vertex> verts = planes[i].FindVerticesInPlane(vertices).ToList();
+                Vertex[] verts = planes[i].FindVerticesInPlane(vertices);
 
                 // Abort when there are no vertices anyway. Something went wrong...
-                if (verts.Count <= 0)
+                if (verts.Length < 3)
                     return null;
 
-                foreach (Vertex v in verts)
+                foreach (Face face in BowyerWatson(verts))
+                {
+                    FixNormal(face, planes[i].Normal);
+                    faces.Add(face);
+                }
+            }
+
+            return faces.ToArray();
+        }
+
+        // Apply Bowyer-Watson algorithm to triangulate all the points in a plane.
+        // Pseudo code taken from related wikipedia page and provided on the side in comments.
+        public static Face[] BowyerWatson(Vertex[] vertices)
+        {
+            List<Face> triangles = new List<Face>();
+
+            // Add super triangle to list.
+            Face superTriangle = FindSuperTriangle(vertices);
+            Vertex[] superVertices = superTriangle.GetVertices();
+            triangles.Add(superTriangle);
+
+            // Add points.
+            foreach (Vertex v in vertices)                                                              // for each point in pointList do
+            {
+                List<Face> badTriangles = new List<Face>();                                             // badTriangles := empty set
+                foreach (Face triangle in triangles)                                                    // for each triangle in triangulation do
                 {
-                    centerX += v.X;
-                    centerY += v.Y;
-                    centerZ += v.Z;
+                    if (InCircumsphere(v, triangle))                                                    // if point is inside circumcircle of triangle
+                        badTriangles.Add(triangle);                                                     // add triangle to badTriangles
                 }
 
-                Vertex center = new Vertex(centerX / verts.Count, centerY / verts.Count, centerZ / verts.Count);
-                center.SetNormal(planes[i].Normal);
+                List<Edge> polygon = new List<Edge>();                                                  // polygon := empty set
 
-                // Calculate faces based on some hackish algorithm that seems to work so far. Might need replacement later.
-                // Algorithm:  1. Find vertex closest to the center point.
-                //             2. Sort the list of vertices based on the distance to the vertex found in 1.
-                //             3. Add faces in the following manner: {0,1,2}, {1,2,3}, {2,3,4}, etc...
-                //             4. Check if the face's normal is in the right direction, otherwise: invert it.
-                if (verts.Count >= 3)
+                foreach (Face triangle in badTriangles)                                                 // for each triangle in badTriangles do
                 {
-                    verts.Sort((el1, el2) => center.Distance(el1).CompareTo(center.Distance(el2)));
-                    verts.Sort((el1, el2) => verts[0].Distance(el1).CompareTo(verts[0].Distance(el2)));
-                    for (int j = 0; j < verts.Count - 2; ++j)
+                    foreach (Edge edge in triangle.GetEdges())                                          // for each edge in triangle do
                     {
-                        Face face = new Face(verts.GetRange(j, 3).ToArray());
-                        FixNormal(face, planes[i].Normal);
-                        faces.Add(face);
+                        bool shared = false;
+                        foreach (Face otherTriangle in badTriangles)                                    // if edge is not shared by any other triangles in badTriangles
+                        {
+                            if (triangle == otherTriangle)
+                                continue;
+                            if (otherTriangle.GetEdges().Contains(edge) || otherTriangle.GetEdges().Contains(edge.GetInverse()))
+                            {
+                                shared = true;
+                                break;
+                            }
+                        }
+                        if (!shared)
+                            polygon.Add(edge);                                                          // add edge to polygon
                     }
                 }
+
+                foreach (Face triangle in badTriangles)                                                 // for each triangle in badTriangles do
+                    triangles.Remove(triangle);                                                         // remove triangle from triangulation
+
+                foreach (Edge e in polygon)                                                             // for each edge in polygon do
+                    triangles.Add(new Face(e.A, e.B, v));                                               // newTri := form a triangle from edge to point + add newTri to triangulation
             }
 
-            return faces.ToArray();
+            List<Face> result = new List<Face>();
+
+            foreach (Face t in triangles)                                                               // for each triangle in triangulation
+            {
+                Vertex[] curVertices = t.GetVertices();
+                if (!curVertices.Contains(superVertices[0])                                             // if triangle contains a vertex from original super-triangle
+                    && !curVertices.Contains(superVertices[1])
+                    && !curVertices.Contains(superVertices[2]))
+                    result.Add(t);                                                                      // remove triangle from triangulation
+            }
+
+            return result.ToArray();                                                                    // return triangulation
+        }
+
+        // Finds the Bowyer-Watson super triangle of a set of vertices.
+        private static Face FindSuperTriangle(Vertex[] vertices)
+        {
+            // Setup super triangle.
+            double minX, maxX, minY, maxY, minZ, maxZ;
+
+            minX = minY = minZ = Double.MaxValue;
+            maxX = maxY = maxZ = Double.MinValue;
+
+            foreach (Vertex v in vertices)
+            {
+                if (v.X < minX)
+                    minX = v.X;
+                if (v.X > maxX)
+                    maxX = v.X;
+
+                if (v.Y < minY)
+                    minY = v.Y;
+                if (v.Y > maxY)
+                    maxY = v.Y;
+
+                if (v.Z < minZ)
+                    minZ = v.Z;
+                if (v.Z > maxZ)
+                    maxZ = v.Z;
+            }
+
+            ClippingPlane plane = new ClippingPlane(vertices[0], vertices[1], vertices[2]);
+
+            Vertex a = new Vertex(minX, minY, minZ);
+            Vertex b = new Vertex(maxX, maxY, maxZ);
+
+            Vector ab = (Vector)(b - a);
+            a -= 10*ab;
+            b += 10*ab;
+
+            Vector triBase = Vector.CrossProduct(ab, plane.Normal).Unit();
+
+            double length = ((Vector)(b - a)).Length();
+
+            Vertex c = a + triBase * length;
+            Vertex d = a - triBase * length;
+
+            return new Face(b, c, d);
+        }
+
+        // Checks if a point lies in the circumsphere of a face.
+        private static bool InCircumsphere(Vertex point, Face face)
+        {
+            Tuple<Vertex, double> cs = face.GetCircumsphere();
+            if (point.Distance(cs.Item1) < cs.Item2)
+                return true;
+            return false;
         }
 
         // Fix normals of faces pointing in the wrong direction.

src/RadiantMapToWavefrontObj.csproj

@@ -108,6 +108,7 @@
     <Compile Include="ArrayExtension.cs" />
     <Compile Include="Brush.cs" />
     <Compile Include="ClippingPlane.cs" />
+    <Compile Include="Edge.cs" />
     <Compile Include="Face.cs" />
     <Compile Include="LineSegment.cs" />
     <Compile Include="ObjObject.cs" />

src/Vector.cs

@@ -33,6 +33,12 @@ public double Length()
             return Math.Sqrt(X*X+Y*Y+Z*Z);
         }
 
+        // Calculates the squared length of a vector.
+        public double SquareLength()
+        {
+            return X*X + Y*Y + Z*Z;
+        }
+
         // Returns a unit version of this vector.
         public Vector Unit()
         {

src/Vertex.cs

@@ -57,12 +57,24 @@ public bool OnPlane(ClippingPlane plane)
             return new Vertex(a.X + b.X, a.Y + b.Y, a.Z + b.Z);
         }
 
+        // Override + operator for vectors.
+        public static Vertex operator +(Vertex a, Vector b)
+        {
+            return new Vertex(a.X + b.X, a.Y + b.Y, a.Z + b.Z);
+        }
+
         // Override - operator.
         public static Vertex operator -(Vertex a, Vertex b)
         {
             return new Vertex(a.X - b.X, a.Y - b.Y, a.Z - b.Z);
         }
 
+        // Override - operator for vectors.
+        public static Vertex operator -(Vertex a, Vector b)
+        {
+            return new Vertex(a.X - b.X, a.Y - b.Y, a.Z - b.Z);
+        }
+
         // Override * operator for two points.
         public static Vertex operator *(Vertex a, Vertex b)
         {