Project: QEM Implementation in Rust

Sunnickel | 15.11.2025

---

Overview

I wanted to explore whether it was possible to automatically create multiple LOD (Level of Detail) versions of a 3D model.
Turns out, it is! By using QEM (Quadric Error Metric), we can determine which vertices of a mesh could have optimized positions for better simplification.

This project implements QEM in Rust to progressively simplify meshes while minimizing geometric error.

If you want to read more yourself: Carnegie Mellon University - Surface Simplification Using Quadric Error Metrics

---

How It Works

The process involves multiple steps to calculate metrics for vertex positions and identify the best candidates for simplification:

1. Plane Equation
2. Quadric Matrix
3. Edge Collapses

---

1. Plane Equation

Determines the plane for each triangle and calculates the corresponding equation.

For the plane equation, we take 3 vertices. These should not be random ones, but vertices that together form a triangle:

• v0 acts as the origin.
• v1 and v2 are distinct points.

We calculate the edges of the triangle as vectors from the origin to the other two points, then cross them and normalize the result to get the plane normal. Finally, we calculate d to complete the plane equation.

let edge1 = v1.sub(&v0);
let edge2 = v2.sub(&v0);
let normal = edge1.cross(edge2).normalize();

let a = normal.x;
let b = normal.y;
let c = normal.z;
let d = -(a * v0.x + b * v0.y + c * v0.z);

(a, b, c, d)

---

2. Quadric Matrix

Builds a matrix for each vertex that encodes how far it is from the planes of its neighboring triangles.

From the plane we just calculated, we construct a quadric matrix, which will be used frequently later.

The matrix is calculated as the outer product of the plane vector with itself:

Q = [a, b, c, d]ᵀ × [a, b, c, d]

Implemented in Rust:

let matrix: [[f32; 4]; 4] = [
    [a * a, a * b, a * c, a * d],
    [b * a, b * b, b * c, b * d],
    [c * a, c * b, c * c, c * d],
    [d * a, d * b, d * c, d * d],
];

---

3. Edge Collapses

Iteratively selects and collapses edges with the lowest error to reduce mesh complexity.

Now we calculate the edge collapses, which is the most complex part of QEM.

For each triangle, we extract all three edges:

let edges: [(Vector3D, Vector3D); 3] = [
    (triangle.0.min(triangle.1), triangle.0.max(triangle.1)),
    (triangle.1.min(triangle.2), triangle.1.max(triangle.2)),
    (triangle.2.min(triangle.0), triangle.2.max(triangle.0)),
];

We then loop through the edges and sum the quadric matrices of the two vertices:

for i in 0..4 {
    for j in 0..4 {
        result_matrix[i][j] = matrix_0.unwrap()[i][j] + matrix_1.unwrap()[i][j];
    }
}

Next, we calculate the determinant of the 3x3 submatrix:

let det = q11 * (q22 * q33 - q23 * q23) - 
          q12 * (q12 * q33 - q23 * q13) + 
          q13 * (q12 * q23 - q22 * q13);

If the matrix is invertible

We calculate the optimal vertex position using Cramer’s rule:

let x = -(  q14 * (q22 * q33 - q23 * q23) - 
            q12 * (q24 * q33 - q23 * q34) + 
            q13 * (q24 * q23 - q22 * q34)
        ) / det;

let y = -(  q11 * (q24 * q33 - q23 * q34) - 
            q14 * (q12 * q33 - q23 * q13) + 
            q13 * (q12 * q34 - q24 * q13)
        ) / det;

let z = -(  q11 * (q22 * q34 - q24 * q23) - 
            q12 * (q12 * q34 - q24 * q13) + 
            q14 * (q12 * q23 - q22 * q13)
        ) / det;

And compute the error at this position:

let error = matrix[0][0] * vector.x * vector.x + 
            matrix[1][1] * vector.y * vector.y +
            matrix[2][2] * vector.z * vector.z +
            matrix[3][3] +
            2.0 * (
                matrix[1][0] * vector.x * vector.y + 
                matrix[0][1] * vector.x * vector.z + 
                matrix[0][2] * vector.x + 
                matrix[1][2] * vector.y * vector.z + 
                matrix[1][3] * vector.y + 
                matrix[2][3] * vector.z
            );

If the matrix is not invertible

We fall back to candidate positions: v0, v1, and the midpoint:

let midpoint = Vector3D::new(
    (v0.x + v1.x) / 2.0,
    (v0.y + v1.y) / 2.0,
    (v0.z + v1.z) / 2.0
);

let candidates = [
    (v0, get_candidate_error(&result_matrix, v0)),
    (v1, get_candidate_error(&result_matrix, v1)),
    (midpoint, get_candidate_error(&result_matrix, midpoint)),
];

*candidates.iter()
    .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal))
    .unwrap()

We then choose the candidate with the lowest error as the new vertex position.

---

AI Assistance

This README.md file has been created with the assistance of AI.

Specifically, AI was used for:

Formatting for this README.md All content has been reviewed and edited to ensure accuracy and clarity. While AI tools were utilized to accelerate the writing process, all information has been manually checked for correctness.