Mesh and Texture Generation Automation

A textured mesh is the product clients react to. It is also the only product in this section that adds no measurement accuracy whatsoever — a mesh is a resampling of the point cloud, and every step of its production discards information. Knowing that is the key to automating it well: the goal is a file that opens quickly and looks right, not one that preserves every point.

The pipeline has four stages and each has one decision that matters. Surface reconstruction turns points into triangles, and the choice of method decides whether holes are filled or left open. Decimation reduces the triangle count, and the error budget decides whether breaklines survive. UV unwrapping assigns each triangle a patch of texture, and the seam policy decides whether the result has visible joins. Export writes a format, and the format decides whether a browser can open it at all.

This page automates all four from a classified cloud produced by classifying point clouds with PDAL and Python.

Audience and prerequisites. Python 3.10+, a dense cloud with colour, and enough RAM to hold the mesh — which for a site-scale reconstruction means tiling. Meshing is the most memory-hungry stage in the whole photogrammetry chain.

# Prerequisites

# Conceptual architecture

Surface reconstruction has two families and they fail in opposite directions. Poisson reconstruction fits an implicit function to oriented points and extracts an isosurface. It produces watertight, smooth meshes and always closes holes — including the ones that represent genuinely missing data, which it fills with confident, invented geometry. Ball pivoting and the Delaunay family interpolate between actual points and leave holes where there are no points, producing meshes that are honest and often ragged.

For a visual deliverable, Poisson is usually right: the client wants a solid-looking model and an invented patch of ground under a tree is not being measured. For anything a measurement might be taken from, ball pivoting is right, because its holes are the truth. Making that choice explicit in the job definition — and recording it — is what stops a viewer’s measurement tool being used on invented geometry.

Figure 1 — The reconstruction choice is a policy about honesty, not a quality setting.

# Step 1: Estimate normals before reconstructing anything

Poisson reconstruction needs oriented normals, and the orientation is the part that goes wrong. A normal estimated from a local neighbourhood is ambiguous in sign; if the signs are inconsistent the surface folds back on itself and the result is a mesh full of spikes.

import numpy as np
import open3d as o3d


def prepare_cloud(points: np.ndarray, colors: np.ndarray,
                  *, radius: float = 0.5, max_nn: int = 30,
                  viewpoint_height: float = 120.0) -> o3d.geometry.PointCloud:
    """Load a cloud and orient its normals consistently upward.

    Orienting toward a synthetic camera high above the site is the right
    heuristic for aerial photogrammetry: every surface the drone saw faces
    broadly upward, so a single viewpoint resolves the sign ambiguity without
    the expensive minimum-spanning-tree propagation.
    """
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(points)
    pcd.colors = o3d.utility.Vector3dVector(colors / 255.0)
    pcd.estimate_normals(
        search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=radius,
                                                          max_nn=max_nn))
    centre = points.mean(axis=0)
    pcd.orient_normals_towards_camera_location(
        np.array([centre[0], centre[1], centre[2] + viewpoint_height]))
    return pcd

The viewpoint trick matters for building facades, which face sideways and are the one place the heuristic is weak. Where facades are part of the deliverable, orient with the minimum-spanning-tree method instead and accept the runtime.

# Step 2: Reconstruct, then trim by density

Poisson’s most useful output is not the mesh but the per-vertex density it reports alongside it. Low-density vertices are exactly the ones the algorithm invented, and trimming them turns Poisson from a hole-filler into something closer to honest.

import numpy as np
import open3d as o3d


def poisson_mesh(pcd, *, depth: int = 11, density_percentile: float = 4.0):
    """Poisson reconstruction with low-confidence vertices trimmed away.

    depth controls the octree resolution: each increment doubles the linear
    resolution and roughly quadruples memory. 11 suits a site at a few
    centimetres; 13 will exhaust most machines on anything site-scale.
    """
    mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
        pcd, depth=depth, width=0, scale=1.1, linear_fit=False)

    d = np.asarray(densities)
    cutoff = np.quantile(d, density_percentile / 100.0)
    mesh.remove_vertices_by_mask(d < cutoff)
    mesh.remove_degenerate_triangles()
    mesh.remove_duplicated_vertices()
    mesh.remove_non_manifold_edges()
    return mesh

Trimming at the fourth percentile removes the confident invention over large gaps while leaving the small interpolations that make the mesh pleasant. Raising it toward ten percent produces a mesh close to ball-pivoting honesty with Poisson’s smoothness.

# Step 3: Decimate against an error budget, not a face count

“Reduce to two million faces” is the instruction most pipelines carry, and it is the wrong shape of instruction: it produces different geometric error on every site. The right one is “reduce until the mesh deviates from the original by more than 3 cm”, which is a statement in survey units that a reader can evaluate.

import numpy as np
import pymeshlab


