Mesh Edge Construction

In this post I will present three equivalent algorithms to compute the edges of a polygonal mesh. The algorithms are subsequent optimization steps to produce the exact same result with progressively increased efficiency. In the first part of this text, I will describe the typical representation of a mesh topology, distinguish between different concepts of edges.

To better define the subject, picture the wireframe of a polygonal mesh. Assuming there is no overlapping geometry in the scene, if you’d have to draw each single segment once, you would be drawing edges. That would be different than drawing polygons outline, because in drawing adjacent polygons you would end up drawing some, if not most of the segments twice.

What this post is not

This post is not about topological acceleration structure, like half-edges, winged-mesh, corner table, vertex-vertex, or similar. I borrow some of the terms though. I will refer to half-edges throughout the text, to distinguish from unique edges, I describe the input mesh in the very common face-vertex mesh format.

Introduction

Most of 3D graphics revolves around polygons in one way or another. The cool kids today work with distance fields and Gaussian splats. But let’s forget about that, and NURBS and parametric surfaces for a moment. Chances are, if you’re working on a graphics project, you’re dealing with polygonal meshes.

The standard way to describe a mesh is through a set of points connected in a specific winding order, often referred to as the topology. The simplest possible representation is called triangle soup, where a set of point triplets is defined, with each triplet representing the vertices of a triangle. That’s it. It’s the simplest but also the most wasteful and impractical way to describe a mesh, so let’s skip a few history pages (forget about fans and strips). The next and more useful way to describe a mesh is via a set of points and indices referencing those points. Let’s define a square:

float3 points[] = {(0, 0, 0), (1, 0, 0), (0, 0 1), (1, 0, 1)};
int   indices[] = {0, 1, 2, 1, 2, 3};

The indices define the connection order between points. In this case, assuming triangles, we have six indices forming two triangles. This gives us a square lying in the XZ plane.

Triangle meshes are the lowest common denominator—they are practical for rendering due to extensive hardware support, but they are not very convenient for modeling or manipulation. One way to improve this is to stop assuming that polygons are just triangles. In fact, the word polygon comes from Greek and means many corners. But how many? To define this explicitly, a third piece of information is added to specify how many polygons there are and how many sides each has. This data structure is called Face-Vertex Mesh. For example:

float3 points[] = {(0, 0, 0), (1, 0, 0), (0, 0, 1), (1, 0, 1)};
int   indices[] = {0, 1, 3, 2};
int faceCount[] = {4};

Here, we have an equivalent square with the same points as before, but this time, faceCount tells us there is a single polygon with four indices.

There are many other interesting properties to consider. For example, what is the winding order of the indices in relation to the position—do they rotate clockwise or counterclockwise? But let’s leave that for now. Instead, let’s add a second polygon, another square adjacent to the first, and see what we can observe.

float3 points[] = {(0, 0, 0), (1, 0, 0), (2, 0, 0),
                   (0, 0, 1), (1, 0, 1), (2, 0, 1)};
int   indices[] = {0, 1, 4, 3, 1, 2, 5, 4};
int faceCount[] = {4, 4};

We have two polygons, 6 points, 7 edges. Where do these come from? Each square has 4 sides, but the two squares are adjacent along one edge. The edge formed between point 1 and 4 is shared.

polygon 0: {0, 1, 4, 3} -> edges {(0, 1), (1, 4), (4, 3), (3, 0)}
polygon 1: {1, 2, 5, 4} -> edges {(1, 2), (2, 5), (5, 4), (4, 1)}

The vertical edge in the middle appears twice in the topology: first as the index pair (1, 4) and second as (4, 1). To distinguish between these, we refer to them as opposite half-edges. This term is borrowed from the widely used half-edge data structure in mesh processing. You can think of it as if each polygon owns half of the edge.

Opposite half-edges have the key property that their index order is swapped. In diagrams, they are often represented as half-arrows, following the winding order of the polygon indices. You can break this convention by inverting the winding order of one of the two polygons. When that happens, you create a type of mesh that most modeling applications—and many mesh-processing algorithms—struggle with. This class of meshes is called non-2-manifold, or simply non-manifold.

