#### delaunay.nut

1. /**
2.  * Usage:
3.  *   towns = GSTownList();
4.  *   towns.Valuate(GSTown.GetLocation);
5.  *   neighbours = Delaunay(towns);
6.  */
7. require("graph.nut");
8. require("algebra.nut");
9.
10. function CircleCondition(graph, edge, triangles) { // returns true, if the circumcircle of one triangle contains a point of another triangle
11.     // calculate perpendicular bisectors
12.     local perp1 = Perpendicular([[edge[0].x, edge[0].y], [edge[1].x, edge[1].y]]);
13.     local perp2 = Perpendicular([[edge[0].x, edge[0].y], [triangles[0].x, triangles[0].y]]);
14.     local center = Intersection2(perp1, perp2);
15.     if (center == null) return true; // should never happen, but apparently it does
16.     return Dist2sq(center, [edge[0].x, edge[0].y]) > Dist2sq(center, [triangles[1].x, triangles[1].y]);
17. }
18.
19. function Flip(graph, edge) {
20.     local triangles = graph.FindTriangles(edge[0], edge[1]);
21.     if (triangles.len() < 2) return;
22.     if (!CircleCondition(graph, edge, triangles)) return;
23.     //GSLog.Info("Flip: <" + edge[0] + ", " + edge[1] + "> - <" + triangles[0] + ", " + triangles[1] + ">");
24.     graph.RemoveEdge(edge[0], edge[1]);
26.     graph.flips++;
27.     // recurse
28.     if (graph.HasEdge(edge[0], triangles[0])) Flip(graph, [edge[0], triangles[0]]);
29.     if (graph.HasEdge(edge[0], triangles[1])) Flip(graph, [edge[0], triangles[1]]);
30.     if (graph.HasEdge(edge[1], triangles[0])) Flip(graph, [edge[1], triangles[0]]);
31.     if (graph.HasEdge(edge[1], triangles[1])) Flip(graph, [edge[1], triangles[1]]);
32. }
33.
34. function Delaunay(town_list) { // a list valuated by tileindex
35.     GSLog.Info("Delaunay Triangulation starting");
36.     local neighbours = Graph(); // undirected graph
37.     local projections = [];     // directed edges
38.     local a = neighbours.AddNode(town_list.Begin()); // initialize with edge between first two towns
39.     //GSLog.Info("Town: " + a + " " + GSTown.GetName(GSTile.GetClosestTown(a.ti)));
40.     if (town_list.IsEnd()) return neighbours; // edge case: one town
42.     //GSLog.Info("Town: " + b + " " + GSTown.GetName(GSTile.GetClosestTown(b.ti)));
44.     projections.append([a,b]);
45.     if (a.y != b.y) projections.append([b,a]);
46.     town_list.RemoveValue(a.ti);
47.     town_list.RemoveValue(b.ti);
48.     foreach (g, _ in town_list) {
50.         //GSLog.Info("Town: " + t + " " + GSTown.GetName(GSTile.GetClosestTown(t.ti)));
51.         local plane = t.y;
52.         local new_projections = [];
53.         local last = null;
54.         local parallel_node = null;
55.         // TODO: binary search of matching projections
56.         foreach (edge in projections) {
57.             local point = Intersection([[edge[0].x, edge[0].y], [edge[1].x, edge[1].y]], [[0,plane],[GSMap.GetMapSizeX(),plane]])
58.             /*if (point != null) {
59.                 GSLog.Info("Intersection: <" + edge[0] + ", " + edge[1] + "> - (" + point[0] + ", " + point[1] + ") " + point[2]);
60.             } else {
61.                 GSLog.Info("Intersection: <" + edge[0] + ", " + edge[1] + "> - parallel");
62.             }*/
63.             if (
64.                 (point == null) && (edge[0].y != plane) || // parallel
65.                 (point != null) && (
66.                     (point[2] <= 0) && (t.x <= point[0]) ||
67.                     (point[2] >= 0) && (t.x >= point[0])
68.                 )
69.             ) {
70.                 parallel_node = null;
72.                 //GSLog.Info("Add: <" + edge[0] + ", " + t + ">");
74.                 //GSLog.Info("Add: <" + edge[1] + ", " + t + ">");
75.                 if (neighbours.HasEdge(edge[0], edge[1])) Flip(neighbours, edge);
76.                 if (neighbours.HasEdge(edge[0], t))       Flip(neighbours, [edge[0], t]);
77.                 if (last == null) {
78.                     //GSLog.Info("Insert: <" + edge[0] + ", " + t + ">");
79.                     new_projections.append([edge[0], t]);
80.                 }
81.                 last = [t, edge[1]];
82.             } else if (point == null) /*&& (edge[0].y == plane)*/ {
83.                 parallel_node = edge[1];
84.                 new_projections.append(edge);
85.             } else {
86.                 parallel_node = null;
87.                 if (last != null) {
88.                     //GSLog.Info("Insert: <" + last[0] + ", " + last[1] + ">");
89.                     new_projections.append(last);
90.                     last = null;
91.                 }
92.                 if ((point[0] > 0) && (point[0] < GSMap.GetMapSizeX())) {
93.                     //GSLog.Info("Keep: <" + edge[0] + ", " + edge[1] + ">");
94.                     new_projections.append(edge);
95.                 }/* else {
96.                     GSLog.Info("Skip: <" + edge[0] + ", " + edge[1] + ">");
97.                 }*/
98.             }
99.         }
100.         if (parallel_node != null) {
102.             //GSLog.Info("Add: <" + parallel_node + ", " + t + ">");
103.             last = [parallel_node, t];
104.         }
105.         if (last != null) {
106.             //GSLog.Info("Insert: <" + last[0] + ", " + last[1] + ">");
107.             new_projections.append(last);
108.         }
109.         projections = new_projections;
110.         /*//local s="";
111.         local prev = null;
112.         foreach (i, edge in projections) {
113.             if ((prev != null) && (prev[1] != edge[0])) GSLog.Error("Error: inconsistent projections");
114.             prev = edge;
115.             //if (i>0) s+=", ";
116.             //s += edge[0];
117.         }*/
118.         //GSLog.Info("Graph: " + neighbours);
119.         //GSLog.Info("Projections: ["+s+"]");
120.     }
121.     GSLog.Info("Delaunay Triangulation finished (" + neighbours.nodes.len() + " Towns, " + neighbours.flips + " Flips)");
122.     return neighbours;
123. }