Voronoi diagram, fortune algorithm special cases

Question

I wrote a program that uses Fortunes algorithm to generate Voronoi diagram. It works in most cases, but i have found some special cases i dont know how to fix.

One of them is that if 2 sites that share the same y-coordinate and x coordinates are close, then the circle event(event that ends 2 edges and starts a new one) does not happen, or somewhere in the calculations it goes wrong.

Other is that is 2 sites share the x coordinate and y coordinates are close, then no arc intersection happens. How do i handle that?

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;
    }
}

SOLVED

Ok so i debugged for few hours and found out that the issues come from the fact that sometimes the values are either NaN or Infinity. So i just have to place check where ever it might give me one of these 2.

Answer 1

I faced the same issues as you did. While implementing https://github.com/fewlinesofcode/FortunesAlgorithm. You can examine the code and find answers there, but i also documented it, so here are some exerpts:

Case: There is only one site in the **Beachline**. Exisiting site shares Y coordinate with new site

In this case two parabolas are degenerating into Rays origining at their event point and pointing up.
They have no intersection point but we know two things about the future of this case:
- The *x* coordinate of the intersection point point will be exactly in the middle
- *y* coordinate will lay somewhere far above the points.
We replace *y* with an arbitrary value big enough to cover our case. (`yVal`)

Also there is a degenerate case when the breakpoint falls exactly between two arcs.

Here is how i handled it (sorry, i used Swift language)

            /// Basically we're connecting **HalfEdge** pairs to create a connection like follows:
            ///     |
            ///     o
            ///    / \
            /// We also maintain proper Doubly linked list structures for each of Three neighbouring cells
            let vertex = Circle(
                p1: prev.point!,
                p2: newArc.point!,
                p3: next.point!
                )!.center

            prev.rightHalfEdge?.origin = vertex
            next.leftHalfEdge?.destination = vertex

            let lhe = diagram.createHalfEdge(newArc.cell!)
            newArc.leftHalfEdge = lhe
            lhe.origin = vertex
            let lTwin = diagram.createHalfEdge(prev.cell!)
            lTwin.destination = vertex
            makeTwins(lhe, lTwin)

            let rhe = diagram.createHalfEdge(newArc.cell!)
            newArc.rightHalfEdge = rhe
            rhe.destination = vertex
            let rTwin = diagram.createHalfEdge(next.cell!)
            rTwin.origin = vertex
            makeTwins(rhe, rTwin)

            connect(prev: lTwin, next: prev.rightHalfEdge)
            connect(prev: next.leftHalfEdge, next: rTwin)
            connect(prev: rhe, next: lhe)

            prev.rightHalfEdge = lTwin
            next.leftHalfEdge = rTwin