For the purpose of this text, I define a half-edge simply as a pair of indices, without assuming any connectivity properties. This distinction helps clarify the difference between an edge and a half-edge. So, the index pairs (1, 4) and (4, 1) represent two half-edges, but they still correspond to the same edge. Intuitively, half-edges are derived from the topology: they are determined by the indices, with each half-edge consisting of an index and the next one in the polygon’s winding order.

float3 points[] = ...
int   indices[] = ...
int faceCount[] = ...
int nFaces = N;

for (int faceIndex = 0, offset = 0; faceIndex < nFaces; ++faceIndex)
{
    const int sides = faceCount[faceIndex];
    const int* polygonIndices = indices + offset;
    offset += sides;

    for (int vertex = 0; vertex < sides; vertex++)
    {
        int nextVertex = (vertex + 1) % sides;
        int indicesPair[2] = { polygonIndices[vertex], polygonIndices[nextVertex] };
        
        // Store pair ...
    }
}

This listing doesn’t perform any operations—it simply illustrates how to interpret indices, faceCount, and winding order to construct index pairs, which serve as the base data for half-edges.

What else can we observe? Points 0, 2, 3, and 5 are each connected to a single face, while points 1 and 4 are connected to two faces. A similar observation applies to edges: points 0, 2, 3, and 5 are each connected to two edges, whereas points 1 and 4 are connected to three. These properties are known as point face-valence and point edge-valence.

Closely related is the concept of adjacency—for example, what is the list of faces connected to a point with valence N? Or what are the edges (or opposite indices) associated with a point of valence M? These properties will prove useful later on and will become key insights in designing efficient algorithms!

Unique Edges

Any given mesh has a number of half-edges equal to the number of indices. In the example mesh above, there are 8 indices, meaning there are 8 half-edges. However, there are only 7 edges—7 unique segments that a user can count.

Now, suppose I want to implement a modeling function in a 3D editor where users can click on edges to perform operations. Or perhaps, given a mesh, I need to generate edge indices for wireframe rendering. Or maybe I am implementing subdivision surfaces and I want the user to specify edge creases. How do I determine the number of unique edges? How do I construct that data?

Defining Unique Edges

The first key observation is that opposite half-edges have their index pairs swapped. For example, the polygons containing (1, 4) and (4, 1) both reference the same edge. If we sort each pair in ascending order, both representations become (1, 4). By consistently defining edges as sorted index pairs, we establish a definitive way to identify unique edges. I refer to the lower of the two values the minor index, the other the major index.

// marked with * the swapped pairs
polygon 0: {0, 1, 4, 3} -> edges {(0, 1), (1, 4), (3, 4)*, (0, 3)*}
polygon 1: {1, 2, 5, 4} -> edges {(1, 2), (2, 5), (4, 5)*, (1, 4)*}

Structuring Edge Data

To efficiently store edges, we can define a supporting structure where each edge is represented as a sorted index pair. This allows for easy comparisons and lookups. Instead of a struct, I use a union to store the indices as a single 64-bit unsigned integer, making comparisons more efficient.

// An edge expressed as two indices to points in some geometry. The two indices are sorted
// in ascending order to allow to identify edges shared between faces.
union Edge
{
    inline Edge() {}
    inline Edge(uint32_t v0, uint32_t v1)
    {
        if (v0 < v1) v[0] = v0, v[1] = v1;
        else         v[0] = v1, v[1] = v0;
    }

    inline bool operator  < (const Edge& other) const { return code < other.code; }
    inline bool operator == (const Edge& other) const { return code == other.code; }
    inline bool operator != (const Edge& other) const { return code != other.code; }
    inline const Edge& operator = (const Edge& other) { code = other.code; return *this; }

    struct
    {
        uint32_t v[2];
    };
    uint64_t code;
};

Eliminating Duplicates & Enumerating Edges

