Heritage GIS

Mesh Generation & Optimization for Ruins

Within the broader Photogrammetry & 3D Site Mapping Pipelines framework, mesh generation for archaeological ruins presents distinct computational challenges. Unlike engineered or contemporary structures, ruin geometries exhibit high topological complexity, fragmented stratigraphy, and frequent occlusion-induced voids. Converting dense photogrammetric point clouds into heritage-compliant, analysis-ready meshes requires deterministic automation, explicit coordinate reference system (CRS) validation, and topology-preserving optimization routines. This cluster details a field-tested, Python-driven workflow for generating, cleaning, and validating ruin meshes prior to archival deposition or spatial analysis.

Pipeline Routing & Architectural Positioning

The mesh optimization stage operates deterministically downstream of dense cloud reconstruction. Following dense cloud extraction via Automated Drone Image Processing Workflows, the pipeline ingests .ply or .obj point clouds and executes a strict sequential routing:

  1. Ingestion & CRS Tagging → Validate input coordinate space against target EPSG.
  2. Surface Reconstruction → Screened Poisson or Ball Pivoting with stratigraphic boundary constraints.
  3. Topology Repair → Non-manifold edge collapse, duplicate vertex removal, and normal unification.
  4. Adaptive Decimation → Quadric edge collapse preserving curvature at architectural thresholds.
  5. Hole Filling & Smoothing → Curvature-preserving Laplacian smoothing for exposed stratigraphic faces.
  6. CRS-Bound Validation & Export → Fail-fast spatial sanity checks before handoff to Texture Mapping & UV Alignment Automation.

Each stage implements fail-fast validation gates. If a topology or spatial threshold is breached, execution halts, logs the artifact signature, and prevents geometric corruption from propagating into downstream GIS environments or heritage archives.

Environment & Version-Pinned Dependencies

Reproducibility in heritage automation requires strict dependency pinning. The following configuration isolates the mesh optimization environment and guarantees deterministic behavior across field deployments and institutional compute clusters.

# requirements.txt
numpy==1.26.4
pyproj==3.6.1
pymeshlab==2023.12.post2

Install via isolated environment:

python -m venv .venv_meshopt
source .venv_meshopt/bin/activate
pip install -r requirements.txt

Reference documentation: PyMeshLab Python Bindings and pyproj API Reference.

Deterministic Python Implementation

The implementation below provides a production-ready class for ruin mesh optimization. It integrates explicit EPSG validation, topology repair, and adaptive decimation with structured logging and fail-fast routing.

import os
import logging
import numpy as np
import pymeshlab as ml
from pyproj import CRS, Transformer
from pathlib import Path

# Structured logging for audit trails
logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s | %(levelname)s | %(name)s | %(message)s"
)
logger = logging.getLogger("ruin_mesh_optimizer")

