Quadric Error Metrics and What an Edge Collapse Actually Costs

Simplification is the inverse of refinement and is harder, because the error introduced by removing geometry has to be estimated without knowing the surface it is being compared against. The quadric error metric answers with a quadratic form: for each candidate vertex position, a single scalar says how far the surrounding planes have been displaced. The estimate is not a bound and not a measurement; it is a linear-algebraic surrogate that happens to correlate well with what the eye notices.

The plane set behind a vertex

A vertex in a triangle mesh has a set of planes associated with it: one per incident face, plus, in the Garland–Heckbert formulation, one per face of the neighbouring ring at some weighting. Each plane is a unit normal $n$ and offset $d$, so the plane is the set of points where $n \cdot x + d = 0$.

The squared distance from a point $x$ to that plane is $(n \cdot x + d)^2$. Expanding in coordinates gives a quadratic form in the four-vector $\tilde{x} = (x, y, z, 1)$:

$$ (n \cdot x + d)^2 = \tilde{x}^\top Q \tilde{x}, \qquad Q = \begin{pmatrix} n n^\top & d n \ d n^\top & d^2 \end{pmatrix}. $$

$Q$ is $4\times4$, symmetric, and has rank one: it is the outer product of the four-vector $(n_x, n_y, n_z, d)$ with itself. Summing the $Q$ matrices over all planes associated with a vertex gives a matrix whose quadratic form is the total squared distance from $x$ to all of those planes. The matrix is still $4\times4$ and still symmetric, so it has 10 distinct entries, which is the storage cost per vertex: 40 bytes in single precision and 80 bytes in double, against 12 bytes for a position.

The reason the sum is the right object rather than an average is that the sum counts planes, so a vertex in a densely tessellated region where many small triangles share a plane-like neighbourhood accumulates a large coefficient along that plane’s normal direction and a small one along the tangential directions. The quadric is an anisotropic measure: it is expensive to move a point off the plane of a flat region and cheap to move it within that plane. That anisotropy is why the method preserves creases and flat regions while flattening curved detail.

Collapsing an edge

An edge collapse takes two vertices $v_1$ and $v_2$ and replaces them with one vertex at position $\bar{v}$, deleting the two triangles that shared the edge and retriangulating the hole implicitly by relinking the remaining neighbours. The new vertex inherits the plane set of both, so its quadric is $Q = Q_1 + Q_2$. That addition is exact in the sense that the resulting quadratic form equals the sum of squared distances to the union of both plane sets, evaluated at the same point. No error is introduced by the merge itself.

The position $\bar{v}$ is chosen to minimize $\bar{v}^\top Q \bar{v}$. Differentiating the quadratic and setting the gradient to zero gives a $3\times3$ linear system from the upper-left block:

$$ \begin{pmatrix} q_{11} & q_{12} & q_{13} \ q_{12} & q_{22} & q_{23} \ q_{13} & q_{23} & q_{33} \end{pmatrix} \bar{v} = -\begin{pmatrix} q_{14} \ q_{24} \ q_{34} \end{pmatrix}. $$

The block is the sum of outer products $n n^\top$, so it is positive semidefinite and singular exactly when the plane normals span fewer than three directions. For a vertex on a flat region all normals are parallel, the matrix has rank one, and the system is singular. The standard fallbacks are to try the two endpoints and the midpoint and take whichever has the smallest quadric value, or to add a small multiple of the identity to the block to make it invertible. The first is more honest, since it reports a value that is actually attained rather than one from an extrapolated position.

import numpy as np

def quadric(n, d):
    v = np.array([n[0], n[1], n[2], d], dtype=float)
    return np.outer(v, v)

def best_position(Q):
    A, b = Q[:3, :3], -Q[:3, 3]
    if abs(np.linalg.det(A)) > 1e-12:
        return np.linalg.solve(A, b)
    # degenerate: all planes parallel, so any point on them ties
    return None

def collapse_cost(Q1, Q2, candidates):
    Q = Q1 + Q2
    sol = best_position(Q)
    if sol is not None:                        # interior optimum, if one exists
        candidates = list(candidates) + [tuple(sol)]
    best = None
    for c in candidates:
        p = np.array([c[0], c[1], c[2], 1.0])
        cost = float(p @ Q @ p)
        if best is None or cost < best[0]:
            best = (cost, tuple(c))
    return best

