Adding edges dynamically in a DCEL/half-edge based graph?

Question

I'm trying to implement a vector graphics "drafting" system, where, in essence, users can draw lines on the screen and interact with the regions created by intersecting lines. I'm struggling with determining/evaluating what these regions are.

I've tried a few different solutions to this problem, chiefly keeping a list of edges and running a BFS to find shortest cycles, but this brought up a myriad of issues where the BFS would shortcut in illegal ways, and holes and degenerate edges caused more problems than I could count, so I moved on to a DCEL, half-edge system.

I've read seemingly everything I can on this topic, including two articles referenced frequently on here: http://kaba.hilvi.org/homepage/blog/halfedge/halfedge.htm and http://www.flipcode.com/archives/The_Half-Edge_Data_Structure.shtml. However, neither of these seem to answer this problem that I have when it comes to dynamically adding edges to the graph.

Let's say I start out with this single edge. Image

The half-edges connect with each other in a cycle, and the global, unbounded "outside face" is connected to one of the half edges. Easy, got that.

Then we add another edge, attached to the center vertex: Image

The new half-edges work fine, and we update the edges that flow into v1's next pointers to be the only other edges available that aren't their twins. Again, makes sense to me.

What confuses me to no end is what happens here, when we add a third edge to the center vertex: Image

I know that's what it's supposed to look like and link up to be, but I am so bewildered on how to achieve that programmatically, because I'm not sure how I can determine whether the edge (4,1) should point to edge (1,2), or edge (1,3) (similarly for what edge should point to (1,4)).

The answer seems obvious when looking at the image, but when you try to rationalize it in a robust, airtight algorithmic way, my brain melts and I can't figure it out. The textbook I'm reading (Computational Geometry, Mark de Berg et al., pg 35), just says

"[to test where the edge] should be in the cyclic order of the edges around vertex v".

The algorithm given in the hilvi.org article for finding the outgoing and incoming edges to link up with doesn't even seem to work, as it will take vertex 1, and follow its outgoing edge's twin until it finds a "free" edge, which in this case, is (2,1) which is wrong. (Unless I'm understanding it incorrectly, I could be understanding this whole problem wrong.)

So I'm absolutely stumped. My only idea now is to create some sort of heading attribute for each half edge, where I measure the angle created by the edge, and choose the edges that way, and maybe that's right, but that seems so against what the half-edge structure seems to support, at least in the articles I'm reading about it, nothing seems to mention anything like that. Any help would be extremely appreciated. I've been on this problem for over a week now and just can't seem to get unstuck.

Answer 1

Right, so I've spent a lot of time thinking about this problem, and to be honest I'm kinda surprised that I can't find a direct answer to this issue. So, in case anyone in the future runs into a similar problem of wanting to populate a half-edge graph from the ground up, here's a solution that works. I don't have a blog, so I'm writing it here.

I have no idea if it's the best answer, but it works in linear time and seems simple to me.

I'll be dealing with the following objects/classes, that vary slightly from the conventional DCEL:

class Vertex {
    x;
    y;

    edges = []; //A list of all Half Edges with their origin at this vertex.
                //Technically speaking this could be calculated as needed, 
                  and you could just keep a single outgoing edge, but I'm not 
                  in crucial need of space in my application so I'm just 
                  using an array of all of them.
}


class HalfEdge {
    origin; //The Vertex which this half-edge emanates from

    twin; // The half-edge pair to this half-edge

    face; // The region/face this half-edge is incident to

    next; // The half-edge that this half-edge points to
    prev; // The half-edge that points to this half-edge

    angle; //The number of degrees this hedge is CW from the segment (0, 0) -> (inf, 0)
}


class Face {
    outer_edge; //An arbitrary half-edge on the outer boundary defining this face.
    inner_edges = []; //A collection of arbitrary half-edges, each defining
                      //A hole in the face.

    global; //A boolean describing if the face is the global face or not.
            //This could also be done by having a single "global face" Face instance. 
            //This is simply how I did it.
}

For initializing a vertex at(x,y):

1. Verify that a vertex with the given (x,y) coordinates does not already exist. If it does, you don't have to do anything (except maybe return this existing vertex if you're using it immediately).

