Among the various mesh simplification techniques available, edge collapse is most widely adopted in practice. This strategy is implemented in many major graphics libraries and tools like CGAL , QSlim, and meshoptimizer. An edge collapse operation merges the two endpoints of an edge into a single new vertex, effectively removing the edge and the two triangles that shared it. Repeating this operation iteratively leads to a simplified mesh that maintains the overall structure of the original. Cost functions help determine which edge to collapse and where to place the resulting vertex in order to best preserve the model’s visual and geometric details. While mesh simplification is a well-studied topic, newcomers to the field often face a steep learning curve when engaging with foundational papers. Many of these works emphasize final equations or high-level algorithmic descriptions, offering little insight into the underlying geometric reasoning or the practical implementation details. As a result, readers may struggle to build an intuitive understanding of how and why edge collapse-based simplification works, or how to translate theory into working code.

This paper aims to bridge that gap by offering a detailed, implementation-aware analysis of edge collapse-based mesh simplification on a manifold mesh. Our contributions are as follows:

• We present a complete, end-to-end simplification pipeline that includes well-chosen data structures for representing mesh connectivity, deep analysis of cost functions presented in foundational papers in this space, and the edge collapse algorithm that binds both of these.

• Unlike many prior works that present only the final cost metrics or optimization functions, we derive and explain them along with the geometric meaning behind these formulations, allowing readers to understand the rationale behind each step.

• Our goal is two-fold: to serve as a conceptual guide for learners who want to understand the inner workings of simplification algorithms, and to act as a practical reference for developers looking to implement their own systems.

In this paper, we first categorize and review different families of mesh simplification algorithms. Since the efficiency of edge-collapse operations depends on fast access to mesh connectivity and rapid local updates, int the following section, we discusse data structures that can be employed to store and manage the mesh connectivity information. We then present a comprehensive edge-collapse algorithm, including detailed programmatic checks to prevent mesh degeneracies. Our most extensive section examines cost computation strategies, explaining the mathematical formulations from foundational papers alongside practical implementations. In the next section, we cover advance edge collapse techinques that account for per-vertex attributes. Finally, we provide supplemental mathematical results and proofs that support these techniques.

Mesh simplification techniques vary widely, but most can be grouped by the strategy they use to reduce geometric complexity while maintaining topology as presented in [Cignoni et al. 1998]. We have supplemented this list with recent advances in the field that leverage modern techniques such as machine learning and neural networks.

An early strategy for mesh decimation focused on detecting coplanar or nearly coplanar surface patches and merging them into larger polygonal regions as presented in [De-Haemer Jr and Zyda 1991] and [Hinker and Hansen 1993]. These regions are subsequently re-triangulated to produce a mesh with fewer faces. Despite its simplicity, the method often degraded geometric detail and introduced topological inconsistencies.

Another method, known as vertex clustering, groups nearby vertices based on spatial proximity and replaces each cluster with a single representative vertex, followed by local re-triangulation as presented in [Rossignac and Borrel 1993] and improved in [Low and Tan 1997]. While faster, this method was again found to compromise detail and topological accuracy.

A more refined and topology-sensitive method is iterative local decimation, which incrementally removes vertices, edges, or faces based on localized geometric evaluations.

These operations are typically guided by cost functions designed to preserve the mesh’s overall structure and appearance [Garland and Heckbert 1997;Lindstrom and Turk 1998;Schroeder et al. 1992]. Extensions such as simplification envelopes [Cohen et al. 1996] presents bounded error control by forcing the resulting simplified mesh to lie between two offset meshes.

In energy-based optimization methods, such as the one presented in [Hoppe et al. 1993], a global cost function evaluates the overall quality of the mesh. Simplification is carried out through iterative edge-based operations such as collapse, swap, or split that aims at minimizing both the local and global cost function. Although this approach with global optimization promises a better overall structural preservation, it is less commonly used in practice due to its computational complexity.

A different strategy is retiling, introduced in [Turk 1992], which begins by randomly placing a user-defined reduced number of new vertices on the original surface which are then adjusted based on areas of high curvature. A new reduced triangulation is built on this vertex set. Although effective in reducing triangle count, this method lacks support for per-vertex attributes, making it less suitable for applications like computer-aided design (CAD) or physical simulations where such data is essential.

Another notable approach to mesh simplification is voxelization, as used by works such as [He et al. 1995] and [He et al. 1996]. Here, the mesh is first sampled into a voxel grid, and a low-pass filter is applied at each grid point to generate a discrete scalar field. A triangulated surface is then extracted from this field using the standard marching cubes algorithm or an adaptive variant of it at an isovalue dictated by the filter. The detail of voxel-based meshes can be adjusted via resolution, but the method sees limited industrial use. It smooths sharp features making itself unsuitable for CAD and is computationally expensive due to volumetric processing, and lacks explicit geometric error control, making output quality difficult to guarantee.