Now, we need a way to create a unique set of edges to remove duplicates. Additionally, we might want to assign each edge a unique index—for example, which of these is edge 3?

A natural and principled way to establish a definite order is to count edges as we encounter them in the topology. This ensures consistency and makes edge indexing straightforward.

The map approach

Let’s start with the laziest possible approach—using an associative container such as a map or hash map. In this example, we use a map to detect duplicates and store the results in a vector called edges, without keeping the map afterward.

The initial algorithm constructs edges by iterating over pairs of point indices for each face. It inserts each edge into a map (of the unordered variety), ensuring only the first occurrence of each edge is retained while duplicates are discarded.

Since insertion into a map has a unit cost of O(log n), the overall complexity of this approach is O(n log n). This makes it relatively expensive compared to other possible methods.

// My initial algorithm makes edges by iterating faces point indices pairs, insert into a map
// (of the unordered variety) to keep only the first instances of those edges and discard duplicates.
// This makes the algorithm O(n log n) in complexity. The unit cost of inserting into the map is high.

// Prepare edges data. An edge is defined as the two indices of the two points forming the edge.
// Edges are unique and shared across faces.
std::vector edges; //< The result

uint32_t nextEdge = 0;
std::map edgesMap;
for (int fIndex = 0, offset = 0; fIndex < nFaces; ++fIndex)
{
    int count = faceCounts[fIndex];
    for (int vIndex = 0; vIndex < count; vIndex++)
    {
        uint32_t v0 = indices[offset + vIndex];
        uint32_t v1 = indices[offset + ((vIndex + 1) % count)];
        Edge edge(v0, v1);

        auto [iterator, newEntry] = edgesMap.insert({ edge, nextEdge });
        if (newEntry)
        {
            edges.push_back(edge);
            nextEdge++;
        }
        else
        {
            // If you need to associate the existing edge index to an half edge (or the topology)
            // that can be accessed as iterator->second.
        }
    }
    offset += count;
}

Why Use a Map Instead of a Set?

Technically, we could use a set instead of a map, but there are cases where we may need to retrieve the index of an edge associated with a half-edge. Using a map allows access to iterator->second, which can be useful when integrating edges with the topology.

For better performance, the map can be replaced with a hash-based container such as std::unordered_map. This would reduce the lookup time from O(log n) to O(1) on average, improving efficiency. I leave this as an exercise for the reader.

Flat vector + sort

Associative containers like map and unordered_map are often inefficient because they allocate memory in a granular way. Even if the insertion algorithm may be efficient in theory, the execution is often not due to memory management. The following algorithm retains the same core concept but implements it with a simpler memory layout. The code itself is longer, but the amount of executed code is smaller, leading to better performance.

// Similar to the algorithm using an unordered_map, but now using flat vectors.
// This approach requires more steps and a higher memory peak but achieves a 50% speed improvement.
// The overall complexity remains O(n log n).
struct EdgeEntry
{
    Edge e;
    uint32_t entryIndex;
    bool operator==(const EdgeEntry& other) const { return e.code == other.e.code; }
    bool operator!=(const EdgeEntry& other) const { return e.code != other.e.code; }
};


const int nIndices = indices.size();
const int nFaces   = faceCounts.size();

std::vector halfEdges;
halfEdges.resize(nIndices);

for (int fIndex = 0, offset = 0; fIndex < nFaces; ++fIndex)
{
    const int count = faceCounts[fIndex];
    for (int vIndex = 0; vIndex < count; ++vIndex)
    {
        uint32_t v0 = indices[offset + vIndex];
        uint32_t v1 = indices[offset + ((vIndex + 1) % count)];
        Edge edge(v0, v1);

        halfEdges[offset + vIndex] = { edge, uint32_t(offset + vIndex) };
    }
    offset += count;
}

// Sort to coalesce consecutive entries of the same edge. Use stable sort to
// preserve the order in which edges were introduced. This allows us to remove
// duplicates and keep the first entries of such edges. This first sort step is the
// expensive portion as the list of edged with duplicates is typically twice as long.
std::stable_sort(halfEdges.begin(), halfEdges.end(),
    [](const EdgeEntry& a, const EdgeEntry& b) { return a.e.code < b.e.code; });