class RuinMeshOptimizer:
    def __init__(self, input_path: str, output_path: str, target_epsg: int = 32633):
        self.input = Path(input_path)
        self.output = Path(output_path)
        self.target_epsg = target_epsg
        self.target_crs = CRS.from_epsg(target_epsg)
        self.ms = ml.MeshSet()
        
    def validate_crs_bounds(self, vertices: np.ndarray) -> bool:
        """Validate that mesh coordinates fall within expected projected CRS extents."""
        try:
            # Sanity threshold: projected coordinates rarely exceed 10^7 meters
            if np.any(np.abs(vertices) > 1e7):
                logger.warning(f"Coordinate magnitudes exceed typical projected CRS bounds for EPSG:{self.target_epsg}. Verify input datum.")
                return False
            
            # Optional: Transform sample bounds to WGS84 for cross-check
            transformer = Transformer.from_crs(self.target_crs, CRS.from_epsg(4326), always_xy=True)
            min_lon, min_lat = transformer.transform(vertices.min(axis=0)[0], vertices.min(axis=0)[1])
            max_lon, max_lat = transformer.transform(vertices.max(axis=0)[0], vertices.max(axis=0)[1])
            
            if not (-180 <= min_lon <= max_lon <= 180 and -90 <= min_lat <= max_lat <= 90):
                logger.error("Transformed bounds fall outside valid WGS84 geographic limits.")
                return False
            return True
        except Exception as e:
            logger.error(f"CRS validation failed: {e}")
            return False

    def run_optimization_pipeline(self) -> bool:
        """Execute deterministic mesh generation and optimization sequence."""
        if not self.input.exists():
            logger.error(f"Input file not found: {self.input}")
            return False

        try:
            logger.info(f"Loading point cloud: {self.input}")
            self.ms.load_new_mesh(str(self.input))
            m = self.ms.current_mesh()
            vertices = np.array(m.vertex_matrix())
            
            # Gate 1: Spatial Validation
            if not self.validate_crs_bounds(vertices):
                logger.error("Pipeline halted: CRS bounds validation failed.")
                return False

            # Gate 2: Topology Repair
            logger.info("Executing non-manifold repair and normal unification...")
            self.ms.apply_filter("remove_duplicate_vertices")
            self.ms.apply_filter("remove_duplicate_faces")
            self.ms.apply_filter("remove_null_faces")
            self.ms.apply_filter("re_orient_all_faces_coherentely")

            # Gate 3: Surface Reconstruction (Screened Poisson)
            logger.info("Generating screened Poisson surface...")
            self.ms.apply_filter("generate_surface_reconstruction_screened_poisson", 
                                 octdepth=10, samplespernode=1.5, fullreconstruction=True)

            # Gate 4: Adaptive Decimation (Quadric Edge Collapse)
            logger.info("Applying curvature-preserving decimation...")
            self.ms.apply_filter("meshing_decimation_quadric_edge_collapse", 
                                 targetfacenum=500000, qualitythr=0.5, preservenormal=True)

            # Gate 5: Hole Filling & Smoothing
            logger.info("Filling stratigraphic voids and smoothing...")
            self.ms.apply_filter("meshing_close_holes", maxholesize=100)
            self.ms.apply_filter("coord_laplacian_smooth", stepsmoothnum=3, cotangentweight=True)

            # Final Validation & Export
            final_mesh = self.ms.current_mesh()
            final_vertices = np.array(final_mesh.vertex_matrix())
            if final_vertices.shape[0] == 0:
                logger.error("Pipeline halted: Resulting mesh is empty.")
                return False

            self.ms.save_current_mesh(str(self.output), save_vertex_color=False, save_vertex_normal=True)
            logger.info(f"Optimization complete. Exported to: {self.output}")
            return True

        except ml.MeshSetException as e:
            logger.error(f"MeshSet operation failed: {e}")
            return False
        except Exception as e:
            logger.error(f"Unexpected pipeline failure: {e}")
            return False

if __name__ == "__main__":
    optimizer = RuinMeshOptimizer(
        input_path="data/site_A_dense_cloud.ply",
        output_path="output/site_A_optimized_mesh.obj",
        target_epsg=32633  # WGS 84 / UTM zone 33N
    )
    success = optimizer.run_optimization_pipeline()
    exit(0 if success else 1)

Validation & Archival Compliance

Heritage spatial data requires strict adherence to geodetic and topological standards before deposition. The pipeline enforces three compliance layers:

  1. EPSG/CRS Integrity: Coordinates are validated against projected bounds (e.g., EPSG:32633, EPSG:25832) and cross-checked via inverse transformation to EPSG:4326. Misaligned datums trigger immediate pipeline termination to prevent silent spatial drift in GIS environments.
  2. Topological Soundness: Non-manifold edges, null faces, and inverted normals are resolved prior to decimation. Euler characteristic checks and manifold status verification ensure compatibility with 3D GIS viewers and finite element analysis tools.
  3. Archival Metadata: Exported meshes should be accompanied by sidecar .prj or GeoPackage 3D extensions documenting the exact CRS, reconstruction parameters, and decimation targets. Reference the PROJ Coordinate Transformation Library for authoritative datum transformation chains and OGC 3D spatial standards for long-term preservation.

For academic and institutional teams, integrate this cluster into CI/CD pipelines with automated regression testing on known ruin typologies (e.g., Roman masonry, Neolithic megaliths, collapsed vaulting). Deterministic routing guarantees that mesh topology remains reproducible across software updates, field seasons, and multi-institutional research collaborations.