Recent work has explored neural methods that either simplify meshes directly or offer implicit representations enabling level-of-detail control. [Potamias et al. 2022] employs a differentiable neural network to select a subset of input vertices using a sparse attention mechanism and re-triangulate the selected vertices, producing simplified meshes in a data-driven, generalizable manner without per-mesh retraining. [Chen et al. 2023] generates a coarse base mesh using QEM, followed by neural remeshing through face splits. A per-face latent feature representation is transmitted and decoded on the client-side to reconstruct finer meshes. This approach implicitly generates simplified representations across multiple LODs. [Park et al. 2019] learns a signed distance field (SDF) representation from a voxelized representation of mesh. [Takikawa et al. 2021] extends it by creating multiscale SDFs giving real-time rendering at various LODs via ray marching. Although simplified triangle meshes can be extracted using methods like marching cubes, this undermines the efficiency of its implicit representation.

Mesh connectivity data structures are designed to efficiently organize and manage the relationships between elements of a mesh such as which faces share an edge, which edges are connected to a vertex, or which vertices make up a face. They allow algorithms to rapidly traverse and manipulate the mesh’s topology. Below are two data structures commonly used to represent mesh connectivity, along with an evaluation of their suitability for supporting edge collapse operations. The Corner Table data structure introduced in [Rossignac 2002] is a compact mesh representation where each triangle’s three “corners” (vertex-triangle associations) are stored in a list. For edge collapse, it efficiently manages the edge-collapse updates and supports fast querying on the mesh. The Half-Edge data structure presented in [McGuire 2000], is widely used due to its intuitive design and broad support across mesh libraries. In this structure, each mesh edge is represented by a pair of half-edges pointing in opposite directions, each storing connectivity to associated elements such as vertices, faces, and neighboring edges. While not the most memory-efficient option, it enables fast mesh queries and local updates, making it ideal for operations like edge collapse. In the code listing below, we present the interfaces that a typical connectivity data structure would support. The queries listed in table. 1 are necessary for the edge collapse-based mesh simplification algorithm, so they must be handled efficiently by the chosen mesh connectivity data structure. As Table 1 illustrates, both the Half-Edge and Corner Table structures are adept at handling these queries with optimal time complexities, making them well-suited for edge-collapse based mesh simplification.

Time complexity (half-edge / corner table )

Get all triangles connected to vertex 𝑣 𝑂(degree(𝑣))

Get all edges connected to vertex 𝑣 𝑂(degree(𝑣))

Get all vertices connected to vertex 𝑣 𝑂(degree(𝑣))

Get all edges connected to vertex 𝑣 𝑂(degree(𝑣))

Table 1. Mesh query operations and their time complexities using different data structures

Edge collapse-based simplification iteratively reduces the number of triangles in a mesh while preserving its overall shape and features. The core algorithm remains largely consistent across different implementations, with key differences lying in the cost metric and vertex placement strategies. The algorithm typically involves the following steps:

Cost assignment and optimal vertex placement calculation 1. A cost is computed for each edge in the mesh to estimate the geometric error introduced by collapsing it. Simultaneously, the optimal position for the resulting merged vertex is determined. This step is critical, as it is where most edge collapse based simplification strategies diverge.

2. The computed cost, along with the edge and its optimal replacement vertex, is stored in a priority queue.

A target triangle count is either defined internally by the program or specified externally by the client code. Then, the following steps are repeated until the target triangle count is reached:

1. Select the edge with the lowest collapse cost from the priority queue.

2. Perform validity checks to ensure that collapsing that edge preserves the mesh’s manifoldness. (The three validity checks we employ are explained below.)

3. If the edge passes all validity checks, collapse it by replacing the edge with the computed vertex and removing the two adjacent triangles.

4. Since the collapse locally alters the mesh, recompute the costs of all the edges connected to the collapsed edge, and update the corresponding entries in the priority queue to maintain accuracy for the next iteration.

5. Update the mesh’s connectivity data structure to reflect the changes made by the collapse.

Checks 1 and 2 follow the criteria established in [Hoppe et al. 1993], while check 3 is derived empirically. These checks are crucial for avoiding degeneracies that may result in invalid or non-manifold mesh structures. 2. Two-neighbor connectivity check: Verify that exactly one pair of edges is merged on each side of the collapsing edge. This condition holds when the two collapsing vertices share exactly two common neighbors. A connectivity-related nonmanifold triangle formation is illustrated in Figure 4.

In edge collapse-based mesh simplification, an error metric is assigned to each edge that estimates the cost of collapsing it. Edges with the lowest error are prioritized for collapse. Additionally, we need effective strategies to determine the best new vertex position that will replace the collapsed edge while minimizing the geometric distortion.

