corrected error in triangulation (edge.equals)

VectoCore/Utils/DelauneyMap.cs

@@ -31,8 +31,10 @@
 
 using System;
 using System.Collections.Generic;
+using System.Diagnostics;
 using System.Linq;
 using Newtonsoft.Json;
+using Org.BouncyCastle.Asn1.Crmf;
 using TUGraz.VectoCore.Exceptions;
 using TUGraz.VectoCore.Models;
 
@@ -71,23 +73,31 @@ namespace TUGraz.VectoCore.Utils
 			var triangles = new List<Triangle> { superTriangle };
 
 			// iteratively add each point into the correct triangle and split up the triangle
-			foreach (var point in Points) {
+			for (var i = 0; i < Points.Count; i++) {
+				var point = Points[i];
+
 				// If the vertex lies inside a triangle, the edges of the triangle are 
 				// added to the edge buffer and the triangle is removed from list.
 				var containerTriangles = triangles.FindAll(t => t.ContainsInCircumcircle(point));
-				triangles.RemoveAll(t => t.ContainsInCircumcircle(point));
+				triangles = triangles.Except(containerTriangles).ToList();
 
 				// Remove duplicate edges. This leaves the convex hull of the edges.
 				// The edges in this convex hull are oriented counterclockwise!
-				var convexHullEdges = containerTriangles.
-					SelectMany(t => t.GetEdges()).
-					GroupBy(edge => edge).
-					Where(group => group.Count() == 1).
-					SelectMany(group => group);
+				var allEdges = containerTriangles.SelectMany(t => t.GetEdges());
+				var groupedEdges = allEdges.GroupBy(edge => edge);
+				var convexHullEdges = groupedEdges.Where(group => group.Count() == 1).Select(group => group.Key);
 
 				var newTriangles = convexHullEdges.Select(edge => new Triangle(edge.P1, edge.P2, point));
 
 				triangles.AddRange(newTriangles);
+
+				// check invariant: m = 2n-2-k
+				// m...triangle count
+				// n...point count (index+1 +3 points on the supertriangle)
+				// k...points on convex hull (exactly 3 --> supertriangle)
+				if (triangles.Count != 2 * (i + 1 + 3) - 2 - 3) {
+					throw new VectoException("Triangulation invariant violated! Triangle count and point count doesn't fit together.");
+				}
 			}
 
 			_convexHull = triangles.FindAll(t => t.SharesVertexWith(superTriangle)).

VectoCore/Utils/VectoMath.cs

@@ -32,9 +32,7 @@
 using System;
 using System.Collections.Generic;
 using System.Diagnostics;
-using System.Diagnostics.Contracts;
 using TUGraz.VectoCore.Exceptions;
-using TUGraz.VectoCore.Models.Declaration;
 
 namespace TUGraz.VectoCore.Utils
 {
@@ -139,14 +137,14 @@ namespace TUGraz.VectoCore.Utils
 		public static List<double> QuadraticEquationSolver(double a, double b, double c)
 		{
 			var retVal = new List<double>();
-			var D = b * b - 4 * a * c;
+			var d = b * b - 4 * a * c;
 
-			if (D < 0) {
+			if (d < 0) {
 				return retVal;
-			} else if (D > 0) {
+			} else if (d > 0) {
 				// two solutions possible
-				retVal.Add((-b + Math.Sqrt(D)) / (2 * a));
-				retVal.Add((-b - Math.Sqrt(D)) / (2 * a));
+				retVal.Add((-b + Math.Sqrt(d)) / (2 * a));
+				retVal.Add((-b - Math.Sqrt(d)) / (2 * a));
 			} else {
 				// only one solution possible
 				retVal.Add((-b / (2 * a)));
@@ -233,7 +231,7 @@ namespace TUGraz.VectoCore.Utils
 
 		private bool Equals(Point other)
 		{
-			return X.Equals(other.X) && Y.Equals(other.Y) && Z.Equals(other.Z);
+			return X.IsEqual(other.X) && Y.IsEqual(other.Y) && Z.IsEqual(other.Z);
 		}
 
 		public override bool Equals(object obj)
@@ -285,30 +283,37 @@ namespace TUGraz.VectoCore.Utils
 
 		public Plane(Triangle tr)
 		{
-			var ab = tr.P2 - tr.P1;
-			var ac = tr.P3 - tr.P1;
+			var abX = tr.P2.X - tr.P1.X;
+			var abY = tr.P2.Y - tr.P1.Y;
+			var abZ = tr.P2.Z - tr.P1.Z;
 
-			var cross = ab.Cross(ac);
+			var acX = tr.P3.X - tr.P1.X;
+			var acY = tr.P3.Y - tr.P1.Y;
+			var acZ = tr.P3.Z - tr.P1.Z;
 
-			X = cross.X;
-			Y = cross.Y;
-			Z = cross.Z;
-			W = tr.P1.X * cross.X + tr.P1.Y * cross.Y + tr.P1.Z * cross.Z;
+			X = abY * acZ - abZ * acY;
+			Y = abZ * acX - abX * acZ;
+			Z = abX * acY - abY * acX;
+			W = tr.P1.X * X + tr.P1.Y * Y + tr.P1.Z * Z;
 		}
 	}
 
 	[DebuggerDisplay("Triangle(({P1.X}, {P1.Y}, {P1.Z}), ({P2.X}, {P2.Y}, {P2.Z}), ({P3.X}, {P3.Y}, {P3.Z}))")]
 	public class Triangle
 	{
-		public Point P1;
-		public Point P2;
-		public Point P3;
+		public readonly Point P1;
+		public readonly Point P2;
+		public readonly Point P3;
 
 		public Triangle(Point p1, Point p2, Point p3)
 		{
 			P1 = p1;
 			P2 = p2;
 			P3 = p3;
+
+			if ((P1.X.IsEqual(P2.X) && P2.X.IsEqual(P3.X)) || (P1.Y.IsEqual(P2.Y) && P2.Y.IsEqual(P3.Y))) {
+				throw new VectoException("triangle is not extrapolatable by a plane.");
+			}
 		}
 
 		/// <summary>
@@ -353,15 +358,15 @@ namespace TUGraz.VectoCore.Utils
 			var p2X = P3.X - p.X;
 			var p2Y = P3.Y - p.Y;
 
-			var p0square = p0X * p0X + p0Y * p0Y;
-			var p1square = p1X * p1X + p1Y * p1Y;
-			var p2square = p2X * p2X + p2Y * p2Y;
+			var p0Square = p0X * p0X + p0Y * p0Y;
+			var p1Square = p1X * p1X + p1Y * p1Y;
+			var p2Square = p2X * p2X + p2Y * p2Y;
 
 			var det01 = p0X * p1Y - p1X * p0Y;
 			var det12 = p1X * p2Y - p2X * p1Y;
 			var det20 = p2X * p0Y - p0X * p2Y;
 
-			var result = p0square * det12 + p1square * det20 + p2square * det01;
+			var result = p0Square * det12 + p1Square * det20 + p2Square * det01;
 
 			return result > 0;
 		}
@@ -420,6 +425,7 @@ namespace TUGraz.VectoCore.Utils
 	{
 		public readonly Point P1;
 		public readonly Point P2;
+
 		private Point _vector;
 
 		public Edge(Point p1, Point p2)
@@ -437,8 +443,8 @@ namespace TUGraz.VectoCore.Utils
 
 		protected bool Equals(Edge other)
 		{
-			return Equals(P1, other.P1) && Equals(P2, other.P2)
-					|| Equals(P1, other.P2) && Equals(P1, other.P2);
+			return (P1.Equals(other.P1) && Equals(P2, other.P2))
+					|| (Equals(P1, other.P2) && Equals(P2, other.P1));
 		}
 
 		public override bool Equals(object obj)
@@ -454,7 +460,7 @@ namespace TUGraz.VectoCore.Utils
 
 		public override int GetHashCode()
 		{
-			return (P1.GetHashCode()) ^ (P2.GetHashCode());
+			return P1.GetHashCode() ^ P2.GetHashCode();
 		}
 
 		#endregion