2. If it doesn't, allocate space for and create a new vertex with the corresponding x,y values, and with its incident edge null.

For initializing an edge from vertex A to vertex B:

1. Similarly to many of the articles about this topic, we create two new instances of HalfEdge, one from vertex A to B, one from B to A. They link to each other in that we set their twin, prev, and next pointers all to the other half-edge (hedge).

2. We also set the angle of the hedge. The angle is calculated clockwise from the positive x-axis. The function I implemented is below. This is super important to making this data structure work properly, and that fact that I haven't read anything in the literature about this being important makes me think there has to be a better way, but I digress.

       setAngle(){
           const dx = this.destination().x - this.origin.x;
           const dy = this.destination().y - this.origin.y;
   
           const l = Math.sqrt(dx * dx + dy * dy);
           if (dy > 0) {
               this.angle = toDeg(Math.acos(dx / l));
           } else {
               this.angle = toDeg(Math.PI * 2 - Math.acos(dx / l));
           }
   
            function toDeg(rads) {
                return 180 * rads / Math.PI;
            }
   
        }
   

3. Next, we pair the vertices with their new edges by adding them to their Vertex's edge list, and then we sort the edge list by the hedge's angles from smallest (0) to largest (359).

4. Then and this is the crucial step, in order to link everything up properly, we grab the closest hedge to the new hedge we're trying to link up in CCW order. Basically, wherever our new hedge ends up in the edge list, it's that index - 1 (if index = 0, we return edges[edges.length - 1]). Take that edge's twin, and that becomes our AIn described in the hivli article above. BOut = AIn.next.

5. We set AIn.next = hedgeAB and similarly, hedgeAB.prev = AIn, then hedgeBA.next = AOut, AOut.prev = hedgeBA. Perform steps 3-5 for the hedgeBA as well, except running the CCW search on vertex B.

6. Then, if both vertex A and B were "old" vertices, meaning their edge lists now have at least 2 elements each, a new face has potentially been added, and we need to find it (the edge case is having two isolated edges and connecting them to create an unbounded bucket or cap shape)

For initializing a Face:

1. We need to find all the cycles in a graph. For my first implementation of this, I recalculated all of the cycles every single time, resetting all the faces. This isn't necessary, but it also isn't too expensive, since we're not running a search, everything is in linear time with respect to the number of cycles and number of vertices in each cycle.

2. To do this, we get a list of all the hedges in the graph. It doesn't really matter how you do this, I decided to keep an array of every hedge that I passed into my cycle-finder function each time.

3. Then we look through that list, while the list isn't empty, we take the first item and run its cycle, removing every hedge we find along the way from the list, and adding it into a new cycle, which we add to another list

4. With this new list of cycles, we need to determine whether the cycle is an inside/outside cycle. There are a lot of ways to do this, and the Computational Geometry book mentioned above has a great section about it. The one that I used is to calculate the area defined by each cycle. If the area is >= 0, the cycle is defined by "inside" hedges. Otherwise, it's defined by "outside" hedges.

5. The last step is to set all the face records, again, the aforementioned textbook has a lot of great details about this, but the basic idea is to basically create a virtual "graph" of these cycles, and connect outside cycles (which are holes in faces), to their corresponding inside cycles, which are the outer boundaries of faces. To do this, you look at the leftmost vertex of the cycle and extend a ray out infinitely to the left, and "connect" the cycle with the first downward-facing hedge of a cycle the ray hits (I'll leave the implementation up to you, I have no idea if my way is the best, in short, I checked every cycle with a leftmost vertex left of the current cycle and calculated the rightmost intersection with the y value of the leftmost vertex of the current cycle, and then checked if that was facing downward).

6. With this graph of cycles, run a BFS/DFS starting from each "inside-hedge" cycle (not the holes), and create a face with an arbitrary hedge from the inside-hedge cycle as the outer edge, (or null if it's the global face), and an arbitrary hedge from each found hole-cycle to the face's inner components.

Hey presto, that's it. If you check everything every time, that handles everything. It handles face splitting like a charm and is very robust and quick. I don't know if it's right, but it works.