Mesh Operations Module

Low-level mesh generation and boundary operations.

This module contains functions for creating and manipulating terrain meshes, extracted from the core Terrain class for better modularity and testability.

Optimized with Numba JIT compilation for performance-critical loops.

Core Mesh Generation

src.terrain.mesh_operations.generate_vertex_positions(dem_data, valid_mask, scale_factor=100.0, height_scale=1.0)

Generate 3D vertex positions from DEM data.

Converts 2D elevation grid into 3D positions for mesh vertices, applying scaling factors for visualization. Only generates vertices for valid (non-NaN) DEM values.

Parameters:

• dem_data (np.ndarray) – 2D array of elevation values (height x width)

• valid_mask (np.ndarray) – Boolean mask indicating valid points (True for non-NaN)

• scale_factor (float) – Horizontal scale divisor for x/y coordinates (default: 100.0). Higher values produce smaller meshes. E.g., 100 means 100 DEM units = 1 unit.

• height_scale (float) – Multiplier for elevation values (default: 1.0). Values > 1 exaggerate terrain, < 1 flatten it.

Returns:

(positions, y_valid, x_valid) where:

• positions: np.ndarray of shape (n_valid, 3) with (x, y, z) coordinates

• y_valid: np.ndarray of y indices for valid points

• x_valid: np.ndarray of x indices for valid points

Return type:

tuple

Generate 3D vertex positions from DEM data.

Converts 2D elevation grid to 3D mesh vertices with scaling.

Example:

from src.terrain.mesh_operations import generate_vertex_positions

vertices = generate_vertex_positions(
    dem_data,
    valid_mask=~np.isnan(dem_data),
    scale_factor=100.0,
    height_scale=30.0
)

src.terrain.mesh_operations.generate_faces(height, width, coord_to_index, batch_size=10000)

Generate mesh faces from a grid of valid points.

Creates quad faces for the mesh by checking each potential quad position and verifying that its corners exist in the coordinate-to-index mapping. If a quad has all 4 corners, creates a quad face. If it has 3 corners, creates a triangle face. Skips quads with fewer than 3 corners.

Parameters:

• height (int) – Height of the DEM grid

• width (int) – Width of the DEM grid

• coord_to_index (dict) – Mapping from (y, x) coordinates to vertex indices

• batch_size (int) – Number of quads to process in each batch (default: 10000)

Returns:

List of face tuples, where each tuple contains vertex indices

Return type:

list

Generate triangle faces for mesh with batch processing.

Uses Numba JIT compilation for ~10x speedup on large meshes.

Example:

from src.terrain.mesh_operations import generate_faces

faces = generate_faces(
    height=dem.shape[0],
    width=dem.shape[1],
    coord_to_index=coord_map,
    batch_size=10000
)

src.terrain.mesh_operations._generate_faces_numba(height, width, index_grid)

Numba-accelerated face generation.

Parameters:

• height – Grid height

• width – Grid width

• index_grid – 2D array where index_grid[y,x] = vertex index, or -1 if invalid

Returns:

Array of face vertex indices (n_faces, 4), padded with -1 for triangles n_faces: Number of valid faces

Return type:

faces

Numba-accelerated face generation (internal).

Performance: ~10x faster than pure Python for 1M+ faces.

Boundary Detection

src.terrain.mesh_operations.find_boundary_points(valid_mask)

Find boundary points using morphological operations.

Identifies points on the edge of valid regions using binary erosion. A point is considered a boundary point if it is valid but has at least one invalid neighbor in a 4-connected neighborhood.

Parameters:

valid_mask (np.ndarray) – Boolean mask indicating valid points (True for valid)

Returns:

List of (y, x) coordinate tuples representing boundary points

Return type:

list

Find boundary points using morphological operations.

Identifies edge pixels where valid DEM data meets invalid/NaN regions.

Example:

from src.terrain.mesh_operations import find_boundary_points

valid_mask = ~np.isnan(dem_data)
boundary = find_boundary_points(valid_mask)
print(f"Found {len(boundary)} boundary points")

Boundary Extension