// Remove duplicates
halfEdges.erase(std::unique(halfEdges.begin(), halfEdges.end()), halfEdges.end());

// Sort again by entry index, this produces a list of edges in the topology visiting order,
// which is what we want.
std::sort(halfEdges.begin(), halfEdges.end(),
    [](const EdgeEntry& a, const EdgeEntry& b) { return a.entryIndex < b.entryIndex; });

std::vector edges; //< The result
edges.resize(halfEdges.size());
for (size_t i = 0; i < halfEdges.size(); ++i)
    geo.edges[i] = halfEdges[i].e;

In this algorithm we initially produce a vector of Edges including any duplicate. We know the size of the vector upfront is exactly the number of indices (same as the number of half edges). We don’t store just edges, but edges along with the half edge index, which at this point takes the form of the sequence [0, n-1]. Such vector is allocated once, reducing the cost of memory allocation, however we know we allocate twice as many entries on average; hence this implementation will produce a higher memory utilization peak.

In the second block we use stable sort to group together duplicate edges. Stable sort guaranteed that the relative order among identical entries don’t change, so for each duplicate group, the first entry among the duplicates is the first that was encountered by visiting the topology.

In the third block we remove duplicates, so we are left only with unique edges and their insertion index. We sort once more, this time using the insertion index, which produce the final sequence of edges in the desired order.

Because we used sorting algorithms from the library, we have to copy the sorted EdgeEntry vector to the result vector of edges. One could consider implementing a tighter sorting algorithm where edges and inserting indices are separate so that this last copy wouldn’t be required… However, for memory compaction, if we need to keep edges around for numerous meshes, we probably want to consider the EdgeEntry as temporaries to release, while allocate a tight buffer for the exact number of unique edges, just like the code does.

The result of this version of the algorithm is identical to the original, but it runs about 50% faster when compared to a std::unordererd_map.

Minor Valence

Both algorithms so far performed a search/sort over the entire mesh, which is the naive take. But we know that points are connected to one another in a specific way. The point edge valence tends to be a small number, typically not much bigger than 3-4. If we find a way to leverage this property, we can design an algorithm with O(n) complexity, which will perform much better.

It is not possible to know upfront the vertex valence. Let’s enumerate what each point in the topology connects to:

// point index -> connected to
point 0 -> [1, 3]
point 1 -> [0, 2, 4]
point 2 -> [1, 5]
point 3 -> [0, 4]
point 4 -> [1, 3, 5]
point 5 -> [2, 4]

Apparently, to store the full vertex valence we need a vector for 14 entries. Which is twice the number of unique edges we expect, which is not a known quantity. We don’t know how many entries we have until we count them. We don’t know which edge is duplicate until we compare duplicates.

We can create a variant of the valence by reasoning with sorted indices pairs where we only need to observe the smallest of the two indices. Let’s call this the minor valence. Observe again this quantity:

polygon 0: {0, 1, 4, 3} -> edges {(0, 1), (1, 4), (3, 4), (0, 3)}
polygon 1: {1, 2, 5, 4} -> edges {(1, 2), (2, 5), (4, 5), (1, 4)}

The goal is to identify that the last entry in polygon 1 is duplicate. To do that we only need to know which other indices point 1 is paired with. Extrapolating: we need to know which major indices are paired with each minor index. Index 5 only appears as major index, we don’t need to construct valence for it. Let’s count again:

// minor index -> connected major indices
point 0 -> [1, 3]
point 1 -> [2, 4]
point 2 -> [5]
point 3 -> [4]
point 4 -> [5]
point 5 -> []

We only have 7 entries now, half as many as before. We can liberally allocate our scratch space with the size of the number of indices and our algorithm will stay well within. We also have even a smaller number of entries for each minor index where to search for duplicates. This is very good.

