Decimating Meshes Without Losing Breaklines

The twenty-million-triangle mesh looked superb. The two-million-triangle version delivered to the client has rounded kerbs, a quarry bench edge that curves where it should be sharp, and a building parapet that has become a gentle ramp. The average geometric error is 1.4 cm — well inside budget — and the model is visibly wrong in every place a viewer looks first.

This is not a bug in the decimator. Quadric edge collapse minimises a sum of squared distances, and the cheapest way to reduce that sum is to collapse edges in regions where many triangles describe a simple shape. A crisp break in the surface is, mathematically, an expensive feature holding a lot of triangles; removing it buys a large face reduction for a small average error. The optimiser is doing exactly what it was asked.

# Why the average is the wrong objective

Consider a flat bench with a sharp 1.2 m drop. In the full mesh the drop is described by a few thousand triangles packed along the edge. Collapsing them moves each affected vertex by a few centimetres — small individually — while removing thousands of faces. Elsewhere, on a gently curving slope, collapsing a similar number of faces moves vertices by comparable amounts but saves the same count. The optimiser sees two equally attractive trades.

To a human they are not equal at all. A few centimetres of error spread over a smooth slope is invisible. The same few centimetres applied at a break turns a discontinuity into a bevel, and the eye detects that instantly because it changes the shape class of the feature rather than its position.

The fix is to tell the decimator that some edges are not ordinary. There are three mechanisms, and using all three together costs nothing.

Figure 1 — Why a global error budget is necessary but not sufficient.

# Minimal reproducible solution

Three settings, applied together, change the outcome completely.

import pymeshlab


def decimate_preserving_breaks(in_path: str, out_path: str,
                               target_faces: int,
                               *, quality_threshold: float = 0.35,
                               boundary_weight: float = 2.0) -> dict:
    """Quadric decimation with the three break-preserving constraints enabled.

    preservenormal keeps collapses from flipping a face's orientation, which
    is what turns a crease into a fold. planarquadric adds a term that makes
    collapses inside flat regions cheaper than collapses across a change of
    plane — which is precisely the discrimination the plain algorithm lacks.
    preserveboundary protects open edges, including the outer rim of a tile.
    """
    ms = pymeshlab.MeshSet()
    ms.load_new_mesh(in_path)
    before = ms.current_mesh().face_number()

    ms.meshing_decimation_quadric_edge_collapse(
        targetfacenum=target_faces,
        qualitythr=quality_threshold,     # reject collapses making sliver faces
        preserveboundary=True,
        boundaryweight=boundary_weight,   # how expensive a boundary collapse is
        preservenormal=True,
        preservetopology=True,
        planarquadric=True,               # the key term for flat-vs-break
        autoclean=True)

    ms.save_current_mesh(out_path)
    after = ms.current_mesh().face_number()
    return {"before": before, "after": after, "reduction": 1 - after / before}

planarquadric is the setting that does the heavy lifting and the one most often left off. It adds a penalty proportional to how far the collapse departs from the local plane, which makes collapses within a flat face nearly free and collapses across a change of plane expensive. That is exactly the distinction the human eye is making.

Where the breaks are known in advance — a kerb line digitised from the orthomosaic, a bench crest from a mine plan — they can be enforced rather than inferred:

import numpy as np
import trimesh


def mark_breakline_vertices(mesh_path: str, lines_xy: list[np.ndarray],
                            tolerance: float = 0.15) -> np.ndarray:
    """Vertices within `tolerance` of a supplied breakline, to be locked.

    Locking beats weighting when the geometry is known: a vertex that must
    not move is a constraint, not a preference, and the decimator will find
    its face reduction somewhere else.
    """
    mesh = trimesh.load(mesh_path, process=False)
    v = mesh.vertices[:, :2]
    locked = np.zeros(len(mesh.vertices), dtype=bool)
    for line in lines_xy:
        for i in range(len(line) - 1):
            a, b = line[i], line[i + 1]
            ab = b - a
            t = np.clip(((v - a) @ ab) / max(ab @ ab, 1e-12), 0, 1)
            proj = a + t[:, None] * ab
            locked |= np.linalg.norm(v - proj, axis=1) < tolerance
    return locked

Figure 3 — The features worth keeping are the ones a plain error metric discards first.

# Edge-case matrix

The vegetation row is the trade-off. Break-preserving settings reduce the achievable face count on organic geometry by ten to twenty percent, because canopy has breaks everywhere. On a mixed site this is a small price; on a purely woodland model the settings can be relaxed.