src.terrain.mesh_operations.create_boundary_extension(positions, boundary_points, coord_to_index, base_depth=0.2, two_tier=False, mid_depth=None, base_material='clay', blend_edge_colors=True, surface_colors=None, smooth_boundary=False, smooth_window_size=5, use_catmull_rom=False, catmull_rom_subdivisions=2, use_rectangle_edges=False, dem_shape=None, terrain=None, edge_sample_spacing=0.33, boundary_winding='counter-clockwise', use_fractional_edges=False, scale_factor=100.0, model_offset=None)

Create boundary extension vertices and faces to close the mesh.

Creates a “skirt” around the terrain by adding bottom vertices at base_depth and connecting them to the top boundary with quad faces. This closes the mesh into a solid object suitable for 3D printing or solid rendering.

Supports two modes: - Single-tier (default): Surface → Base (one jump) - Two-tier: Surface → Mid → Base (two-tier with color separation)

Parameters:

• positions (np.ndarray) – Array of (n, 3) vertex positions

• boundary_points (list) – List of (y, x) tuples representing ordered boundary points

• coord_to_index (dict) – Mapping from (y, x) coordinates to vertex indices

• base_depth (float) – Positive depth offset below minimum surface elevation (default: 0.2). Creates a flat base plane at: min_surface_z - base_depth. Positive values extend below surface, negative extend above.

• two_tier (bool) – Enable two-tier mode (default: False)

• mid_depth (float, optional) – Positive depth offset below surface for mid tier (default: base_depth * 0.25, typically 0.05). Positive values extend below surface, negative extend above.

• base_material (str | tuple) – Material for base layer - either preset name (“clay”, “obsidian”, “chrome”, “plastic”, “gold”, “ivory”) or RGB tuple (0-1 range). Default: “clay”

• blend_edge_colors (bool) – Blend surface colors to mid tier (default: True) If False, mid tier uses base_material color for sharp transition

• surface_colors (np.ndarray, optional) – Surface vertex colors (n_vertices, 3) uint8

• smooth_boundary (bool) – Apply smoothing to boundary to eliminate stair-step edges (default: False)

• smooth_window_size (int) – Window size for boundary smoothing (default: 5). Larger values produce smoother curves.

• use_catmull_rom (bool) – Use Catmull-Rom curve fitting for smooth boundary instead of pixel-grid topology (default: False). When enabled, eliminates staircase pattern entirely. NOTE: Computationally expensive (~0.3-2s per terrain). Provides true smooth curves vs simple smoothing.

• catmull_rom_subdivisions (int) – Number of interpolated points per boundary segment when using Catmull-Rom curves (default: 2). Higher values = smoother curve but MORE COMPUTATION. Recommended: 2 (fast) or 3-4 (very smooth).

• use_rectangle_edges (bool) – Use rectangle-edge sampling instead of morphological boundary detection (default: False). ~150x faster than morphological detection. Ideal for rectangular DEMs from raster sources.

• dem_shape (tuple, optional) – DEPRECATED - DEM shape (height, width) for legacy rectangle-edge sampling. Use terrain= parameter instead for transform-aware edges (avoids NaN margins).

• terrain (Terrain, optional) – Terrain object for transform-aware rectangle-edge sampling. Provides original DEM shape and transform pipeline for accurate coordinate mapping without NaN margins. Improves edge coverage from 0.6% (legacy) to ~100% (transform-aware) for downsampled DEMs.

• edge_sample_spacing (float) – Pixel spacing for edge sampling at original DEM resolution (default: 1.0). Lower values = denser sampling, more edge pixels.

• use_fractional_edges (bool) – Use fractional coordinates that preserve projection curvature (default: False). When True, creates smooth curved edge by: 1. Surface tier aligned with mesh boundary (bilinear interpolation, no gap) 2. Mid tier at fractional X,Y positions with offset Z (smooth curve below surface) 3. Base tier at fractional X,Y positions with flat Z (smooth curved base) This eliminates gaps while preserving smooth projection-aware edge curves. Requires terrain= parameter.

Returns:

When two_tier=False (backwards compatible):

(boundary_vertices, boundary_faces)

tuple: When two_tier=True:

(boundary_vertices, boundary_faces, boundary_colors)

Where:

• boundary_vertices: np.ndarray of vertex positions

  Single-tier: (n_boundary, 3) Two-tier: (2*n_boundary, 3) - mid + base vertices

• boundary_faces: list of tuples defining side face quad connectivity

  Single-tier: N quads (surface→base) Two-tier: 2*N quads (surface→mid + mid→base)

• boundary_colors: np.ndarray of (2*n_boundary, 3) uint8 colors (two-tier only)

Return type:

tuple

Create two-tier edge extrusion for terrain mesh.

Main function for edge generation - supports multiple modes:

• Single-tier vs two-tier

• Morphological vs rectangle edge detection

• Standard vs fractional edge positioning

• Catmull-Rom curve smoothing

Example:

from src.terrain.mesh_operations import create_boundary_extension

# Basic two-tier edge
edge_data = create_boundary_extension(
    processed_dem,
    terrain.y_valid,
    terrain.x_valid,
    valid_mask,
    scale_factor=100.0,
    height_scale=30.0,
    two_tier=True,
    edge_color=(0.5, 0.48, 0.45),  # Clay
    base_color=(0.3, 0.3, 0.3)
)

# With Catmull-Rom smoothing
edge_data = create_boundary_extension(
    ...,
    use_catmull_rom=True,
    catmull_rom_subdivisions=20
)

# With fractional edges (smooth projection curves)
edge_data = create_boundary_extension(
    ...,
    use_fractional_edges=True,
    use_rectangle_edges=True
)

See Combined Render: Full-Featured Example for usage examples.

Catmull-Rom Curve Fitting

src.terrain.mesh_operations.catmull_rom_curve(p0, p1, p2, p3, t)

Catmull-Rom spline interpolation between two points.

Evaluates a Catmull-Rom spline at parameter t, using four control points. The curve passes through p1 and p2, and is influenced by p0 and p3.

Parameters:

• p0 – Control points (numpy arrays or tuples)

• p1 – Control points (numpy arrays or tuples)

• p2 – Control points (numpy arrays or tuples)

• p3 – Control points (numpy arrays or tuples)

• t – Parameter in [0, 1] where 0=p1 and 1=p2

Returns:

Point on the curve at parameter t

Evaluate Catmull-Rom spline at parameter t.

Catmull-Rom properties:

• Passes through all control points

• C¹ continuous (smooth first derivative)

• Local control (changing p1/p2 only affects nearby curve)

src.terrain.mesh_operations.fit_catmull_rom_boundary_curve(boundary_points, subdivisions=10, closed_loop=True)

Fit a Catmull-Rom spline curve through boundary points.

Creates a smooth curve that passes through all boundary points by fitting Catmull-Rom spline segments between consecutive points.

Parameters:

• boundary_points – List of (y, x) or (x, y) tuples representing the boundary

• subdivisions – Number of interpolated points per segment (higher = smoother)

• closed_loop – If True, treat the boundary as a closed loop

Returns:

List of interpolated points along the smooth curve

Fit smooth Catmull-Rom curve through boundary points.

Eliminates pixel-grid staircase pattern from morphological boundaries.

Example:

from src.terrain.mesh_operations import fit_catmull_rom_boundary_curve

# Smooth boundary with 20 points per segment
smooth_boundary = fit_catmull_rom_boundary_curve(
    boundary_points,
    subdivisions=20,
    closed_loop=True
)

Boundary Processing

src.terrain.mesh_operations.smooth_boundary_points(boundary_coords, window_size=3, closed_loop=True)

Smooth boundary points using moving average to eliminate stair-step edges.

Applies a moving average filter to boundary coordinates to create smoother curves instead of following pixel grid exactly. This reduces the jagged appearance on curved edges while preserving overall shape.

Parameters:

• boundary_coords – List of (y, x) coordinate tuples representing boundary points

• window_size – Size of smoothing window (must be odd, default: 3). Larger values produce more smoothing.

• closed_loop – If True, treat boundary as closed loop (wrap edges). If False, treat as open path (endpoints less smoothed).

Returns:

Smoothed boundary points as list of (y, x) float tuples

Return type:

list

Examples

>>> boundary = [(0, 0), (0, 1), (1, 1), (1, 2)]
>>> smoothed = smooth_boundary_points(boundary, window_size=3)
>>> # Returns smoothed coordinates with reduced stair-stepping

Apply moving average smoothing to boundary.

src.terrain.mesh_operations.deduplicate_boundary_points(boundary_coords)

Remove duplicate points while preserving the original order.

After coordinate transformations, many boundary points map to the same pixel coordinates, creating duplicates. This function removes duplicates while preserving the original perimeter traversal order.

Parameters:

boundary_coords – List of (y, x) coordinate tuples

Returns:

Deduplicated boundary points in original order

Return type:

list

Remove duplicate consecutive points.

src.terrain.mesh_operations.sort_boundary_points_angular(boundary_coords)

Sort boundary points by angle from centroid to form a closed loop.

This is much faster than nearest-neighbor sorting and works well for dense boundaries (>10K points). Computes the centroid of all boundary points, then sorts by angle, creating a natural closed loop around the perimeter.

After angular sorting, rotates the list so the largest gap between consecutive points becomes the start/end, preventing diagonal faces across the mesh.

Parameters:

boundary_coords – List of (y, x) coordinate tuples representing boundary points

Returns:

Sorted boundary points forming a continuous closed loop

Return type:

list

Sort boundary points by angular position (for closed loops).

src.terrain.mesh_operations.sort_boundary_points(boundary_coords)

Sort boundary points efficiently using spatial relationships.

Uses a KD-tree for efficient nearest neighbor queries to create a continuous path along the boundary points. This is useful for creating side faces that close a terrain mesh into a solid object.

Parameters:

boundary_coords – List of (y, x) coordinate tuples representing boundary points

Returns:

Sorted boundary points forming a continuous path around the perimeter

Return type:

list

Sort boundary points for ordered traversal.

Rectangle Edge Generation

src.terrain.mesh_operations.generate_rectangle_edge_pixels(dem_shape, edge_sample_spacing=1.0)

Generate boundary pixel coordinates by sampling rectangle edges.

This creates ordered (y, x) pixel coordinates around the DEM boundary, forming a simple rectangle. This approach is much faster than morphological boundary detection and works well for rectangular DEMs.

Algorithm: Sample the rectangle boundary edges at given spacing, tracing counterclockwise: top edge → right edge → bottom edge → left edge

Parameters:

• dem_shape (tuple) – DEM shape (height, width)

• edge_sample_spacing (float) – Pixel spacing for edge sampling (default: 1.0)

Returns:

Ordered list of (y, x) pixel coordinates forming the rectangle boundary

Return type:

list

Generate edge pixels using rectangle sampling (fast alternative to morphological).

Performance: ~150x faster than morphological edge detection.

src.terrain.mesh_operations.generate_rectangle_edge_vertices(dem_shape, dem_data, original_transform, transforms_list, edge_sample_spacing=1.0, base_depth=-0.2)

Generate boundary vertices by sampling rectangle edges and applying geographic transforms.

This approach leverages the same transform pipeline used for the DEM to create naturally smooth, curved edges that perfectly match the geographic projection.

Algorithm: 1. Sample the rectangle boundary in original DEM pixel space 2. For each edge vertex, apply the sequence of geographic transforms 3. Creates BOTH surface vertices (at DEM elevation) AND base vertices (at base_depth) 4. Generates quad faces forming vertical walls (“skirt”) around the terrain edge

Vertex layout: - Indices 0 to n-1: Surface vertices (at DEM elevation) - Indices n to 2n-1: Base vertices (at base_depth, same x,y as surface)

Parameters:

• dem_shape (tuple) – Original DEM shape (height, width)

• dem_data (np.ndarray) – Original DEM data array

• original_transform (Affine) – Original affine transform (pixel → geographic)

• transforms_list (list) – List of transform functions to apply sequentially

• edge_sample_spacing (float) – Pixel spacing for edge sampling (default: 1.0)

• base_depth (float) – Z-coordinate for base vertices (default: -0.2)