Building the minor valence table is done in three steps: first we count how many times the minor index appears in all non-unique edges, we can do that iterating over the topology. We sore those in a vector named valence. In this fist pass we don’t know where to store the major indices yet

std::vector valence(geo.points.size(), 0);

// Count how many times the minor edge index appears. Differently from vertex valence,
// some points will be counted zero times.
for (int faceIndex = 0, faceOffset = 0; faceIndex < nFaces; ++faceIndex)
{
    const int sides      =  faceCounts[faceIndex ];
    const uint32_t* face = &indices   [faceOffset];

    // Here we loop over count-1 and extract the last iteration to avoid the modulo over
    // count
    for (int vertexIndex = 0; vertexIndex < sides; vertexIndex++)
    {
        const Edge edge(face[vertexIndex], face[(vertexIndex+1) % sides]);
        valence[edge.v[0]]++;
    }
    faceOffset += sides;
}

How many major indices we have per point is entirely derived from the topology. A way make space for that in a buffer is to compute the prefix sum of the minor valence. This gives us a table of offsets to the adjacency buffer for each major index list.

// Build a prefix-sum in offsets
std::vector offsets(points .size(), 0);
for (int index = 0, offset = 0; index < points.size(); ++index)
{
    offsets[index] = offset;
    offset += valence[index];
    valence[index] = 0; //< reset valence to use it as a counter for the number of candidates
}

The third step is where we build the final minor adjacency table. However, we actually don’t need that table, we can produce the final edges as we fill the table. The table only serves us to check if we encountered a connected vertex or not. In a nutshell. We have an empty table of adjacency; we have offsets to space where to store major indices. We can iterate again, accumulate candidate major indices to compare against and produce the new edges.

Once more let’s iterate over the topology: every time we construct an edge, we read the for the minor index the offset into the adjacency table and the current number of major indices candidates. In the inner loop we go over the candidates. If we can’t find the major index of the current edge, then we have a new edge! As said, this inner loop is very short, 1.5 elements on average. This makes the algorithm linear in complexity!

// Prepare edges data. An edge is defined as the two indices of the two points forming the edge.
std::vector adjacency(geo.indices.size(), 0);
for (int faceIndex = 0, faceOffset = 0; faceIndex < nFaces; ++faceIndex)
{
    const int sides      =  geo.faceCounts[faceIndex ];
    const uint32_t* face = &geo.indices   [faceOffset];

    // Here we loop over count-1 and extract the last iteration to avoid the modulo over count
    for (int vertexIndex = 0; vertexIndex < sides-1; vertexIndex++)
    {
        const Edge edge(face[vertexIndex], face[(vertexIndex + 1) % sides]);
        const int count  = valence[edge.v[0]];
        const int offset = offsets[edge.v[0]];
        int* candidates = &adjacency[offset];

        bool found = false;
        for (int i = 0; (!found) & (i < count); ++i)
            found = (edge.v[1] == candidates[i]);

        if (!found)                         //< Record a new edge:
        {
            valence[edge.v[0]]++;           //< increment the number of candidates
            candidates[count] = edge.v[1];  //< record the new candidate
            edges.push_back(edge);
        }
    }
    
    faceOffset += sides;
}

Profiling, it stands out that a not negligible portion of the time is spent in imul instruction because of the modulo operations… we are at that level of optimization. Without stepping into micro-optimization territory, I take one last step to unroll the loops over polygons size. Here is the final listing:

// Algorithm based on vertex valence. We leverage the mesh connectivity to eliminate
// the sorting step. The algorithm searches for duplicate edge in the vertex valence,
// which is very small set. This algorithm has O(n) complexity. It is over an order
// of magnitude faster than the original.
std::vector valence(points.size(), 0);