The IConstraint class below defines an interface for error metrics. The cost function classes implementing this interface compute the cost ℇ(𝑣) of collapsing an edge for a candidate vertex position 𝑣. Implementations of this class compute the cost as ℇ(𝑣) = 𝑣 𝑇 𝐻𝑣 + 2𝑐 𝑇 𝑣 + 𝑘, and store the entities {𝐻, 𝑐, 𝑘} in this equation in m_H, m_c, m_k.

These will be used to obtain the optimal vertex placement as well, as detailed in section 7.

This method, described in [Garland and Heckbert 1997], defines error as the sum of distances from the new vertex to the planes of surrounding triangles, treating each vertex as their intersection. This captures how much the new vertex deviates from the original geometry, reflecting the introduced distortion.

Referring to Figure 6, we collapse the edge (𝑣 1 , 𝑣 2 ) into a new vertex 𝑣, removing the two adjacent triangles(in planes 𝑃 1 and 𝑃 6 ) and forming a local geometric approximation. The error is measured as the sum of distances from 𝑣 to the original surrounding planes 𝑃 1 through 𝑃 10 . Let ℙ = {𝑃 1 , 𝑃 2 , … , 𝑃 𝑚 }, the set of planes that surround the edge being collapsed. Let ℇ be the error introduced by the newly added vertex 𝑣, given by:

where (𝑛, 𝑑) represent the unit normal 𝑛 and scalar 𝑑 in the equation 𝑛 • 𝑟 + 𝑑 = 0 of each plane 𝑃 ∈ ℙ.

The term 𝑛 𝑇 𝑣 + 𝑑 gives the signed distance from the vertex to the plane. Squaring it ensures that the error is always non-negative, penalizing both positive and negative deviations equally.

Expanding,

To minimize the error, we set its gradient to zero and solve for 𝑣:

Note: When ℇ takes this form, the matrix 𝐻 is its Hessian matrix. In this specific case, 𝐻 turns out to be positive semidefinite. So, the point that makes ∇ℇ = 0 corresponds to a minimum point rather than a maximum or a saddle point.

If the matrix 𝐻 is non-invertible (i.e., det(𝐻) = 0), the optimal vertex position cannot be computed this way. In such cases, fallback strategies or alternative constraints are used. For this cost function, a non-invertible 𝐻 indicates that the surface surrounding the edge collapse is flat, as explained below:

in block notation. This means that each term of the sum is itself a non-invertible matrix, as all its columns are parallel to 𝑛. So, when 𝐻 is non-invertible, all the normals of the planes forming 𝐻 are parallel. This occurs if the local surface is flat.

The “quadric” in the name of this method is derived from the form this error takes when 𝑣 is represented in homogeneous 4-dimensional coordinates as the vector ( 𝑣 1 ). In that case, the error ℇ is expressed as follows:

where the authors define the 4x4 matrix 𝑄 = ( 𝐻 𝑐 𝑐 𝑇 𝑘 ) above as the total error quadric for this edge. It can further be decomposed as the sum of fundamental error quadrics 𝐾 𝑃 for each plane 𝑃 ∈ ℙ:

However, the same authors in [Garland and Heckbert 1998] found this formulation impractical because it requires computationally expensive matrix operations on higherdimensional matrices, like inversion. For this reason, it won’t be discussed further.

The standard QEM method struggles with boundary edges -those with only one adjacent face. As noted in [Garland and Heckbert 1998], a modified QEM was proposed to address this and preserve boundary edges. Consider the red boundary edge between 𝑣 1 and 𝑣 2 , selected for collapse under two distinct surrounding geometries. The new vertex 𝑣 is computed as the intersection of the adjacent planes 𝑝 1 , 𝑝 2 and 𝑝 3 because the distance of that point from all these planes is zero -resulting in the minimum possible quadric error. Depending on their configuration, this intersection may lie above or below the original boundary.

In conventional QEM, no explicit constraint anchors the new vertex to the boundary. Consequently, collapsing a boundary edge tends to displace the vertex away from the boundary, a deviation that compounds as more boundary edges are collapsed. This progressive drift results in noticeable degradation of mesh quality, as seen in Figure 8. To counteract this effect, [Garland and Heckbert 1998] introduces an imaginary plane 𝑝 ′ as shown in Figure 9 in addition to the actual planes adjacent to the collapsing edge. 𝑝 ′ is defined as the plane perpendicular to the mesh plane containing edge (𝑣 1 , 𝑣 2 ). A new term, 𝑑 2 (𝑣, 𝑝 ′ ), representing the squared distance between the vertex and 𝑝 ′ , is incorporated into the error metric. As 𝑣 moves away from the boundary, this term increases, exerting a corrective pull toward the boundary. To strengthen this constraint, the quadric for 𝑝 ′ is scaled by a large constant before being added to the quadrics of the edge endpoints. The method is further extended to treat edges separating faces with different attribute values (e.g., material indices) as boundaries. This ensures that such attribute boundaries are preserved during simplification, concentrating edges and faces along these divisions for improved alignment. An example of this extension is shown in Figure 10.

