我编写了一个使用Fortunes算法生成Voronoi图的程序。它在大多数情况下都有效，但是我发现一些特殊情况我不知道如何解决。
[其中之一是，如果两个共享相同y坐标和x坐标的站点关闭，则不会发生圆事件（结束2条边并开始一个新边的事件），或者在计算中的某个地方发生错误。
其他是2个站点共享x坐标和y坐标很近，因此没有弧交会发生。我该如何处理？
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class Fortune{
Voronoi diagram;
SortedQueue<AEvent> events;
Arc root;
float ly;
public Voronoi MakeDiagram()
{
// Broken seeds: 14, 15,
Random.InitState(16);
diagram = new Voronoi();
events = new SortedQueue<AEvent>(sortByY);
createSites(6);
constructDiagram();
return diagram;
}
private void createSites(int count)
{
for(int i = 0; i < count; i++)
{
events.Enqueue(new SiteEvent(new Vector2(Random.Range(0, 100), Random.Range(0, 100))));
}
}
private void constructDiagram()
{
while(events.Size() > 0)
{
events.Dequeue().handleEvent(this);
}
FinishEdges(root);
for(int i = 0; i < diagram.edges.Count; i++)
{
VEdge e = diagram.edges[i];
if(e.neighbour != null)
e.start = e.neighbour.end;
}
}
//Inser new parabol to beachline
public void AddParabol(Vector2 point)
{
ly = point.y;
Debug.DrawLine(point, point + Vector2.up, Color.blue, 100f);
if (root == null)
{
root = new Arc(point);
return;
}
if(root.isLeaf && root.site.y - point.y < 1)
{
root.isLeaf = false;
root.SetLeft(new Arc(root.site.point));
root.SetRight(new Arc(point));
VPoint mid2 = new VPoint((root.site.x + point.x) / 2, 100);
if(root.site.x > point.x)
{
root.edge = new VEdge(mid2, point, root.site.point);
}else
{
root.edge = new VEdge(mid2, root.site.point, point);
}
diagram.AddEdge(root.edge);
return;
}
Arc a = GetArcUnderPoint(root, point);
if(a.cevent != null)
{
events.Remove(a.cevent);
}
VPoint mid = new VPoint(point.x, GetY(a, point.x));
Debug.DrawLine(mid.point, mid.point + Vector2.up, Color.cyan, 100f);
VEdge l = new VEdge(mid, a.site.point, point);
VEdge r = new VEdge(mid, point, a.site.point);
diagram.AddEdge(l);
l.neighbour = r;
Arc lc = new Arc(a.site.point);
Arc m = new Arc(point);
Arc rc = new Arc(a.site.point);
a.SetLeft(lc);
Arc inter = new Arc();
inter.SetLeft(m);
inter.SetRight(rc);
a.SetRight(inter);
a.isLeaf = false;
a.edge = l;
inter.edge = r;
CheckCircleEvent(lc);
CheckCircleEvent(rc);
}
//Remove arc from tree
public void RemoveParabol(CircleEvent e)
{
ly = e.y;
Arc a = e.a;
Arc lp = a.leftParent;
Arc rp = a.rightParent;
Arc lc = lp.leftChild;
if (lc.cevent != null) events.Remove(lc.cevent);
Arc rc = rp.rightChild;
if (rc.cevent != null) events.Remove(rc.cevent);
VPoint end = new VPoint(e.x, GetY(a,e.x));
lp.edge.end = end;
rp.edge.end = end;
VEdge ne = new VEdge(end, lc.site.point, rc.site.point);
Arc higher = new Arc();
Arc cur = a;
while (cur != root)
{
cur = cur.parent;
if (cur == lp) higher = lp;
if (cur == rp) higher = rp;
}
higher.edge = ne;
diagram.AddEdge(ne);
Arc gp = a.parent.parent;
if(a.parent.left == a)
{
if (gp.left == a.parent) gp.SetLeft(a.parent.right);
if (gp.right == a.parent) gp.SetRight(a.parent.right);
}
else
{
if (gp.left == a.parent) gp.SetLeft(a.parent.left);
if (gp.right == a.parent) gp.SetRight(a.parent.left);
}
CheckCircleEvent(lc);
CheckCircleEvent(rc);
}
//Check if there is any circle events to come
public void CheckCircleEvent(Arc a)
{
Arc lp = a.leftParent;
Arc rp = a.rightParent;
if (lp == null || rp == null) return;
Arc lc = lp.leftChild;
Arc rc = rp.rightChild;
if (lc.site == rc.site) return;
VPoint s = Intersection(lp.edge, rp.edge);
if (s == null) return;
float r = Vector2.Distance(s.point, a.site.point);
if (s.y - r > ly) return;
CircleEvent e = new CircleEvent(a);
e.y = s.y - r;
e.x = s.x;
e.a = a;
a.cevent = e;
events.Enqueue(e);
}
//Find intersection of 2 edges and check if edge directions point to intersection point
public VPoint Intersection(VEdge left, VEdge right)
{
/*
* y = f1*x + g1
* y = f2*x + g2
* f1*x +g1 = f2*x+g2
* f1*x - f2*x = g2-g1
* x = (g2-g1)/(f1-f2)
*/
float x = (right.g - left.g) / (left.f - right.f);
float y = left.f * x + left.g;
Vector2 s = new Vector2(x, y);
if ((x - left.start.x) / left.dir.x < 0) return null;
if ((y - left.start.y) / left.dir.y < 0) return null;
if ((x - right.start.x) / right.dir.x < 0) return null;
if ((y - right.start.y) / right.dir.y < 0) return null;
return new VPoint(s);
}
// What if no intersection
public Arc GetArcUnderPoint(Arc root, Vector2 point)
{
// drawTree(root);
Arc curr = root;
while (!curr.isLeaf)
{
float x = ArcIntersectionX(curr, point.x);
Debug.Log("Returned =" + x);
if (float.IsNaN(x) || float.IsInfinity(x))
{
// drawTree(curr);
}
if (x > point.x) curr = curr.left;
else curr = curr.right;
}
return curr;
}
//Draws parabols to see where issueas are
public void drawTree(Arc a)
{
if(a.isLeaf)
{
for (int i = 0; i < 100; i++)
{
float y = GetY(a, i);
Vector2 t = new Vector2(i, y);
Debug.DrawLine(t, t + Vector2.up, Color.yellow, 100f);
}
}
if (!a.isLeaf)
{
drawTree(a.left);
drawTree(a.right);
}
}
// Find x coordinate of 2 arc intersection. Arcs are thel eft and right child of given node(Arc in tree)
public float ArcIntersectionX(Arc a, float x)
{
Arc lc = a.leftChild;
Arc rc = a.rightChild;
VPoint ls = lc.site;
VPoint rs = rc.site;
/*Intersection of 2 parabols
*
* y1 = 1 / 4f1 x ^ 2 - v11 / 2f1 x + v11 ^ 2 / 4f1 + v12
* y2 = 1 / 4f2 x ^ 2 - v21 / 2f2 x + v21 ^ 2 / 4f2 + v22
*
* 1 / 4f1 x ^ 2 - 1 / 4f2 x ^ 2 - v11 / 2f1 x + v21 / 2f2 x + v11 ^ 2 / 4f1 - v21 ^ 2 / 4f2 + v12 - v22 =
* (1 / 4f1 - 1 / 4f2) x ^ 2 + (-v11 / 2f1 + v21 / 2f2) x + (v11 ^ 2 / 4f1 - v21 ^ 2 / 4f2 + c12 - v22) = 0
*
* v11 = ls.x
* v21 = rs.x
*
* v12 = ls.y - f1
* v22 = rs.y - f2
*
* f1 = (ls.y - ly) / 2
* f2 = (rs.y - ly) / 2
*/
float f1 = (ls.y - ly) / 2;
float f2 = (rs.y - ly) / 2;
float a1 = (1/(4*f1) - 1/(4*f2));
float b1 = -ls.x / (2 * f1) + rs.x / (2 * f2);
float c1 = ls.x * ls.x / (4 * f1) - rs.x * rs.x / (4 * f2) + (ls.y - f1) - (rs.y - f2);
float D = Mathf.Sqrt(b1 * b1 - 4 * a1 * c1);
float x1 = (-b1 + D) / (2 * a1);
float x2 = (-b1 - D) / (2 * a1);
Vector2 p1 = new Vector2(x1, GetY(lc, x1));
Vector2 p2 = new Vector2(x2, GetY(lc, x2));
if(ls.y > rs.y)
{
return Mathf.Min(x1, x2);
}
else
{
return Mathf.Max(x1, x2);
}
}
public bool HasValue(float val)
{
return (!float.IsNaN(val) && !float.IsInfinity(val));
}
// Find where x = cons and arc intersect
public float GetY(Arc a, float x)
{
/*
* y = 1 / 4f (x - v1)^2 + v2
*
* f = (a.site.y - ly) / 2
* v1 = a.site.x
* v2 = a.site.y - f
*/
float f = (a.site.y - ly) / 2;
return 1 / (4 * f) * (x - a.site.x) * (x - a.site.x) + (a.site.y-f);
}
public void FinishEdges(Arc root)
{
if (root.isLeaf) return;
VEdge e = root.edge;
if(e.end == null)
{
float mx;
if (e.dir.x > 0) mx = Mathf.Max(100, e.start.x + 10);
else mx = Mathf.Min(0, e.start.x - 10);
float my = mx * e.f + e.g;
/*
if(my > 100)
{
my = 100;
mx = (100 - e.g) / e.f;
}
else if(my < 0)
{
my = 0;
mx = -e.g / e.f;
}
*/
VPoint end = new VPoint(mx, my);
e.end = end;
}
FinishEdges(root.left);
FinishEdges(root.right);
}
public static int sortByY(AEvent a, AEvent b)
{
if (a.y > b.y) return -1;
if (a.y < b.y) return 1;
if (a.x < b.x) return -1;
if (a.x > b.x) return 1;
return 0;
}
}
已解决
好，所以我调试了几个小时，发现问题出在以下事实：有时值是NaN或Infinity。因此，我只需要检查可能会给我其中之一的2。
我和您一样面临同样的问题。在执行https://github.com/fewlinesofcode/FortunesAlgorithm时。您可以检查代码并在那里找到答案，但是我也记录了代码，因此这里有一些摘录：