Fixing Out-of-Memory Failures in City-Scale Mesh Decimation

This guide prevents the MemoryError and OOM-killer terminations that strike when automated mesh decimation and point-cloud processing scale from a few buildings to a whole city, using tile-by-tile streaming instead of whole-survey loads, chunked laspy reads, a bounded ProcessPoolExecutor, explicit freeing of intermediates, and resident-set measurement with psutil. It sits in the cross-section failure modes area because the failure spans the boundary between the processing pipeline and the LOD tiler that consumes its output.

You hit this the first time a pipeline that ran fine on a test block is pointed at a full survey: a job that decimated ten buildings in seconds tries to load a 40 GB .laz or a directory of thousands of meshes at once, resident memory climbs past the machine’s RAM, and the process is killed with no traceback or a bare MemoryError. The fix is architectural — never hold the whole survey — and it is the same discipline whether the payload is points or triangles.

Whole-survey load versus streamed tile-by-tile processing The top path loads the entire survey into memory at once and crosses the RAM budget, triggering the out-of-memory killer. The bottom path streams one chunk or tile at a time through a bounded pool of workers, holding resident memory flat below the budget line regardless of survey size. Whole survey40 GB .laz load all pointsRSS > RAM budget OOM kill (137) chunk_iteratorone chunk / tile bounded pooldel + gc.collect flat RSSunder budget Peak memory tracks one chunk times worker count — not the survey size
Loading the whole survey crosses the RAM budget and invites the OOM killer; streaming one chunk or tile through a memory-bounded pool holds resident set flat no matter how large the survey grows.

Prerequisites

  • Python 3.10+ with laspy 2.5+ (with the lazrs or laszip backend for .laz), numpy 1.26+, trimesh 4.4+, and psutil 5.9+: pip install "laspy[lazrs]>=2.5" "numpy>=1.26" "trimesh>=4.4" "psutil>=5.9".
  • Source data as classified .laz/.las point clouds or per-tile .glb/.ply meshes in a projected metric CRS — the examples use EPSG:32618 (UTM 18N) so that chunk extents and point counts are in metres and are directly comparable across tiles.
  • A machine whose RAM you know: the whole point of the exercise is to keep peak resident-set size (RSS) below that ceiling with headroom for the OS, so measure it (psutil.virtual_memory().total) rather than assuming.

Step-by-Step

1. Measure peak RSS so “out of memory” becomes a number

Before changing anything, instrument the job. Without a measured RSS you are guessing, and the fix cannot be verified. psutil reads the process’s resident set; sample it around the heavy call.

python
import os, psutil

proc = psutil.Process(os.getpid())

def rss_mb():
    return proc.memory_info().rss / 1024**2

budget_mb = psutil.virtual_memory().total / 1024**2 * 0.6   # leave 40% headroom
print(f"RSS at start {rss_mb():.0f} MB, soft budget {budget_mb:.0f} MB")

Keep the soft budget well below total RAM — decimation allocates temporary index and normal arrays several times the size of the geometry, so 60% of physical memory is a safe working ceiling.

2. Stream a large point cloud in fixed-size chunks

The core fix for point data: never call laspy.read, which materialises every point. Open the file and iterate fixed-size chunks so resident memory is bounded by the chunk, not the survey.

python
import laspy
import numpy as np

CHUNK = 5_000_000   # points per chunk; ~120 MB for XYZ float64 + attributes

def stream_ground_points(path):
    """Yield decimated ground points chunk by chunk, holding one chunk at a time."""
    with laspy.open(path) as reader:          # header only — no points loaded yet
        assert reader.header.parse_crs().to_epsg() == 32618, "expected EPSG:32618"
        for chunk in reader.chunk_iterator(CHUNK):
            xyz = np.vstack((chunk.x, chunk.y, chunk.z)).T      # this chunk only
            ground = xyz[chunk.classification == 2]             # ASPRS ground class
            keep = ground[::4]                                  # 4x voxel-free thin
            yield keep
            del xyz, ground                                    # free before next chunk

Each yield hands downstream a small array; the del drops the chunk’s working buffers so the next iteration starts from the header cursor, not a growing accumulation.

Resident memory over the life of a job Loading the whole survey produces a resident set that climbs monotonically until it crosses the machine limit and the process is killed. Streaming one chunk at a time produces a sawtooth that peaks at one chunk plus the worker count, and never approaches the limit however large the survey grows. RAM limit killed — exit 137 one chunk × workers, flat forever elapsed time RSS The two shapes are diagnostic on their own — a climbing line is never a tuning problem Log RSS every few seconds and the distinction costs nothing; without it, both failures read as "the job died"
A monotonic climb means something is being accumulated that should not be. A sawtooth that is simply too tall means the chunk is too large. The remedies are unrelated.

3. Decimate one tile at a time, freeing intermediates

For meshes, the same rule applies per file. Load one tile, decimate, export, and explicitly drop every intermediate before the next tile — do not build a list of loaded meshes.

python
import gc, trimesh
from pathlib import Path

def decimate_tile(glb_path: Path, out_dir: Path, ratio=0.125):
    mesh = trimesh.load(glb_path, force="mesh", process=False)
    target = max(500, int(len(mesh.faces) * ratio))
    simplified = mesh.simplify_quadric_decimation(target)
    out = out_dir / (glb_path.stem + "_lod1.glb")
    simplified.export(out)
    faces_in, faces_out = len(mesh.faces), len(simplified.faces)
    del mesh, simplified          # drop both meshes' vertex/face/index buffers
    gc.collect()                  # reclaim before the next tile is loaded
    return out, faces_in, faces_out