This constraint, introduced in [Lindstrom and Turk 1998] helps preserve the mesh volume. If the new vertex replacing the collapsed edge isn’t chosen carefully, it can distort the model. For instance, using the edge midpoint as the new vertex might increase the volume in concave areas or decrease it in convex ones. The goal of this constraint is to preserve volume locally at each collapse, thereby minimizing the overall volume change across the whole model.

Neither boundary nor volume preservation guarantee geometric integrity; boundaries may deform, and surfaces can lose detail. However, these constraints serve as useful heuristics. Preserving simple, quantifiable properties like area and volume helps reduce extreme distortions, even if local features like sharp edges or curves are lost. While these constraints don’t capture fine geometric details, they provide an efficient way to maintain overall structure, balancing accuracy and performance without the complexity of exact boundary or volume preservation.

When an edge 𝑒 is collapsed, it sweeps out a tetrahedral volume as illustrated in Let 𝑡 = [𝑣 1 , 𝑣 2 , 𝑣 3 ], 𝑡 ′ = [𝑣, 𝑣 2 , 𝑣 3 ], and the volume swept by 𝑡 as 𝑣 1 moves linearly to 𝑣 be 𝑉(𝑣, 𝑣 1 , 𝑣 2 , 𝑣 3 ). 𝑉 is positive if 𝑣 is above the plane of 𝑡 and negative otherwise.

Thus, to preserve the local volume at the site of an edge collapse, the sum of volumes of tetrahedra swept with all triangles 𝑇 = {𝑡 1 , 𝑡 2 , …, 𝑡 𝑛 } connected to edge 𝑒 are considered. The change in volume is given by (the superscript 𝑡 indicates the vertices belonging to triangle 𝑡):

Solving for ℇ = 0 and expanding the determinant along the fourth row, we get,

Representing the determinants that include 𝑣 as scalar triple products, we get,

Simplifying the term

, we get:

where 𝑛 𝑡 is the normal of the plane containing

with magnitude equal to triangle area.

Substituting the above simplified term in Equation 1, we get,

The above equation is of the form 𝑣 • 𝑁 = 𝐷 which defines a plane. This implies that the vector 𝑣 is restricted to lie on a plane. So, any point on that plane will satisfy the equation above, implying that volume preservation alone is not enough to fully determine 𝑣: we need 2 other constraint equations to do so. By contrast, volume preservation only ensures that the total volume added and removed balances out to zero. And, as we know, it forces the vertex to lie on a plane but doesn’t tell us exactly where on that plane to place it, so it leaves some freedom. Moreover, it can lead to local distortions if large volumes are added and subtracted in different areas.

From the volume preservation constraint formulation, we know that the change of volume induced by an edge collapse is:

where:

• 𝑉(𝑣, 𝑣 1 , 𝑣 2 , 𝑣 3 ) is the volume swept out by 𝑡 when 𝑣 1 moves in a linear path to 𝑣 If the vertex 𝑣 is above the plane of a triangle 𝑡, the signed volume 𝑉 of the tetrahedron is positive. If below, it’s negative. But for optimization, we care about how much the volume changes, not the direction. So, we use the unsigned volume change. To get unsigned volume, we could use |𝑉| or 𝑉 2 . We use 𝑉 2 because it is differentiable everywhere. This matters because optimization algorithms rely on gradients and |𝑉| has a kink at zero where the gradient is undefined.

So we express ℇ as the sum of squares of volumes instead. So we get,

and we get,

This constraint uniquely determines 𝑣, except in degenerate cases where det(𝐻) = 0. Just like the case of QEM, this happens in locally flat regions of the geometry, since we know that 𝐻 is defined as:

which has the form used in QEM. So, 𝐻 becomes non-invertible in flat regions, where all 𝑛 𝑡 are parallel and the sum reduces to scaled rank-1 terms. In such situations, alternative constraints or fallback strategies are required for vertex placement.

This constraint is discussed in [Lindstrom and Turk 1998]. It helps compute optimal vertex placement and edge collapse error by preserving the area of boundaries. In Figure 12, the image on the left shows a boundary edge (in red) on a planar hole, while the image on the right shows it being replaced by a new vertex (red dot) after edge collapse. As a result, although the total shaded area is preserved (red area loss offset by blue area gain), the boundary’s shape and structure are visibly altered. Thus, the constraint preserves area, not the boundary itself.

As per Figure 13, let edge (𝑣 1 , 𝑣 2 ) be collapsed into vertex 𝑣 and let ℇ be the net area change. Collapsing a boundary edge connects the new vertex 𝑣 to two other boundary Collapsing a boundary edge connects two others, forming three triangles (𝑛 = 3) as shown in Figure 14. The squared norm is used for computational simplicity. The error is then defined as the squared magnitude of the total area change induced by vertex 𝑣 as:

