From 62b7153a8a8bc0622fa88cc5968a8e6a604ad0e4 Mon Sep 17 00:00:00 2001
From: wo80 <christian.woltering@tu-dortmund.de>
Date: Fri, 18 Feb 2022 16:27:20 +0100
Subject: [PATCH] Fix Voronoi failing for meshes with undead vertices.

---
 .../Smoothing/SimpleSmootherTest.cs           | 71 +++++++++++++++++++
 src/Triangle/Mesh.cs                          |  2 +-
 src/Triangle/Voronoi/VoronoiBase.cs           | 30 +++++---
 3 files changed, 94 insertions(+), 9 deletions(-)
 create mode 100644 src/Triangle.Tests/Smoothing/SimpleSmootherTest.cs

diff --git a/src/Triangle.Tests/Smoothing/SimpleSmootherTest.cs b/src/Triangle.Tests/Smoothing/SimpleSmootherTest.cs
new file mode 100644
index 0000000..318cb87
--- /dev/null
+++ b/src/Triangle.Tests/Smoothing/SimpleSmootherTest.cs
@@ -0,0 +1,71 @@
+using NUnit.Framework;
+using System.Collections.Generic;
+using System.Linq;
+using TriangleNet.Geometry;
+using TriangleNet.Meshing;
+using TriangleNet.Smoothing;
+
+namespace TriangleNet.Tests.Smoothing
+{
+    public class SimpleSmootherTest
+    {
+        [Test]
+        public void TestSmoothWithDuplicate()
+        {
+            var poly = GetPolygon();
+
+            var options = new ConstraintOptions() { ConformingDelaunay = true };
+
+            var quality = new QualityOptions()
+            {
+                MinimumAngle = 30.0
+            };
+
+            var mesh = poly.Triangulate(options, quality);
+
+            Assert.AreEqual(1, mesh.Vertices
+                .Where(v => v.Type == VertexType.UndeadVertex)
+                .Count());
+
+            quality.MaximumArea = 0.2;
+
+            mesh.Refine(quality, true);
+
+            Assert.AreEqual(1, mesh.Vertices
+                .Where(v => v.Type == VertexType.UndeadVertex)
+                .Count());
+
+            var smoother = new SimpleSmoother();
+
+            // Smooth mesh.
+            smoother.Smooth(mesh, 25);
+        }
+
+        private Polygon GetPolygon()
+        {
+            var poly = new Polygon(7);
+
+            var p = new List<Vertex>()
+            {
+                new Vertex(0.0, 0.0, 1),
+                new Vertex(2.0, 0.0, 1),
+                new Vertex(2.0, 1.0, 1),
+                new Vertex(2.0, 2.0, 1),
+                new Vertex(0.0, 2.0, 1),
+                new Vertex(1.5, 1.6, 1),
+                new Vertex(2.0, 0.0, 1) // duplicate
+            };
+
+            poly.Points.AddRange(p);
+
+            poly.Add(new Segment(p[0], p[1], 1));
+            poly.Add(new Segment(p[1], p[2], 1));
+            poly.Add(new Segment(p[2], p[3], 1));
+            poly.Add(new Segment(p[3], p[4], 1));
+            poly.Add(new Segment(p[4], p[0], 1));
+            poly.Add(new Segment(p[2], p[5], 2));
+
+            return poly;
+        }
+    }
+}
diff --git a/src/Triangle/Mesh.cs b/src/Triangle/Mesh.cs
index 07214c0..98147f6 100644
--- a/src/Triangle/Mesh.cs
+++ b/src/Triangle/Mesh.cs
@@ -383,7 +383,7 @@ namespace TriangleNet
         {
             numbering = NodeNumbering.None;
 
-            undeads = 0;               // No eliminated input vertices yet.
+            //undeads = 0;               // No eliminated input vertices yet.
             checksegments = false;     // There are no segments in the triangulation yet.
             checkquality = false;      // The quality triangulation stage has not begun.
 
diff --git a/src/Triangle/Voronoi/VoronoiBase.cs b/src/Triangle/Voronoi/VoronoiBase.cs
index 37978e2..ba328df 100644
--- a/src/Triangle/Voronoi/VoronoiBase.cs
+++ b/src/Triangle/Voronoi/VoronoiBase.cs
@@ -52,17 +52,17 @@ namespace TriangleNet.Voronoi
         /// Generate the Voronoi diagram from given triangle mesh..
         /// </summary>
         /// <param name="mesh"></param>
-        /// <param name="bounded"></param>
         protected void Generate(Mesh mesh)
         {
-            mesh.Renumber();
-
             base.edges = new List<HalfEdge>();
             this.rays = new List<HalfEdge>();
 
+            // Undead vertices cannot be Voronoi cell generators.
+            int count = mesh.vertices.Count - mesh.undeads;
+
             // Allocate space for Voronoi diagram.
             var vertices = new Vertex[mesh.triangles.Count + mesh.hullsize];
-            var faces = new Face[mesh.vertices.Count];
+            var faces = new Face[count];
 
             if (factory == null)
             {
@@ -74,10 +74,21 @@ namespace TriangleNet.Voronoi
             // Compute triangles circumcenters.
             var map = ComputeVertices(mesh, vertices);
 
-            // Create all Voronoi faces.
+            // Ensure linear numbering of vertices (excluding undeads).
+            int vid = 0;
+
+            // Create all Voronoi faces, skipping undead vertices.
             foreach (var vertex in mesh.vertices.Values)
             {
-                faces[vertex.id] = factory.CreateFace(vertex);
+                if (vertex.type == VertexType.UndeadVertex)
+                {
+                    vertex.id = count++;
+                }
+                else
+                {
+                    vertex.id = vid++;
+                    faces[vertex.id] = factory.CreateFace(vertex);
+                }
             }
 
             ComputeEdges(mesh, vertices, faces, map);
@@ -93,13 +104,16 @@ namespace TriangleNet.Voronoi
         /// Compute the Voronoi vertices (the circumcenters of the triangles).
         /// </summary>
         /// <returns>An empty map, which will map all vertices to a list of leaving edges.</returns>
+        /// <remarks>
+        /// This method will also change triangle ids (to ensure linear numbering of triangles).
+        /// </remarks>
         protected List<HalfEdge>[] ComputeVertices(Mesh mesh, Vertex[] vertices)
         {
             Otri tri = default(Otri);
             double xi = 0, eta = 0;
             Vertex vertex;
             Point pt;
-            int id;
+            int id, i = 0;
 
             // Maps all vertices to a list of leaving edges.
             var map = new List<HalfEdge>[mesh.triangles.Count];
@@ -107,7 +121,7 @@ namespace TriangleNet.Voronoi
             // Compute triangle circumcenters
             foreach (var t in mesh.triangles)
             {
-                id = t.id;
+                t.id = id = i++;
                 tri.tri = t;
 
                 pt = predicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);