def decimate_to_error(in_path: str, out_path: str, *,
                      max_error_m: float = 0.03,
                      floor_faces: int = 50_000) -> dict:
    """Binary-search the face count that just meets a geometric error budget."""
    ms = pymeshlab.MeshSet()
    ms.load_new_mesh(in_path)
    original_faces = ms.current_mesh().face_number()

    lo, hi = floor_faces, original_faces
    best = None
    while lo < hi:
        target = (lo + hi) // 2
        ms.load_new_mesh(in_path)
        ms.meshing_decimation_quadric_edge_collapse(
            targetfacenum=target, preserveboundary=True, preservenormal=True,
            preservetopology=True, planarquadric=True)
        ms.apply_filter("compute_scalar_by_distance_from_another_mesh_per_vertex",
                        measuremesh=ms.current_mesh_id(), refmesh=0)
        err = float(np.percentile(np.abs(ms.current_mesh().vertex_scalar_array()), 95))
        if err <= max_error_m:
            best = (target, err)
            hi = target
        else:
            lo = target + 1

    if best is None:
        raise ValueError(f"cannot meet {max_error_m} m even at {floor_faces} faces")
    ms.save_current_mesh(out_path)
    return {"faces": best[0], "p95_error_m": best[1],
            "reduction": 1 - best[0] / original_faces}

preserveboundary and planarquadric together are what keep breaklines — the top of a kerb, the edge of a bench — from being smoothed away. Decimation without them produces a mesh that is dimensionally correct on average and visibly wrong exactly where a viewer looks. Decimating meshes without losing breaklines covers the failure in detail.

Figure 2 — The budget is the instruction; the face count is the answer. Reversing the two is why meshes ship at arbitrary sizes.

# Step 4: Export a format the client can actually open

The format decision is usually made badly because OBJ is familiar. OBJ is a text format: a five-million-triangle mesh becomes a 600 MB file that takes a browser thirty seconds to parse. glTF binary (.glb) stores the same mesh in about 90 MB and loads in a second, with the texture embedded rather than shipped as a sidecar that gets lost.

from pathlib import Path

import numpy as np
import trimesh


def export_glb(mesh_path: str, texture_path: str, out_path: str,
               *, origin_shift: tuple[float, float, float] | None = None) -> dict:
    """Export a textured mesh as a single self-contained .glb.

    The origin shift is not cosmetic: glTF stores positions as 32-bit floats,
    and a UTM easting of 500 000 leaves about 3 cm of resolution. Shifting to
    a local origin and recording the offset preserves millimetres.
    """
    mesh = trimesh.load(mesh_path, process=False)
    if origin_shift is None:
        origin_shift = tuple(np.asarray(mesh.bounds).mean(axis=0))
    mesh.apply_translation(-np.asarray(origin_shift))

    image = trimesh.load_image = None  # placeholder: material assembled below
    material = trimesh.visual.material.SimpleMaterial(
        image=trimesh.visual.texture.Image.open(texture_path))
    mesh.visual = trimesh.visual.TextureVisuals(uv=mesh.visual.uv,
                                                material=material)
    mesh.export(out_path)
    return {"origin_shift": list(map(float, origin_shift)),
            "faces": int(len(mesh.faces)),
            "bytes": Path(out_path).stat().st_size}

The float32 problem in the comment is the one that produces “my model looks fine but jitters when I zoom in”. glTF has no double-precision position type, so a mesh exported in projected coordinates is quantised to a few centimetres before it ever reaches the viewer. Shifting to a local origin and recording the offset in the metadata is the standard fix, and exporting OBJ and glTF for web viewers covers the rest of the format’s constraints.

# Step 5: Unwrap and bake the texture

Between decimation and export sits the stage that decides how the model actually looks. Texturing has two halves: assigning every triangle a region of a two-dimensional image (unwrapping), and filling that image with colour sampled from the source photographs (baking).

Unwrapping is a packing problem. The mesh is cut into patches, each patch is flattened, and the flattened patches are arranged in an atlas. Three parameters govern the outcome. Patch count trades distortion against seams: fewer, larger patches mean fewer visible joins but more stretching. Atlas resolution sets the texel density, and 8192 pixels square is the practical ceiling — many mobile viewers refuse larger textures outright. Gutter width is the padding between packed patches, and it is the parameter that prevents the most common visible artefact.

import numpy as np
import xatlas


def unwrap(vertices: np.ndarray, faces: np.ndarray,
           *, resolution: int = 8192, gutter_px: int = 4) -> dict:
    """Cut, flatten and pack the mesh into a UV atlas.

    The gutter matters more than it sounds. Mip-mapping averages neighbouring
    texels, so patches packed edge to edge bleed into each other at every
    zoom level below full resolution — which is most of the time a viewer
    spends. Four pixels of padding costs about one percent of the atlas and
    removes the artefact entirely.
    """
    atlas = xatlas.Atlas()
    atlas.add_mesh(vertices, faces)
    opts = xatlas.PackOptions()
    opts.resolution = resolution
    opts.padding = gutter_px
    opts.bruteForce = False
    atlas.generate(pack_options=opts)

    vmap, indices, uvs = atlas[0]
    return {"vertices": vertices[vmap], "faces": indices, "uv": uvs,
            "atlas_px": resolution, "charts": atlas.chart_count,
            "utilisation": float(atlas.utilization[0])}

