/**
* Usage:
* towns = GSTownList();
* towns.Valuate(GSTown.GetLocation);
* neighbours = Delaunay(towns);
*/
require("graph.nut");
require("algebra.nut");
function CircleCondition(graph, edge, triangles) { // returns true, if the circumcircle of one triangle contains a point of another triangle
// calculate perpendicular bisectors
local perp1 = Perpendicular([[edge[0].x, edge[0].y], [edge[1].x, edge[1].y]]);
local perp2 = Perpendicular([[edge[0].x, edge[0].y], [triangles[0].x, triangles[0].y]]);
local center = Intersection2(perp1, perp2);
if (center == null) return true; // should never happen, but apparently it does
return Dist2sq(center, [edge[0].x, edge[0].y]) > Dist2sq(center, [triangles[1].x, triangles[1].y]);
}
function Flip(graph, edge) {
local triangles = graph.FindTriangles(edge[0], edge[1]);
if (triangles.len() < 2) return;
if (!CircleCondition(graph, edge, triangles)) return;
//GSLog.Info("Flip: <" + edge[0] + ", " + edge[1] + "> - <" + triangles[0] + ", " + triangles[1] + ">");
graph.RemoveEdge(edge[0], edge[1]);
graph.AddEdge(triangles[0], triangles[1]);
graph.flips++;
// recurse
if (graph.HasEdge(edge[0], triangles[0])) Flip(graph, [edge[0], triangles[0]]);
if (graph.HasEdge(edge[0], triangles[1])) Flip(graph, [edge[0], triangles[1]]);
if (graph.HasEdge(edge[1], triangles[0])) Flip(graph, [edge[1], triangles[0]]);
if (graph.HasEdge(edge[1], triangles[1])) Flip(graph, [edge[1], triangles[1]]);
}
function Delaunay(town_list) { // a list valuated by tileindex
GSLog.Info("Delaunay Triangulation starting");
local neighbours = Graph(); // undirected graph
local projections = []; // directed edges
local a = neighbours.AddNode(town_list.Begin()); // initialize with edge between first two towns
//GSLog.Info("Town: " + a + " " + GSTown.GetName(GSTile.GetClosestTown(a.ti)));
if (town_list.IsEnd()) return neighbours; // edge case: one town
local b = neighbours.AddNode(town_list.Next());
//GSLog.Info("Town: " + b + " " + GSTown.GetName(GSTile.GetClosestTown(b.ti)));
neighbours.AddEdge(a, b);
projections.append([a,b]);
if (a.y != b.y) projections.append([b,a]);
town_list.RemoveValue(a.ti);
town_list.RemoveValue(b.ti);
foreach (g, _ in town_list) {
local t = neighbours.AddNode(g);
//GSLog.Info("Town: " + t + " " + GSTown.GetName(GSTile.GetClosestTown(t.ti)));
local plane = t.y;
local new_projections = [];
local last = null;
local parallel_node = null;
// TODO: binary search of matching projections
foreach (edge in projections) {
local point = Intersection([[edge[0].x, edge[0].y], [edge[1].x, edge[1].y]], [[0,plane],[GSMap.GetMapSizeX(),plane]])
/*if (point != null) {
GSLog.Info("Intersection: <" + edge[0] + ", " + edge[1] + "> - (" + point[0] + ", " + point[1] + ") " + point[2]);
} else {
GSLog.Info("Intersection: <" + edge[0] + ", " + edge[1] + "> - parallel");
}*/
if (
(point == null) && (edge[0].y != plane) || // parallel
(point != null) && (
(point[2] <= 0) && (t.x <= point[0]) ||
(point[2] >= 0) && (t.x >= point[0])
)
) {
parallel_node = null;
neighbours.AddEdge(edge[0], t);
//GSLog.Info("Add: <" + edge[0] + ", " + t + ">");
neighbours.AddEdge(edge[1], t);
//GSLog.Info("Add: <" + edge[1] + ", " + t + ">");
if (neighbours.HasEdge(edge[0], edge[1])) Flip(neighbours, edge);
if (neighbours.HasEdge(edge[0], t)) Flip(neighbours, [edge[0], t]);
if (last == null) {
//GSLog.Info("Insert: <" + edge[0] + ", " + t + ">");
new_projections.append([edge[0], t]);
}
last = [t, edge[1]];
} else if (point == null) /*&& (edge[0].y == plane)*/ {
parallel_node = edge[1];
new_projections.append(edge);
} else {
parallel_node = null;
if (last != null) {
//GSLog.Info("Insert: <" + last[0] + ", " + last[1] + ">");
new_projections.append(last);
last = null;
}
if ((point[0] > 0) && (point[0] < GSMap.GetMapSizeX())) {
//GSLog.Info("Keep: <" + edge[0] + ", " + edge[1] + ">");
new_projections.append(edge);
}/* else {
GSLog.Info("Skip: <" + edge[0] + ", " + edge[1] + ">");
}*/
}
}
if (parallel_node != null) {
neighbours.AddEdge(parallel_node, t);
//GSLog.Info("Add: <" + parallel_node + ", " + t + ">");
last = [parallel_node, t];
}
if (last != null) {
//GSLog.Info("Insert: <" + last[0] + ", " + last[1] + ">");
new_projections.append(last);
}
projections = new_projections;
/*//local s="";
local prev = null;
foreach (i, edge in projections) {
if ((prev != null) && (prev[1] != edge[0])) GSLog.Error("Error: inconsistent projections");
prev = edge;
//if (i>0) s+=", ";
//s += edge[0];
}*/
//GSLog.Info("Graph: " + neighbours);
//GSLog.Info("Projections: ["+s+"]");
}
GSLog.Info("Delaunay Triangulation finished (" + neighbours.nodes.len() + " Towns, " + neighbours.flips + " Flips)");
return neighbours;
}