// Count how many times the smallest point index in a edge appears. Some point indices
// will be counted zero times (as we only count the smallest of the two indices), which
// is a bit different from a vertex valence count.
for (int faceIndex = 0, faceOffset = 0; faceIndex < nFaces; ++faceIndex)
{
    const int sides      =  faceCounts[faceIndex ];
    const uint32_t* face = &indices   [faceOffset];

    // Here we loop over count-1 and extract the last iteration to avoid the modulo over
    // count
    for (int vertexIndex = 0; vertexIndex < (sides-1); vertexIndex++)
    {
        const Edge edge(face[vertexIndex], face[vertexIndex+1]);
        valence[edge.v[0]]++;
    }
    // Extracted last entry to avoid modulo over count
    {
        const Edge edge(face[sides-1], face[0]);
        valence[edge.v[0]]++;
    }
    faceOffset += sides;
}

// Build a prefix-sum in offsets
std::vector offsets(points.size(), 0);
for (int index = 0, offset = 0; index < points.size(); ++index)
{
    offsets[index] = offset;
    offset += valence[index];
    valence[index] = 0; //< reset valence to use it as a counter for the number of candidates
}

// Prepare edges data. An edge is defined as the two indices of the two points forming the edge.
std::vector adjacency(indices.size(), 0);
for (int faceIndex = 0, faceOffset = 0; faceIndex < nFaces; ++faceIndex)
{
    const int sides      =  faceCounts[faceIndex ];
    const uint32_t* face = &indices   [faceOffset];

    // The inner loop body is implemented as a lambda to not repeat it twice (to avoid the
    // modulo as we did in the for loop to produce the initial valence). This is the fastest
    // option I have testes so far.
    auto edgeIteration = [&](int vertexIndex, int vertexIndexNext)
    {
        const Edge edge(face[vertexIndex], face[vertexIndexNext]);
        const int count  = valence[edge.v[0]];
        const int offset = offsets[edge.v[0]];
        int* candidates  = &adjacency[offset];

        bool found = false;
        for (int i = 0; (!found) & (i < count); ++i)
            found = (edge.v[1] == candidates[i]);

        if (!found)                            //< Record a new edge:
        {
            valence   [edge.v[0]]++;           //< increment the number of candidates
            candidates[count    ] = edge.v[1]; //< record the new candidate
            geo.edges.push_back(edge);
        }
    };

    // Here we loop over count-1 and extract the last iteration to avoid the modulo over count
    for (int vertexIndex = 0; vertexIndex < (sides-1); vertexIndex++)
    {
        edgeIteration(vertexIndex, vertexIndex+1);
    }
    {
        edgeIteration(sides-1, 0);
    }
    
    faceOffset += sides;
}

Conclusions

I presented three different algorithms that produces the same deterministic result. An interesting new challenge would be to implement a parallel variant, perhaps to run on the GPU. The first two passes of the third algorithm “minor valence” algorithm are trivially parallel, but the last part isn’t unless one can discount the requirement of deterministic edges order in the output.

To the best of my knowledge, the minor valence algorithm I have introduced is novel. But as things go, it is not unlikely there are many GFX developers out there who implemented a similar algorithm. Likely I may have overlooked the literature too, it wouldn’t be the first time I think I invented something new, only to discover after the fact somebody else got there first. If you know this algorithm already has a name, please let me know and I will add a link to it.

[Correction] I have been contacted by other developers who created similar solutions to my minor valence algorithm, but no papers seem to have been published. Perhaps most notably:

• Don Williamson posted his implementation back in 2020. According to him, precursors of his implementation have been used for 20 years, starting from games such as Advent Shadow, Splinter Cell, Fable, and more.
• Aytek Aman & Simon Fenney resorted on a technique similar to my flat vector + sort algorithm on their Fast Triangle Pairing for Raytracing HPG publication.
• Stéphane Ginier described something in between traditional vertex valence with a fixed max count, supported by a fall back to a hash map for points exceeding the max valence.

If I’ll have the time, I will come back to this post to add a table of measurements in respect to some known mesh data. For now, I just say I observed speed improvements over an order of magnitude of my last solution compared to the initial map-based algorithm with std::unordered_map. I didn’t try to compare variants leveraging fast hash maps, mostly because I prefer specific minimalist algorithms over complex containers. If you try let me know that too!