Every collapse candidate in the mesh is scored this way, the cheapest is performed, and the quadrics and costs of all affected edges are recomputed. A priority queue makes the selection $O(\log n)$ amortized per collapse, and the number of collapses needed to reach a target face count is at most the number of faces removed.

Greedy order and the queue problem

The algorithm is greedy: it takes the locally cheapest collapse each time, which means the cost of an edge is a function of the current mesh, not of the original. After a collapse, every edge incident to the changed vertex has a different quadric and therefore a stale priority. The standard implementation pushes updated entries into a heap and discards stale ones on pop by checking a version counter. Without that check, the heap fills with entries for collapses that no longer exist, and a collapse can be applied twice.

Greedy order also produces a specific failure that the metric does not detect: a collapse can invert a triangle. If the new vertex position lies on the far side of a neighbouring triangle’s plane, that triangle flips orientation. A flipped triangle is worse than a slightly higher error, so the usual guard tests each surviving face’s normal before and after and rejects the collapse if the dot product goes negative. Some implementations also reject collapses that produce a face with an angle below a threshold, which prevents slivers.

Topological validity is a second guard. An edge collapse changes the link of a vertex, and in a closed manifold that can change the genus — collapsing an edge in a thin tube can pinch it closed. The link condition, stated for a manifold triangle mesh, requires that the intersection of the two vertex links be exactly the two opposite vertices of the faces sharing the edge. Enforcing it prevents non-manifold results at the cost of rejecting some geometrically valid collapses.

Boundary and attribute quadrics

Interior surfaces have planes on both sides; a mesh with a boundary has vertices whose incident faces only cover part of the neighbourhood, and the accumulated quadric is rank-deficient in the direction pointing off the boundary. The consequence is that the boundary can drift and shrink during simplification, because moving a boundary vertex along the boundary curve costs nothing. The standard fix is to add constraint planes perpendicular to the boundary: for each boundary edge, a plane containing the edge and perpendicular to the adjacent face. These penalize movement off the boundary curve.

Attributes need their own quadrics. A vertex carrying a colour, a texture coordinate, or a normal lives in a higher-dimensional space, and the appropriate error measure is the displacement in that space, not just in position. The usual extension appends the attribute dimensions to the vector and builds quadrics over them, which grows the matrix to $7\times7$ for position plus a three-component colour and $6\times6$ for position plus a two-component texture coordinate, since the constant term remains. The matrix is still symmetric and still a sum of outer products, so the formulation does not change, but the storage per vertex and the solve cost both grow.

The subtlety is that attribute error and position error are not commensurate. A displacement of one unit in position and a displacement of one unit in a normalized texture coordinate have entirely different visual significance, so the composite quadric needs per-attribute scaling, and the scaling factors are the one part of the method with no principled derivation. They are tuned by whoever builds the pipeline.

Texture coordinates make this worse than the other attributes, because their scale depends on the parameterization and its distortion. A chart that maps a large surface region into a small part of the atlas has a large texel-per-unit-area ratio, so a small movement in UV space corresponds to a large movement on the surface, and a fixed weight for the UV quadric is wrong by that ratio. The honest approach weights each attribute quadric by the inverse of its local scale, which means the weights have to be recomputed when the parameterization changes — that is, after any re-unwrap, not once at asset build time.

What the metric cannot see

The quadric measures displacement of planes, and planes are local. A collapse that removes a feature entirely — a spike, a thin fin, a small hole — can have a low quadric cost if the feature’s planes happen to be nearly coaxial with those of the surrounding surface. The metric is blind to silhouette and to any error measured against a viewpoint.

It is also blind to the difference between an error that is hidden and one that is visible. A collapse inside a crevice and a collapse on a convex silhouette can score identically; the rendering consequence is very different. View-dependent variants exist that weight the quadric by projected screen-space error, and they are what a level-of-detail system uses when it needs a continuous range of detail rather than a single decimated mesh.

Finally, the quadric is computed from the input mesh, so the error it measures is relative to that mesh and to nothing else. A decimated mesh compared against the original reads as accurate if the quadric is small, and that says nothing about how far the original was from the surface it was meant to represent. The plane set encodes the input’s own simplified view of geometry, and simplification compounds approximations rather than replacing them.

The result of decimation is an indexed triangle mesh and feeds directly into the structures described under indexed meshes and the vertex cache; simplification order and vertex ordering interact, since a greedy collapse sequence produces an index stream with poor locality unless it is reordered afterwards.