as 𝐸 1 and

as 𝐸 2 . This gives,

To simplify the cross-product term, 𝐸 1 × 𝑣 can be written as 𝒮𝑣, where 𝒮 is the skewsymmetric matrix for the cross product with 𝐸 1 :

Simplifying the squared norm gives:

which has the same form as in QEM, suggesting that the solution should be the same: 𝑣 = -𝐻 -1 𝑐.

However, for this constraint, 𝐻 is non-invertible as 𝒮 is non-invertible, being a skewsymmetric matrix. Therefore, 𝑣 cannot be fully determined. Below is a demonstration of the constraints that can actually be extracted from setting the gradient to zero.

= 𝐸 1 × 𝐸 2 -𝐸 1 (𝐸 1 • 𝑣) + 𝑣(𝐸 1 • 𝐸 1 ) … using vector triple product

This is a 3D vector equation in 𝑣, which can be split into 3 scalar equations to solve for its components. We can choose any basis for this, but for convenience, we choose:

Projecting Equation 2 onto 𝐾 gives us:

Projecting Equation 2 onto 𝐸 1 gives us:

This simplifies to 0 = 0, which is a degenerate result. This indicates that the component of 𝑣 parallel to 𝐸 1 is not determined by this minimization problem, as it does not affect the value of the error function.

Projecting Equation 2 onto 𝐸 1 × 𝐾 gives us:

Thus, the solution space of this optimization lies in the intersection of two planes defined by:

Even with this approach of minimizing ℇ, 𝑣 remains undetermined. Thus, additional constraints are needed alongside the two equations to solve for 𝑣.

The boundary optimization constraint introduced in [Lindstrom and Turk 1998] minimizes boundary triangle area like boundary preservation, but focuses on unsigned area. The error is expressed as a sum of squared signed areas, giving:

where,

• 𝑣 ∈ ℝ 3 Expanding Equation 3gives:

, gives:

To simplify the cross-product term, 𝑣 × 𝑒 1 can be written as 𝒮𝑣, where 𝒮 is the skewsymmetric matrix for the cross product with 𝑒 1 . Also, the scalar triple product identity can be used to rearrange 2(𝑣 × 𝑒 1 ) 𝑇 𝑒 2 = 2(𝑒 1 × 𝑒 2 ) 𝑇 𝑣.

which is of the same form as earlier constraints, so we obtain 𝑣 as 𝑣 = -𝐻 -1 𝑐.

As in earlier constraints, if det(𝐻) = 0, the constraint becomes degenerate and we need to use other constraints or use fallback strategies for vertex placement.

Triangle shape optimization tries to improve triangle quality. Skinny or stretched triangles can cause shading issues, while more even, equilateral triangles make the mesh look and work better. As shown in Figure 15, the vertex placement on the left side leads to a cleaner triangle structure thanks to its more regular, evenly shaped triangles. The placement on the right side contains long, stretched triangles, which make the mesh look less tidy and visually less appealing.

Before analyzing this further, an important triangle shape quality metric needs to be introduced: the area-to-perimeter ratio. Regular triangles (like equilateral ones) have a higher area-to-perimeter ratio, which means they’re more compact and less stretched.

In the triangle shape preservation constraint, the goal is to maximize the area-toperimeter ratio. We make an assumption that the region around the collapsing edge is nearly flat. This means the total area of the nearby triangles doesn’t change much after the edge collapse. So, to improve their area-to-perimeter ratio, we can focus on reducing the perimeter alone, which is determined by the edge lengths.

That’s why we minimize the sum of the squared lengths of the edges connected to the new vertex. This pulls the vertex into a position where the edges are more evenly distributed and shorter, leading to more balanced, less skinny triangles.

We formulate this error ℇ as:

where 𝑣 𝑖 refers to each neighboring vertex {𝑣 1 , … , 𝑣 𝑛 } connected to 𝑣 in the mesh.

On expanding this equation for ℇ, we get,

𝐼 is an identity matrix = 𝑣 𝑇 𝐻𝑣 + 2𝑐 𝑇 𝑣 + 𝑘 which can be solved in the same way as in other constraints. Here, it leads to this solution:

Thus, optimizing triangle shape will lead to choosing the centroid of the neighboring vertices.

Skinny triangles are avoided because they cause shading artifacts. Rasterization interpolates per-vertex data (e.g., normals) using barycentric coordinates, defined as:

In skinny triangles, the denominator 𝐴𝑟𝑒𝑎(⧍𝐴𝐵𝐶) becomes very small, making the coordinates numerically unstable. Small floating-point errors in the vertex positions or sub-areas can then cause large interpolation errors, producing shading artifacts.

When the matrix used to compute 𝑣 is non-invertible, fallback strategies are needed. A simple and common fallback is to place the new vertex at the midpoint of the edge (𝑣 1 , 𝑣 2 ) being collapsed:

vec3 GetFallbackVertex(IEdge * collapse_edge) { auto verts = collapse_edge->GetVertices(); return 0.5 * (verts[0]->GetPosition() + verts[1]->GetPosition()); }

To solve for the new vertex 𝑣, a system of linear equations is formed from several constraints as discussed in the earlier sections. Each linear equation is of the form:

Here, each 𝑎 represents the normal of a constraint plane, and 𝑏 is the corresponding offset. Geometrically, we are finding the point 𝑣 that lies at the intersection of all these planes in ℝ 3 . In theory, only three linearly independent constraints (planes) are needed to uniquely determine a point in 3D. However, the algorithm includes more than three constraints to ensure robustness. That’s because:

• some constraints may become redundant (linearly dependent).

• some may become degenerate in flat or symmetric regions.

The final vertex placement includes only the best three by checking for linear independence and stability using the following criteria while adding constraints one by one: First Constraint(𝐚 𝟏 ): Is it valid?

return length(A) > 0; } Second Constraint(𝐚 𝟐 ): Is it linearly independent from the first?

This checks if the angle 𝜃 between the first and second constraint normals is sufficiently large, that is, that they are not almost parallel. 𝛼 is a threshold angle used to determine acceptable linear independence.

This ensures that the third constraint’s normal 𝑎 3 does not lie in the plane formed by the normals of the first two constraints, again up to a threshold angle 𝛼. This guarantees that the three planes intersect at a single point in 3D space, defining a unique solution for the new vertex.

Note that, here we use sin(𝛼) and not cos(𝛼) when computing the dot product. We define 𝛼 as the angle between the plane formed by 𝑎 1 and 𝑎 2 and the new vector 𝑎 3 . So as Figure 16 suggests, the angle between the vectors 𝑎 1 × 𝑎 2 and 𝑎 3 will be 90°-𝛼. So we have,

Putting it all together, we revisit the overall edge collapse algorithm and implement all its components. We begin by computing the constraints described above in an order that best suits our use case. Next, we apply the constraint selection criteria to form a set of solvable constraints. Finally, we compute the optimal vertex position based on this set, resorting to our fallback strategy if the system remains unsolvable. The resulting optimal position and its associated cost are then returned for further handling in the priority queue.

