From 96ed124d72488f389ee2e1ea2b47eaa18378c192 Mon Sep 17 00:00:00 2001 From: Michael Krisper <michael.krisper@tugraz.at> Date: Fri, 26 Feb 2016 17:32:01 +0100 Subject: [PATCH] corrected error in triangulation (edge.equals) --- VectoCore/Utils/DelauneyMap.cs | 24 ++++++++++----- VectoCore/Utils/VectoMath.cs | 56 +++++++++++++++++++--------------- 2 files changed, 48 insertions(+), 32 deletions(-) diff --git a/VectoCore/Utils/DelauneyMap.cs b/VectoCore/Utils/DelauneyMap.cs index bd1965528a..48a75f6326 100644 --- a/VectoCore/Utils/DelauneyMap.cs +++ b/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)). diff --git a/VectoCore/Utils/VectoMath.cs b/VectoCore/Utils/VectoMath.cs index 2f792a3e1d..fd0434295c 100644 --- a/VectoCore/Utils/VectoMath.cs +++ b/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 -- GitLab