Calling gc.collect() after del forces trimesh’s numpy-backed buffers to be reclaimed immediately rather than at the next allocation pressure, which is what keeps RSS flat across thousands of tiles.

Where the copies come from in one decimation step A single decimation of a two gigabyte mesh holds the source arrays, the library's internal copy, the output arrays and, until the reference is released, the previous iteration's result. Four live copies of nearly the same data is the normal peak, not a leak. source arrays2.0 GB library copy2.0 GB output arrays0.6 GB previous resultstill referenced0.6 GB Peak is 5.2 GB for a 2.0 GB mesh — and the only one you can delete is the fourth del the previous result before the next call, not after it, or both live at once for the duration Budget from the peak, never from the file size
Sizing a job from the input file is the most common way to under-provision it by a factor of two or three. The peak is what has to fit.

Both figures point at the same discipline: measure rather than estimate. Peak resident set is one call to resource.getrusage at the end of a run, or a background thread sampling psutil.Process().memory_info().rss every few seconds if you want the shape as well as the maximum. Logging that number alongside the chunk size and worker count turns capacity planning into interpolation, and turns “it OOMs on the big city” into a ratio you can act on.

4. Parallelise with a bounded process pool

Decimation is CPU-bound and embarrassingly parallel per tile, but an unbounded pool multiplies peak memory by the worker count and re-triggers the OOM you just fixed. Size the pool to fit within the budget: each worker holds roughly one tile’s footprint, so workers ≤ budget / per_tile_peak.

python
from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path
import os

PER_TILE_PEAK_MB = 700          # measured peak RSS of one decimate_tile call
max_workers = max(1, min(os.cpu_count(), int(budget_mb // PER_TILE_PEAK_MB)))

def run(tiles, out_dir):
    out_dir = Path(out_dir); out_dir.mkdir(exist_ok=True)
    results = []
    with ProcessPoolExecutor(max_workers=max_workers) as pool:
        futures = {pool.submit(decimate_tile, t, out_dir): t for t in tiles}
        for fut in as_completed(futures):        # collect as they finish, not all at once
            results.append(fut.result())
    return results

print(f"using {max_workers} workers within {budget_mb:.0f} MB budget")

Submitting all futures up front is fine — they queue — but bounding max_workers by the memory budget, not just the core count, is what prevents the pool from over-committing RAM.

5. Aggregate streamed results without reloading

When a downstream step needs all the thinned points (for example to build one tileset extent), accumulate running statistics or write to a memory-mapped array, never a growing Python list of full chunks.

python
import numpy as np

def survey_bounds(paths):
    """Fold streamed chunks into a bounds accumulator — O(1) memory in survey size."""
    lo = np.full(3, np.inf)
    hi = np.full(3, -np.inf)
    total = 0
    for path in paths:
        for keep in stream_ground_points(path):
            lo = np.minimum(lo, keep.min(axis=0))
            hi = np.maximum(hi, keep.max(axis=0))
            total += len(keep)
    return lo, hi, total

lo, hi, total = survey_bounds([Path("survey") / f for f in ("n.laz", "s.laz")])
print(f"kept {total:,} points, extent {(hi - lo)[:2]} m (EPSG:32618)")

Expected Output & Verification

A correctly streamed run holds RSS roughly flat regardless of survey size — resident memory tracks one chunk or one tile plus the worker multiple, not the total data volume.

text
RSS at start 92 MB, soft budget 19660 MB
using 12 workers within 19660 MB budget
tile block_0003.glb: 214880 -> 26860 faces
tile block_0004.glb: 198122 -> 24765 faces
kept 41,203,588 points, extent [1188.4 1421.7] m (EPSG:32618)
peak RSS 8140 MB (budget 19660 MB) — OK

Verify the fix by asserting peak RSS stayed under budget across the whole run, and by confirming peak is independent of survey size — process a survey twice as large and peak RSS should barely move.

python
peak = rss_mb()
assert peak < budget_mb, f"peak RSS {peak:.0f} MB exceeded budget {budget_mb:.0f} MB"
print(f"peak RSS {peak:.0f} MB (budget {budget_mb:.0f} MB) — OK")

If peak scales with survey size, an intermediate is still being accumulated — a list of chunks, an un-del-eted mesh, or an unbounded pool — and the offending step is the one whose output grows.

Log the peak resident set alongside the chunk size and worker count on every run. Capacity planning then becomes interpolation from your own numbers rather than a guess, and a regression in memory behaviour shows up as a trend instead of as an outage.

Common Errors

MemoryError or a silent kill with exit code 137. Exit 137 is the Linux OOM killer (128 + SIGKILL), which leaves no Python traceback, so the job just vanishes from the scheduler. It means resident memory crossed the physical ceiling, almost always because laspy.read or trimesh.load was called on the whole survey. Switch to laspy.open(...).chunk_iterator and per-tile loading as in steps 2 and 3, and cap the pool by memory in step 4.

laspy raises LaspyException: laszip is not available on .laz. The compressed reader backend is missing, so chunk_iterator cannot decode the file. Install a backend with pip install "laspy[lazrs]" (pure-Rust, no system dependency) or laspy[laszip], and confirm with laspy.LazBackend.detect_available() before the run. Reading .las uncompressed masks this but multiplies disk and memory cost.

Peak RSS climbs slowly across thousands of tiles despite per-tile loading. numpy-backed buffers from trimesh are not always reclaimed at the moment the Python reference drops, so memory creeps until an allocation forces collection. Add an explicit del of every mesh and intermediate followed by gc.collect() at the end of each tile (step 3); without the forced collection, the drift eventually crosses the budget on a long enough run.

Back to Cross-Section Failure Modes in Digital Twins.