void GetCollapseVertex(const IMesh * mesh, const IEdge &collapse_edge , vec3 &collapse_vertex, float &collapse_error) { vector constraints;

If needed, they can be determined using the Gram-Schmidt process, which begins with 𝑛 linearly independent vectors and iteratively removes components parallel to previously computed basis vectors.

Let the plane be determined by three points 𝑝, 𝑞, 𝑟 ∈ ℝ 𝑛 . As shown in Figure 17, two orthogonal axes 𝑒 1 and 𝑒 2 can be defined on the plane spanned by (𝑝, 𝑞, 𝑟):

• 𝑒 1 is the unit vector in the direction 𝑞 -𝑝, i.e., 𝑒 1 = 𝑞 -𝑝 ‖ ‖ ‖ ‖ 𝑞 -𝑝 ‖ ‖ ‖ ‖

• Using the Gram-Schmidt process, 𝑒 2 can be obtained by taking the vector 𝑟 -𝑝, removing its projection onto 𝑒 1 (to ensure orthogonality with 𝑒 1 ), and then normalizing. In other words,

Next, a point 𝑝 lying on the plane is chosen, and the vector from 𝑝 to 𝑣, i.e 𝑤 = 𝑣 -𝑝, is expressed as the sum of its components along the orthonormal basis {𝑒 1 , 𝑒 2 , … , 𝑒 𝑛 }:

Next, the components along the plane i.e (𝑒 1 , 𝑒 2 ) are removed from 𝑤, giving us a vector 𝑢:

Here, 𝑢 is the perpendicular from 𝑣 to the plane, whose squared length is precisely ℇ:

Expanding this out further using 𝑤 = 𝑣 -𝑝, we get:

= 𝑣 𝑇 𝐻𝑣 + 2𝑐 𝑇 𝑣 + 𝑘 Note that this expression matches the form of the error used in the original QEM method. This means that, aside from differences in the values and dimensions of the entities {𝐻, 𝑐, 𝑘}, the procedure for calculating the cost and determining the optimal 𝑣 remains unchanged: the optimal 𝑣 is still given by -𝐻 -1 𝑐.

Furthermore, once computed, 𝑣 not only represents the optimal vertex position but also encodes the optimal values for all associated scalar attributes. [Hoppe 1996] introduced progressive meshes -sequences of meshes representing varying levels of detail of an input mesh, each created through successive edge collapse operations. They propose an alternative method to compute the cost of each edge collapse by defining it as the difference in an energy function. The cost reflects how much the energy function of the mesh changes before and after the collapse, with edges causing smaller energy differences considered better candidates for collapse. Their approach also handles discrete face attributes and scalar vertex attributes at each level of detail.

Before introducing the cost function, we define key geometric entities and setup steps needed to compute the cost of an edge collapse. Figure 18 shows an example.

The original mesh (before any edge collapse operations) is denoted as:

• Vertices: V = {𝑣 1 , … , 𝑣 𝑛 } , where each 𝑣 𝑖 ∈ ℝ 3 These points will serve to approximate M in the energy functions.

For each 𝑥 𝑖 , compute its scalar attribute 𝑥 𝑖 ∈ ℝ 𝑑 via barycentric interpolation on its containing face yielding the attribute set 𝑋 corresponds to 𝑋:

Consider the mesh M after some edge collapses. Let 𝑀 ↑ be the state before collapsing another edge, and 𝑀 the state after edge collapse. The cost for this collapse is formulated as:

Here, the spring constant 𝜅 weights this term.

𝐸 𝑠𝑐𝑎𝑙𝑎𝑟 is similar to 𝐸 𝑑𝑖𝑠𝑡 , but for scalar attributes. It measures the squared difference between the attributes of each sampled point 𝑥 𝑖 and those of its projection 𝑝(𝑥 𝑖 ):

where 𝑐 𝑠𝑐𝑎𝑙𝑎𝑟 weights this term.

𝐷 𝑑𝑖𝑠𝑐 penalizes collapsing a sharp edge 𝑒 that affects discontinuities tracked by 𝑋 ′ :

Here, numProject(𝑋 ′ , 𝑒) denotes the number of points in 𝑋 ′ projecting onto 𝑒. This term is applied only if the collapse alters the discontinuity connectivity, based on criteria presented in [Hoppe 1996]. To forbid the collapses entirely, set 𝐷 𝑑𝑖𝑠𝑐 = ∞.

Recall that 𝑀 ↑ is the state of the mesh before collapsing the current edge, and 𝑀 the mesh after the collapse. The new vertex is placed at position 𝑣 with attributes 𝑣. We need to compute 𝑣 and 𝑣 that minimize the cost ℇ.

As per Equation 4, the terms in ℇ depend on 𝑣 and 𝑣 as follows:

• 𝐸 𝑑𝑖𝑠𝑡 (𝑀) and 𝐸 𝑠𝑝𝑟𝑖𝑛𝑔 (𝑀) depend only on 𝑣.

• 𝐸 𝑠𝑐𝑎𝑙𝑎𝑟 (𝑀) depends only on 𝑣.

• 𝐷 𝑑𝑖𝑠𝑐 and the energy terms for 𝑀 ↑ are independent of both 𝑣 or 𝑣. So, they only affect the cost ℇ but not the optimization.

Thus, here are the steps to compute ℇ, 𝑣 and 𝑣 for a given edge collapse:

1. Minimize ∆𝐸 𝑑𝑖𝑠𝑡 + ∆𝐸 𝑠𝑝𝑟𝑖𝑛𝑔 over 𝑣.

2. Minimize ∆𝐸 𝑠𝑐𝑎𝑙𝑎𝑟 over 𝑣.

3. If the discontinuity criteria hold, compute 𝐷 𝑑𝑖𝑠𝑐 .

4. Return ℇ = ∆𝐸 𝑑𝑖𝑠𝑡 + ∆𝐸 𝑠𝑝𝑟𝑖𝑛𝑔 + ∆𝐸 𝑠𝑐𝑎𝑙𝑎𝑟 + 𝐷 𝑑𝑖𝑠𝑐 along with the optimal 𝑣 and 𝑣.

Minimizing the error ℇ is more complex than previous error functions and requires a staged optimization process. We describe that process in detail below, and the data required for computing and optimizing the error is shown in an example in Figure 19.

We can observe that any edge with unchanged endpoints in 𝑀 and 𝑀 ↑ cancels out in ∆𝐸 𝑠𝑝𝑟𝑖𝑛𝑔 . Since differences occur only near the collapsed edge, only two edge sets contribute nonzero terms:

1. 𝑁: edges in 𝑀 incident on the new vertex at 𝑣.

2. 𝑁 ↑ : edges in 𝑀 ↑ incident on the collapsed edge.

Thus, we can simplify ∆𝐸 𝑠𝑝𝑟𝑖𝑛𝑔 as follows:

where the projection of the same sampled point 𝑥 𝑖 is 𝑝(𝑥 𝑖 ) in 𝑀 and 𝑝 ↑ (𝑥 𝑖 ) in 𝑀 ↑ . Since 𝑀 and 𝑀 ↑ differ only in the vicinity of the one edge being collapsed, most points project identically in both states and cancel out. Thus, only points 𝑌 ⊆ 𝑋 projecting onto the neighborhood of the edge contribute:

We define the neighborhood of an edge as the set of faces connected to either endpoint of that edge in 𝑀 or 𝑀 ↑ . This simplification is based on a locality assumption: the new vertex 𝑣 stays near the collapsed edge, so the projections of distant points in 𝑋 remain unchanged. Although allowing 𝑣 farther away could lower the cost, it would require recomputing over all 𝑋. In practice, restricting 𝑣 and using the simplified ∆𝐸 𝑑𝑖𝑠𝑡 works well.

Since where 𝑝 lies on some face in 𝑀 with vertices (𝑣 𝑎 , 𝑣 𝑏 , 𝑣 𝑐 ). We use barycentric coordinates 𝛽 = (𝛽 𝑎 , 𝛽 𝑏 , 𝛽 𝑐 ) to express 𝑝:

We can define a function 𝛽(𝑦) to denote the projected barycentric coordinates for any point 𝑦: 𝛽(𝑦) = arg min 𝛽 ‖𝑦 -𝑉𝛽‖ 2 , 𝑝(𝑦) = 𝑉𝛽(𝑦)

Now, since both 𝑉 and 𝛽 are face-specific, we extend 𝑉 to the full mesh 𝑀 with 𝑛 vertices {𝑣 1 , … , 𝑣, … , 𝑣 𝑛 } and any given 𝛽 to an 𝑛-dimensional vector, zeroing out the entries for all vertices except (𝑣 𝑎 , 𝑣 𝑏 , 𝑣 𝑐 )

Thus, 𝐸 𝑑𝑖𝑠𝑡 (𝑀) can be rewritten in terms of 𝑣 (as contained within 𝐕):

We can now minimize 𝐸 𝑑𝑖𝑠𝑡 (𝑀) over 𝑣. Since evaluating it requires projecting each 𝑦 𝑖 via an inner minimization in β space, the problem becomes nested: an outer minimization over 𝑣 and inner minimizations over all β(𝑦 𝑖 ).

We solve the nested minimization iteratively as shown in Figure 20, starting with an initial guess for 𝑣 and alternating between: optimizing 𝑣 with fixed β(𝑦 𝑖 ), then updating β(𝑦 𝑖 ) with fixed 𝑣. This repeats until convergence, i.e., when the values of 𝑣 and each β(𝑦 𝑖 ) don’t change much between iterations. In practice, a small number of iterations is sufficient.

The inner minimization -over each β(𝑦 𝑖 ) with fixed 𝑣 -is called the projection subproblem, i.e., projecting all points in 𝑌 onto 𝑀. A brute-force method is to try projecting every 𝑦 𝑖 on every face of 𝑀 and compute β(𝑦 𝑖 ) corresponding to the closest face. But, [Hoppe 1996] adds two speedups to this approach:

1. Use a spatial partitioning structure to find candidate faces in 𝑂(1) time per point, especially useful early on or after edge collapses in new regions.

2. If 𝑦 𝑖 was projected onto 𝑀 ↑ , limit its projection on 𝑀 to the faces neighboring the previous one, leveraging locality.

The outer minimization -over 𝑣 while keeping all β(𝑦 𝑖 ) constant -is now solved by rewriting 𝐸 𝑑𝑖𝑠𝑡 (𝑀) using that constancy.

Here, 𝑌 ⊆ 𝑋 is a local subset of sample attribute vectors. We define this neighborhood as the set of all sample points on the faces adjacent to the edge being collapsed, consistent with the approach used for 𝐸 𝑑𝑖𝑠𝑡 .

Since 𝐸 𝑠𝑐𝑎𝑙𝑎𝑟 ( 𝑀 ↑ ) is independent of 𝑣, it is ignored in the optimization process and added only to the final cost.

The optimization is simplified by reusing the barycentric coordinate sets computation β and β ↑ from the previous optimization of ∆𝐸 𝑑𝑖𝑠𝑡 . 𝑝

In this paper, we investigated mesh simplification using edge collapse in detail by performing a deep dive into four important papers providing variations on this algorithm: [Garland and Heckbert 1997;Lindstrom and Turk 1998;Garland and Heckbert 1998;Hoppe 1996].

We started by discussing the basics of mesh simplification and the different categories of algorithms that are used for that purpose. We then focused on edge collapse and the half-edge data structure typically used to implement it. Next, we outlined the general algorithm used for simplification via edge collapses, including important edge cases that need to be considered. We then performed an elaborate analysis of the process of computing the error introduced by a candidate vertex placement through a variety of metrics, and discussed how the associated constraints can be assembled to form a solvable system of linear equations that yield the final, optimal vertex placement for a given edge collapse.

In the process, we also dealt with other important considerations while performing mesh simplification, such as handling boundary edges and vertex/face attributes.

We believe this work can help people interested in geometry processing and mesh simplification to understand these potent algorithms and metrics in depth, implement them for their use cases, and inspire further work in this field.