Atlas utilisation below about sixty percent is a signal worth acting on: it means the packer could not fit the patches efficiently, usually because the mesh was cut into many small charts by a noisy surface. Smoothing lightly before unwrapping, or reducing the chart count, recovers both texture quality and file size.

Baking then samples each texel from the source imagery. The decision that matters here is which photograph to sample from where several saw the same triangle, and the naive answer — the first one — produces a patchwork of exposures. Choosing by viewing angle and blending across the boundary produces a texture that looks continuous.

import numpy as np


def choose_source_image(triangle_normal: np.ndarray,
                        camera_dirs: np.ndarray,
                        camera_distances: np.ndarray,
                        max_angle_deg: float = 60.0) -> int | None:
    """Pick the best photograph to texture one triangle from.

    Score by how squarely the camera faced the surface, penalised by
    distance. Rejecting oblique views outright avoids the smeared texture
    that comes from sampling a surface seen almost edge-on.
    """
    cos = camera_dirs @ (triangle_normal / np.linalg.norm(triangle_normal))
    ok = cos > np.cos(np.radians(max_angle_deg))
    if not ok.any():
        return None                     # no acceptable view: leave a hole to fill
    score = np.where(ok, cos / np.maximum(camera_distances, 1e-6), -np.inf)
    return int(np.argmax(score))

Returning None rather than falling back to a bad view is deliberate. A triangle with no square-on observation is better left to a hole-filling pass, which will interpolate from its neighbours, than textured from a photograph that saw it at eighty degrees and recorded four pixels of smear.

Two operational notes close the stage. Bake at the atlas resolution you will ship, not higher — downsampling a 16k atlas to 8k afterwards blurs the seam gutters back into each other. And record the atlas resolution, chart count and utilisation in the run metadata, because a mesh that looks worse than last month’s is almost always one whose packing degraded rather than one whose geometry did.

# Parameter deep-dive

# Verification and output inspection

Mesh quality is checkable without opening a viewer, and every check below has caught a real defect.

import numpy as np
import trimesh


def inspect_mesh(path: str, *, expect_watertight: bool) -> dict:
    """Structural checks a delivered mesh must pass."""
    m = trimesh.load(path, process=False)
    report = {
        "faces": int(len(m.faces)),
        "vertices": int(len(m.vertices)),
        "watertight": bool(m.is_watertight),
        "winding_consistent": bool(m.is_winding_consistent),
        "has_uv": m.visual.uv is not None if hasattr(m.visual, "uv") else False,
        "degenerate_faces": int(np.count_nonzero(m.area_faces < 1e-12)),
        "bounds": m.bounds.tolist(),
    }
    problems = []
    if expect_watertight and not report["watertight"]:
        problems.append("expected a watertight mesh and got an open one")
    if not report["winding_consistent"]:
        problems.append("inconsistent winding — the viewer will show black patches")
    if not report["has_uv"]:
        problems.append("no UV coordinates — the texture cannot be applied")
    if report["degenerate_faces"]:
        problems.append(f"{report['degenerate_faces']} zero-area faces")
    report["problems"] = problems
    return report

Inconsistent winding is the defect worth checking hardest, because its symptom — patches that render black from one angle and correctly from another — is easily mistaken for a texture problem and sends people down the wrong path for an afternoon.

Figure 3 — Five stages, each of which can produce a defect the next one cannot fix.

# Troubleshooting

The mesh has spikes all over it. Normal orientation is inconsistent. Re-orient toward a viewpoint above the site, and check that the cloud was denoised first — a single outlier drags a spike through the surface.

Poisson produced a bubble around the whole site. The octree extended beyond the data and closed the surface at its edge. Trim by density, and crop the result to the cloud’s own extent.

The mesh looks right and the texture is grey. Either the UV coordinates were lost during decimation, or the material references a texture file that was not exported alongside it. Export .glb, which embeds it.

The model jitters or shimmers when zoomed in. Positions are in projected coordinates and have been quantised to float32. Shift to a local origin before export.

Meshing runs out of memory. Poisson at depth 12 on a site-scale cloud needs tens of gigabytes. Tile the cloud with a buffer, mesh per tile, and merge — the same discipline as tiled classification.

Black patches that move as the camera moves. Inconsistent face winding. Run a winding repair before export rather than adjusting the viewer.

# Related

• Decimating Meshes Without Losing Breaklines
• Fixing Texture Seams and Baking Artifacts
• Exporting OBJ and glTF for Web Viewers
• Classifying Point Clouds with PDAL and Python
• Memory Management for Large Point Clouds

← Point Cloud Processing & 3D Deliverables