# Verification snippet

Verify against a metric that reflects the failure. A ninety-fifth-percentile error catches gross damage; the useful test measures how much edge sharpness was lost.

import numpy as np
import trimesh


def sharpness_retention(original_path: str, decimated_path: str,
                        sharp_angle_deg: float = 40.0) -> dict:
    """Share of sharp dihedral edges in the original that survive decimation.

    A dihedral angle above the threshold is a break. Measuring the total
    length of such edges before and after gives a single number that tracks
    exactly the defect the average error misses.
    """
    def sharp_edge_length(path: str) -> float:
        m = trimesh.load(path, process=False)
        angles = m.face_adjacency_angles
        sharp = angles > np.radians(sharp_angle_deg)
        edges = m.face_adjacency_edges[sharp]
        if len(edges) == 0:
            return 0.0
        v = m.vertices
        return float(np.linalg.norm(v[edges[:, 0]] - v[edges[:, 1]], axis=1).sum())

    before = sharp_edge_length(original_path)
    after = sharp_edge_length(decimated_path)
    ratio = after / before if before else float("nan")
    return {"sharp_length_before_m": before, "sharp_length_after_m": after,
            "retention": ratio,
            "verdict": "acceptable" if ratio > 0.75 else "breaklines lost"}

A retention figure above about 0.75 corresponds to a model that reads correctly; below 0.5 the bevelling is obvious to anyone who looks. Tracking this number across releases is far more informative than tracking the face count, and it takes seconds.

Figure 2 — Retention, not face count, is the number worth putting in a release note.

# Choosing where the detail goes

Once retention is being measured, a better strategy than uniform decimation becomes available: spend the face budget where it matters. A site model does not need the same triangle density on a car park as on a retaining wall, and most viewers handle a mesh with varying density perfectly well.

The practical approach is a two-pass decimation driven by a region mask. Areas of interest — structures, edges of excavation, anything the client named in the brief — are decimated to a tight error budget; everything else is decimated hard. The combined mesh is smaller than a uniform one at the same perceived quality, often by half.

import numpy as np
import pymeshlab


def decimate_by_region(in_path: str, out_path: str,
                       detail_error_m: float = 0.015,
                       bulk_error_m: float = 0.08) -> dict:
    """Two-tier decimation: tight budget inside the detail mask, loose outside.

    The masks are applied as vertex quality, which pymeshlab's decimator
    reads as a per-vertex weight. A vertex with high quality is expensive to
    collapse, so the reduction happens preferentially in the bulk regions.
    """
    ms = pymeshlab.MeshSet()
    ms.load_new_mesh(in_path)
    ratio = bulk_error_m / detail_error_m
    ms.compute_scalar_by_function_per_vertex(
        q=f"(q > 0.5) ? {ratio} : 1.0")          # detail vertices weighted up
    ms.meshing_decimation_quadric_edge_collapse(
        targetfacenum=ms.current_mesh().face_number() // 8,
        qualityweight=True, preserveboundary=True,
        preservenormal=True, planarquadric=True)
    ms.save_current_mesh(out_path)
    return {"faces": ms.current_mesh().face_number()}

Two things make this worth the extra step on a regular delivery. The face budget is usually set by the client’s viewer rather than by the geometry, so it is fixed — and within a fixed budget, uniform decimation is strictly worse than weighted decimation whenever the site has areas of differing importance, which is every site. And the mask itself is cheap to produce: a buffer around the structures already classified in separating buildings from vegetation in point clouds is usually good enough, and where the client named specific features, a digitised polygon is better still.

Record the mask alongside the mesh. A model that is sharp in the wrong places is a question about the brief, and answering it a month later is far easier when the region definition shipped with the file.

# When to escalate

• Retention is poor even with the settings on. The breaks may not exist in the source mesh — a Poisson reconstruction at too low an octree depth rounds them before decimation ever runs. Check the full-resolution mesh first.
• The client needs measurable edges, not just visible ones. A mesh is the wrong deliverable. Extract the breaklines as vectors from the point cloud and deliver them alongside a coarser mesh.
• The target face count cannot be met with retention above 0.75. The site has more genuine detail than the budget allows. Either raise the budget or split the model into a coarse whole-site mesh and high-detail tiles for the areas that matter.

# Related

• Mesh and Texture Generation Automation
• Exporting OBJ and glTF for Web Viewers
• Fixing Texture Seams and Baking Artifacts

← Mesh and Texture Generation Automation