Returns:

(boundary_vertices, boundary_faces) where:

• boundary_vertices: (2*n)x3 array of vertex positions (surface + base)

• boundary_faces: List of quad faces forming vertical walls

Return type:

tuple

Generate edge vertices from rectangle edge pixels.

src.terrain.mesh_operations.generate_transform_aware_rectangle_edges(terrain, coord_to_index, edge_sample_spacing=1.0)

Generate rectangle edge pixels by sampling original DEM perimeter and mapping through transform pipeline.

This function solves the NaN margin problem by sampling edges at the original DEM resolution (where all perimeter pixels are valid) and mapping them through the transform pipeline to final mesh coordinates.

Uses affine transforms to map coordinates:

original pixel → geographic → final transformed pixel

Parameters:

• terrain – Terrain object with dem_shape, dem_transform, data_layers

• coord_to_index – Dict mapping (y, x) final pixels to vertex indices

• edge_sample_spacing – Pixel spacing for edge sampling at original resolution (default 1.0)

Returns:

Edge pixel coordinates in final transformed space,

filtered to only valid mesh vertices

Return type:

list

Raises:

ValueError – If terrain is None or lacks required transform data

Generate rectangle edges with geographic transform awareness.

src.terrain.mesh_operations.generate_transform_aware_rectangle_edges_fractional(terrain, edge_sample_spacing=1.0)

Generate rectangle edge vertices with FRACTIONAL coordinates preserving projection curvature.

Unlike generate_transform_aware_rectangle_edges() which rounds to integers and filters to existing mesh vertices, this function returns the true fractional coordinates from the non-linear projection transformation.

This preserves the curved boundary that results from: - WGS84 → UTM Transverse Mercator projection (non-linear, causes curvature) - Horizontal flip transform - Downsampling

Parameters:

• terrain – Terrain object with dem_shape, dem_transform, data_layers

• edge_sample_spacing – Pixel spacing for edge sampling at original resolution (default 1.0)

Returns:

Fractional edge coordinates in final mesh space.

These coordinates preserve the true curved boundary and may extend slightly beyond the integer grid bounds.

Return type:

list of (y, x) tuples

Raises:

ValueError – If terrain is None or lacks required transform data

Generate fractional rectangle edges preserving projection curvature.

Fractional edges: Preserve subtle curves from WGS84→UTM transformation.

src.terrain.mesh_operations.diagnose_rectangle_edge_coverage(dem_shape, coord_to_index)

Diagnose how well rectangle edge sampling will work for this DEM.

Checks what percentage of the rectangle perimeter has valid mesh vertices. Helps determine if rectangle-edge sampling is appropriate for this dataset.

Parameters:

• dem_shape (tuple) – DEM shape (height, width)

• coord_to_index (dict) – Mapping from (y, x) to vertex indices

Returns:

Diagnostic information including:

• total_edge_pixels: Total pixels on rectangle perimeter

• valid_edge_pixels: How many have valid mesh vertices

• coverage_percent: Percentage of edge that’s valid

• edge_validity: Per-edge breakdown (top, right, bottom, left)

• recommendation: Whether to use rectangle edges or morphological

Return type:

dict

Diagnostic function to analyze rectangle edge coverage.

Performance Optimizations

Numba JIT compilation:

• Face generation: ~10x speedup

• Boundary processing: ~5x speedup

• Requires: numba package (pip install numba)

• Automatic fallback to Python if numba unavailable

Typical performance:

Operation	Array Size	Time
Generate vertices Generate faces (Numba) Find boundary (morph) Find boundary (rectangle) Catmull-Rom fit	1M points 2M triangles 4096² 4096² 1000 pts, sub=20	~50ms ~200ms ~500ms ~3ms (150x) ~100ms

Memory usage:

• Vertices: 12 bytes per vertex (3 float32s)

• Faces: 12 bytes per triangle (3 int32s)

• 1M vertices + 2M triangles ≈ 36MB

See Also

• Core Module - Uses these functions for mesh creation

• Blender Integration Module - Applies vertex colors to generated meshes

• Combined Render: Full-Featured Example - Example usage with two-tier edges