Finding the Twins when Implementing Catmull-Clark subdivision using Half-Edge mesh

Question

Note: The description became a little longer than expected. Do you know a readable implementation of this algorithm using this mesh? Please let me know!

---

I'm trying to implement Catmull-Clark subdivision using Matlab, because later on the results have to be compared with some other stuff already implemented in Matlab. First try was with a Vertex-Face mesh, the algorithm works but it is not very efficient, since you need neighbouring information for edges and faces. Therefore, I'm now using a half-edge mesh. See also Designing a Data Structure for Polyhedral Surfaces, by Lutz Kettner (PDF link on that page).

My problem lies in finding the Twin HalfEdges, I'm just not sure how to do this. Below I'm describing my thoughts on the implementation, trying to keep it concise.

---

Half-Edge mesh (using indices to Vertices/HalfEdges/Faces):

Vertex (x,y,z,Outgoing_HalfEdge)
HalfEdge (HeadVertex (or TailVertex, which one should I use), Next, Face, Twin).
Face (HalfEdge)

To keep it simple for now, assume that every face is a quadrilateral. The actual mesh is a list of Vertices, HalfEdges and Faces. The new mesh will consist of NewVertices, NewHalfEdges and NewFaces, like this (note: Number_... is the number of ...):

NumberNewVertices: Number_Faces + Number_HalfEdges/2 + Number_Vertices
NumberNewHalfEdges: 4 * 4 * NumberFaces
NumberNewfaces: 4 * NumberFaces

Catmull-Clark:

Find the FacePoint (centroid) of each Face:
--> Just average the x,y,z values of the vertices, save as a NewVertex. 
Find the EdgePoint of each HalfEdge:
--> To prevent duplicates (each HalfEdge has a Twin which would result in the same HalfEdge)
--> Only calculate EdgePoints of the HalfEdge which has the lowest index of the Pair.
Update old Vertices

Ok, now all the new Vertices are calculated (however, their Outgoing_HalfEdge is still unknown). Next step to save the new HalfEdges and Faces. This is the part causing me problems!

Loop through each old Face, there are 4 new Faces to be created
  (because of the quadrilateral assumption)
First create the 4 new HalfEdges per New Face, starting at the FacePoint to the Edgepoint
Next a new HalfEdge from the EdgePoint to an Updated Vertex
Another new one from the Updated Vertex to the next EdgePoint
Finally the fourth new HalfEdge from the EdgePoint back to the FacePoint.

The HeadVertex of each new HalfEdge is known, the Next HalfEdge too. The Face is also known (since it is the new face you're creating!). Only the Twin HalfEdge is unknown, how should I know this?

By the way, while looping through the Vertices of the new Face, assign the Outgoing_HalfEdge to the Vertices. This is probably the place to find out which HalfEdge is the Twin.

Finally, after the 4 new HalfEdges are created, save the Face with the HalfVertex index the last newly created HalfVertex.

---

I hope this is clear, if needed I can post my (obviously not-yet-finished) Matlab code.

---

Edit: Thanks for moving my thread. I posted a link to the source code in a comment, please notice that this implementation considers general polygonal meshes, so not only quadrilateral meshes.

Furthermore, the twins of new HalfEdges 1 and 4 (red numbers in each new Face) are rather easy to find if you consider an old quadrilateral face divided into 4 new Faces (green numbers):

So, how to find the Twins of the 2- and 3 HalfEdges?

Answer 1

It seems the conceptual problem you've had is that you're trying to add half-edges one at a time and then wondering how they add up. Edges, though, are the real unit of modification, so you should always add them in pairs.

To implement (a single pass of) the algorithm, I'll annotate each of the model elements with a "newly created" flag, which indicates whether the element was created as a result of the algorithm. The top-level loop is to iterate on non-modified faces.

1. First ensure that each of the original edges of the face has been split. While doing this, create a list of each "new" vertex indicent to the face; these are the midpoints. -- a. To split an edge, we first locate the corresponding half-edge. Create a new vertex. We insert a new pair of half-edges into each linked list, adjusting the endpoints to the new vertex. Mark all four of the half-edges as new as well as the new vertex.

2. The first stop of subdividing the face is different than the others. Create a new vertex V to be in the middle of the old face and select a new vertex W incident to the face. We'll connect them as follows. Suppose the linked list of edges near W looks like ..aWb... Create a new pair of half edges c and d. Replace W in the linked list with WcVdW to make the list '..aWcVdWb..'. This creates a "floating edge" in the middle of the face. The data structure, however, ensures that we have a linked list of half-edges that represent the inner perimeter of the polygon. Mark vertex W and the half edges c and d as new.

3. Now for each remaining new vertex, we'll create a new pair of half-edges, but this time each pair will also create a new face. Select the previous ..cVdWb.. linked list sequence. Because all the original edge have been subdivided, that list extends to ..cVdWbXeYf... X is an old vertex and Y is a new one. Create a new pair of half edges g and g that will connect vertices V and Y. Extract the sequence VdWbXeY from the linked list and add g to it to create a new face [VdWbXeYg]. Add the half edge h to connect V and Y in the old face to make ..cVhYf... Mark the new face as new. If there are not more new vertices, we're done. If not, map the names ..cVhYf.. to ..cVdWb.. for iteration.

The notation is kind of nasty, but conceptually its pretty easy. Each of the three steps adds half-edges in pairs; in step 1 by dividing an edge, in steps 2 and 3 by adding them. Each of these additions keeps the incidence invariants of the polyhedron representation intact, which means you get improved locality of modification in the code.