"""CPU point-cloud understanding for existing reconstruction outputs.
|
|
This B capability uses Open3D for point-cloud reading, downsampling, ground
|
plane fitting, and mesh reconstruction. ``geoai.masks_to_vector`` converts
|
the elevated-object raster to GeoJSON. It does not run photo-based SfM/MVS or
|
claim semantic object labels.
|
"""
|
|
from __future__ import annotations
|
|
import argparse
|
import json
|
import math
|
import time
|
from dataclasses import dataclass
|
from datetime import UTC, datetime
|
from importlib import metadata as importlib_metadata
|
from pathlib import Path
|
from typing import Any
|
|
import geopandas as gpd
|
import laspy
|
import numpy as np
|
import open3d as o3d
|
import pandas as pd
|
import rasterio
|
from PIL import Image
|
from rasterio.features import shapes
|
from rasterio.transform import from_origin
|
from scipy import ndimage
|
from scipy.spatial import cKDTree
|
from shapely.geometry import shape
|
|
|
SUPPORTED_SUFFIXES = {".ply", ".pcd", ".xyz", ".xyzn", ".xyzrgb", ".las", ".laz"}
|
|
# LAS class codes keep the result useful outside this workbench. The labels are
|
# deliberately conservative: uncertain elevated points remain unclassified.
|
SEMANTIC_CLASSES = {
|
1: {"key": "other_unknown", "label": "Other / unknown", "color": (128, 128, 128)},
|
2: {"key": "ground", "label": "Ground", "color": (151, 111, 51)},
|
5: {"key": "vegetation", "label": "Vegetation", "color": (59, 163, 87)},
|
6: {"key": "building_structure", "label": "Building / structure", "color": (224, 115, 55)},
|
15: {"key": "pole_tower", "label": "Pole / tower candidate", "color": (149, 89, 210)},
|
16: {"key": "power_line", "label": "Power-line candidate", "color": (231, 196, 61)},
|
}
|
SEMANTIC_PRIORITY = {1: 1, 2: 2, 5: 3, 6: 4, 15: 5, 16: 6}
|
|
|
@dataclass(frozen=True)
|
class PointCloudData:
|
points: np.ndarray
|
colors: np.ndarray
|
source_has_rgb: bool
|
|
|
def _normalise_rgb(colors: np.ndarray) -> np.ndarray:
|
values = np.asarray(colors, dtype=np.float64)
|
if values.size == 0:
|
return values.reshape((-1, 3))
|
high = float(np.nanpercentile(values, 99.5))
|
divisor = 65_535.0 if high > 255 else 255.0
|
return np.clip(values / divisor, 0.0, 1.0)
|
|
|
def read_point_cloud_data(path: Path) -> PointCloudData:
|
if path.suffix.lower() in {".las", ".laz"}:
|
source = laspy.read(path)
|
points = np.column_stack((source.x, source.y, source.z)).astype(np.float64)
|
dimensions = set(source.point_format.dimension_names)
|
has_rgb = {"red", "green", "blue"}.issubset(dimensions)
|
colors = _normalise_rgb(np.column_stack((source.red, source.green, source.blue))) if has_rgb else np.full((len(points), 3), 0.72)
|
else:
|
cloud = o3d.io.read_point_cloud(str(path))
|
points = np.asarray(cloud.points, dtype=np.float64)
|
has_rgb = cloud.has_colors()
|
colors = np.asarray(cloud.colors, dtype=np.float64) if has_rgb else np.full((len(points), 3), 0.72)
|
if len(points) < 50:
|
raise ValueError(f"Point cloud needs at least 50 finite points: {path.name}")
|
finite = np.isfinite(points).all(axis=1) & np.isfinite(colors).all(axis=1)
|
if not finite.any():
|
raise ValueError(f"Point cloud has no finite points: {path.name}")
|
return PointCloudData(points=points[finite], colors=np.clip(colors[finite], 0.0, 1.0), source_has_rgb=has_rgb)
|
|
|
def read_point_cloud(path: Path) -> o3d.geometry.PointCloud:
|
data = read_point_cloud_data(path)
|
cloud = o3d.geometry.PointCloud()
|
cloud.points = o3d.utility.Vector3dVector(data.points)
|
cloud.colors = o3d.utility.Vector3dVector(data.colors)
|
return cloud
|
|
|
def voxel_downsample_data(data: PointCloudData, voxel_size: float) -> PointCloudData:
|
"""Keep one deterministic representative per voxel while retaining RGB."""
|
origin = data.points.min(axis=0)
|
cells = np.floor((data.points - origin) / voxel_size).astype(np.int64)
|
_, indices = np.unique(cells, axis=0, return_index=True)
|
indices.sort()
|
return PointCloudData(points=data.points[indices], colors=data.colors[indices], source_has_rgb=data.source_has_rgb)
|
|
|
def _plane_height(plane: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
|
a, b, c, d = plane
|
return -(a * x + b * y + d) / c
|
|
|
def fit_ground(cloud: o3d.geometry.PointCloud, distance_threshold: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
plane, inliers = cloud.segment_plane(distance_threshold=distance_threshold, ransac_n=3, num_iterations=1_000)
|
plane_values = np.asarray(plane, dtype=float)
|
if abs(plane_values[2]) < 0.7:
|
raise ValueError("Dominant RANSAC plane is too steep to be used as ground.")
|
if plane_values[2] < 0:
|
plane_values *= -1
|
values = np.asarray(cloud.points)
|
ground_z = _plane_height(plane_values, values[:, 0], values[:, 1])
|
height_above_ground = values[:, 2] - ground_z
|
return plane_values, np.asarray(inliers, dtype=int), height_above_ground
|
|
|
def orient_ground_up(cloud: o3d.geometry.PointCloud, distance_threshold: float, up_axis: str) -> tuple[o3d.geometry.PointCloud, dict[str, Any]]:
|
"""Optionally rotate a source-local cloud so its dominant ground plane uses Z as up."""
|
if up_axis == "z":
|
return cloud, {"method": "source_z_axis", "applied": False}
|
plane, _ = cloud.segment_plane(distance_threshold=distance_threshold, ransac_n=3, num_iterations=1_000)
|
normal = np.asarray(plane[:3], dtype=float)
|
normal /= np.linalg.norm(normal)
|
if normal[2] < 0:
|
normal *= -1
|
target = np.array([0.0, 0.0, 1.0])
|
axis = np.cross(normal, target)
|
axis_norm = np.linalg.norm(axis)
|
angle = math.atan2(axis_norm, float(np.dot(normal, target)))
|
if axis_norm > 1e-12:
|
rotation = o3d.geometry.get_rotation_matrix_from_axis_angle(axis / axis_norm * angle)
|
cloud = o3d.geometry.PointCloud(cloud)
|
cloud.rotate(rotation, center=(0.0, 0.0, 0.0))
|
return cloud, {
|
"method": "auto_align_dominant_plane_to_z",
|
"applied": bool(axis_norm > 1e-12),
|
"source_plane_normal": [round(float(value), 8) for value in normal],
|
"rotation_degrees": round(math.degrees(angle), 4),
|
}
|
|
|
def _grid_bounds(values: np.ndarray, cell_size: float) -> tuple[float, float, int, int]:
|
minimum = values[:, :2].min(axis=0)
|
maximum = values[:, :2].max(axis=0)
|
width = max(1, int(math.ceil((maximum[0] - minimum[0]) / cell_size)) + 1)
|
height = max(1, int(math.ceil((maximum[1] - minimum[1]) / cell_size)) + 1)
|
if width > 4_096 or height > 4_096:
|
raise ValueError("Point-cloud extent and cell size would create a raster larger than 4096 by 4096.")
|
return float(minimum[0]), float(maximum[1]), width, height
|
|
|
def estimate_local_ground(values: np.ndarray, cell_size: float) -> tuple[np.ndarray, np.ndarray, Any]:
|
"""Estimate a conservative local ground surface from low points in each XY cell.
|
|
This is intentionally a CPU rule baseline, not a bare-earth DEM algorithm.
|
It behaves more safely than one global plane on a sloped power-line corridor.
|
"""
|
origin_x, origin_y, width, height = _grid_bounds(values, cell_size)
|
transform = from_origin(origin_x, origin_y, cell_size, cell_size)
|
columns = np.clip(((values[:, 0] - origin_x) / cell_size).astype(int), 0, width - 1)
|
rows = np.clip(((origin_y - values[:, 1]) / cell_size).astype(int), 0, height - 1)
|
surface = np.full((height, width), np.inf, dtype=np.float64)
|
np.minimum.at(surface, (rows, columns), values[:, 2])
|
valid = surface != np.inf
|
if valid.sum() < 8:
|
raise ValueError("Too few occupied XY cells to estimate a local ground surface.")
|
nearest = ndimage.distance_transform_edt(~valid, return_distances=False, return_indices=True)
|
filled = surface[tuple(nearest)]
|
# A small median window suppresses isolated low outliers while retaining
|
# terrain changes at the metre-scale used by this first baseline.
|
smoothed = ndimage.median_filter(filled, size=3, mode="nearest")
|
point_ground = smoothed[rows, columns]
|
return point_ground, smoothed, transform
|
|
|
def _cell_metrics(values: np.ndarray, cell_size: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
"""Return per-point count/min/max metrics of its horizontal neighbourhood."""
|
minimum = values[:, :2].min(axis=0)
|
cells = np.floor((values[:, :2] - minimum) / cell_size).astype(np.int64)
|
_, inverse = np.unique(cells, axis=0, return_inverse=True)
|
counts = np.bincount(inverse)
|
mins = np.full(len(counts), np.inf)
|
maxs = np.full(len(counts), -np.inf)
|
np.minimum.at(mins, inverse, values[:, 2])
|
np.maximum.at(maxs, inverse, values[:, 2])
|
return counts[inverse], mins[inverse], maxs[inverse]
|
|
|
def continuous_wire_candidates(
|
values: np.ndarray,
|
heights: np.ndarray,
|
green_dominant: np.ndarray,
|
compact_counts: np.ndarray,
|
compact_span: np.ndarray,
|
) -> np.ndarray:
|
"""Keep sparse high line candidates for review without using RGB as a veto."""
|
cell_size = 1.0
|
minimum = values[:, :2].min(axis=0)
|
cells = np.floor((values[:, :2] - minimum) / cell_size).astype(int)
|
width = int(cells[:, 0].max()) + 1
|
height = int(cells[:, 1].max()) + 1
|
flat = cells[:, 1] * width + cells[:, 0]
|
cell_count = width * height
|
# LiDAR RGB is not a material classifier. In the supplied corridor scan
|
# conductors and tower parts are often green, so RGB must not reject them.
|
candidate = (heights >= 4.0) & (compact_span <= 1.20) & (compact_counts <= 100)
|
cell_candidate = np.zeros(cell_count, dtype=bool)
|
np.logical_or.at(cell_candidate, flat, candidate)
|
# Diagonal conductors may leave a one-cell gap after voxel sampling.
|
connected = ndimage.binary_dilation(cell_candidate.reshape(height, width), iterations=2)
|
labels, component_count = ndimage.label(connected, structure=np.ones((3, 3), dtype=np.uint8))
|
point_components = labels[cells[:, 1], cells[:, 0]]
|
accepted = np.zeros(component_count + 1, dtype=bool)
|
for component in range(1, component_count + 1):
|
original_cells = np.flatnonzero((labels.ravel() == component) & cell_candidate)
|
if len(original_cells) < 3:
|
continue
|
xy = np.column_stack((original_cells % width, original_cells // width)).astype(float)
|
centered = xy - xy.mean(axis=0)
|
eigenvalues = np.linalg.eigvalsh(centered.T @ centered)
|
if eigenvalues[-1] <= 0:
|
continue
|
linearity = 1.0 - eigenvalues[0] / eigenvalues[-1]
|
line_length = float(np.ptp(centered @ np.linalg.eigh(centered.T @ centered)[1][:, -1]))
|
if linearity >= 0.94 and line_length >= 5.0:
|
accepted[component] = True
|
# Broken returns and oblique viewing often split one conductor into many
|
# short components. Preserve all thin high samples as *candidates*; the
|
# continuity gate remains a diagnostic for later model/label work.
|
return candidate
|
|
|
def local_shape_features(values: np.ndarray, neighbours: int = 16, batch_size: int = 100_000) -> tuple[np.ndarray, np.ndarray]:
|
"""Measure local linearity and principal-axis verticality without RGB.
|
|
The work is batched so multi-million-point LAS inputs stay bounded in RAM.
|
"""
|
tree = cKDTree(values)
|
linearity = np.zeros(len(values), dtype=np.float32)
|
principal_verticality = np.zeros(len(values), dtype=np.float32)
|
for start in range(0, len(values), batch_size):
|
stop = min(start + batch_size, len(values))
|
_, indices = tree.query(values[start:stop], k=min(neighbours, len(values)), workers=-1)
|
neighbours_xyz = values[np.atleast_2d(indices)]
|
centred = neighbours_xyz - neighbours_xyz.mean(axis=1, keepdims=True)
|
covariance = np.einsum("nij,nik->njk", centred, centred) / max(neighbours_xyz.shape[1] - 1, 1)
|
eigenvalues, eigenvectors = np.linalg.eigh(covariance)
|
largest = np.maximum(eigenvalues[:, 2], 1e-9)
|
linearity[start:stop] = np.clip((eigenvalues[:, 2] - eigenvalues[:, 1]) / largest, 0.0, 1.0)
|
principal_verticality[start:stop] = np.abs(eigenvectors[:, 2, 2])
|
return linearity, principal_verticality
|
|
|
def expand_pole_tower_components(values: np.ndarray, heights: np.ndarray, pole_seeds: np.ndarray) -> np.ndarray:
|
"""Promote the horizontal and diagonal members around a tower's vertical seeds.
|
|
Lattice towers are not locally vertical at every point. This object-level
|
step is deliberately limited to compact groups of high-confidence vertical
|
members, so a long vegetation corridor is not promoted wholesale.
|
"""
|
expanded = pole_seeds.copy()
|
seed_indices = np.flatnonzero(pole_seeds)
|
if len(seed_indices) < 6:
|
return expanded
|
seed_xy = values[seed_indices, :2]
|
pairs = cKDTree(seed_xy).query_pairs(r=4.5, output_type="ndarray")
|
parent = np.arange(len(seed_indices))
|
|
def find(index: int) -> int:
|
while parent[index] != index:
|
parent[index] = parent[parent[index]]
|
index = int(parent[index])
|
return index
|
|
for left, right in pairs:
|
root_left, root_right = find(int(left)), find(int(right))
|
if root_left != root_right:
|
parent[root_right] = root_left
|
groups: dict[int, list[int]] = {}
|
for index in range(len(seed_indices)):
|
groups.setdefault(find(index), []).append(index)
|
for members in groups.values():
|
if len(members) < 8:
|
continue
|
indices = seed_indices[np.asarray(members)]
|
seed_heights = heights[indices]
|
if float(seed_heights.max() - seed_heights.min()) < 5.0:
|
continue
|
center = np.median(values[indices, :2], axis=0)
|
distances = np.linalg.norm(values[indices, :2] - center, axis=1)
|
if float(np.percentile(distances, 95)) > 4.5:
|
continue
|
radius = min(4.0, max(2.0, float(np.percentile(distances, 90)) + 1.0))
|
lower = max(1.5, float(np.percentile(seed_heights, 5)) - 0.75)
|
upper = float(np.percentile(seed_heights, 95)) + 0.75
|
nearby = np.linalg.norm(values[:, :2] - center, axis=1) <= radius
|
expanded |= nearby & (heights >= lower) & (heights <= upper)
|
return expanded
|
|
|
def classify_semantic_points(values: np.ndarray, colors: np.ndarray, heights: np.ndarray, voxel_size: float) -> np.ndarray:
|
"""Classify points with transparent RGB/geometry rules on CPU.
|
|
The narrow wire and pole rules deliberately require strong evidence. Points
|
that do not satisfy a class remain ``other_unknown`` for user review.
|
"""
|
classes = np.full(len(values), 1, dtype=np.uint8)
|
ground_limit = max(0.30, voxel_size * 1.5)
|
ground = heights <= ground_limit
|
classes[ground] = 2
|
|
red, green, blue = colors.T
|
green_dominant = (green > red * 1.08) & (green > blue * 1.05) & (green > 0.16)
|
compact_counts, compact_min, compact_max = _cell_metrics(values, max(1.0, voxel_size * 5))
|
compact_span = compact_max - compact_min
|
tower_counts, tower_min, tower_max = _cell_metrics(values, max(1.5, voxel_size * 7.5))
|
tower_span = tower_max - tower_min
|
elevated = heights > 1.0
|
|
local_linearity, local_verticality = local_shape_features(values)
|
# Conductors are horizontally linear; colour must not suppress this signal.
|
wire = (heights >= 4.0) & (local_linearity >= 0.86) & (local_verticality <= 0.40) & (compact_counts <= 180)
|
# Pole/tower members are vertically linear. This excludes most foliage,
|
# which has no stable local principal direction.
|
pole_seeds = elevated & (heights >= 3.0) & (local_linearity >= 0.78) & (local_verticality >= 0.72) & (tower_counts <= 500) & ~wire
|
pole = expand_pole_tower_components(values, heights, pole_seeds) & ~wire
|
# Dense locally planar elevated surfaces are structure candidates. This
|
# avoids labelling rough foliage as a building merely from its height.
|
structure = elevated & (compact_span <= 0.65) & (compact_counts >= 15) & ~green_dominant & ~wire & ~pole
|
vegetation = elevated & ~wire & ~pole & ~structure & (green_dominant | (compact_span >= 1.15))
|
classes[vegetation] = 5
|
classes[structure] = 6
|
classes[pole] = 15
|
classes[wire] = 16
|
return classes
|
|
|
def semantic_raster(values: np.ndarray, classes: np.ndarray, cell_size: float) -> tuple[np.ndarray, Any]:
|
origin_x, origin_y, width, height = _grid_bounds(values, cell_size)
|
transform = from_origin(origin_x, origin_y, cell_size, cell_size)
|
columns = np.clip(((values[:, 0] - origin_x) / cell_size).astype(int), 0, width - 1)
|
rows = np.clip(((origin_y - values[:, 1]) / cell_size).astype(int), 0, height - 1)
|
raster = np.zeros((height, width), dtype=np.uint8)
|
priorities = np.asarray([SEMANTIC_PRIORITY.get(int(value), 0) for value in classes], dtype=np.uint8)
|
flat = rows * width + columns
|
winner = np.zeros(height * width, dtype=np.uint8)
|
np.maximum.at(winner, flat, priorities)
|
priority_to_code = {priority: code for code, priority in SEMANTIC_PRIORITY.items()}
|
for priority, code in priority_to_code.items():
|
raster.flat[winner == priority] = code
|
return raster, transform
|
|
|
def write_semantic_preview(raster: np.ndarray, path: Path) -> None:
|
preview = np.full((*raster.shape, 3), 244, dtype=np.uint8)
|
for code, details in SEMANTIC_CLASSES.items():
|
preview[raster == code] = details["color"]
|
Image.fromarray(preview).save(path)
|
|
|
def write_semantic_vectors(raster: np.ndarray, transform: Any, path: Path) -> int:
|
records: list[dict[str, Any]] = []
|
for geometry, value in shapes(raster, mask=raster > 0, transform=transform):
|
code = int(value)
|
details = SEMANTIC_CLASSES.get(code)
|
if not details or code == 1:
|
continue
|
polygon = shape(geometry)
|
if polygon.area <= 0:
|
continue
|
records.append({"class_key": details["key"], "class_label": details["label"], "las_class_code": code, "area_local_units2": round(float(polygon.area), 3), "geometry": polygon})
|
frame = gpd.GeoDataFrame(records, geometry="geometry")
|
if frame.empty:
|
frame = gpd.GeoDataFrame({"class_key": [], "class_label": [], "las_class_code": [], "area_local_units2": []}, geometry=[])
|
frame.to_file(path, driver="GeoJSON")
|
return len(frame)
|
|
|
def write_semantic_las(values: np.ndarray, colors: np.ndarray, classes: np.ndarray, path: Path) -> None:
|
header = laspy.LasHeader(point_format=3, version="1.2")
|
header.scales = np.array([0.001, 0.001, 0.001])
|
header.offsets = np.floor(values.min(axis=0))
|
output = laspy.LasData(header)
|
output.x, output.y, output.z = values.T
|
output.red, output.green, output.blue = (np.clip(colors, 0.0, 1.0) * 65535).astype(np.uint16).T
|
output.classification = classes
|
output.write(path)
|
|
|
def rasterize(
|
values: np.ndarray,
|
heights: np.ndarray,
|
cell_size: float,
|
elevated_threshold: float,
|
) -> tuple[np.ndarray, np.ndarray, np.ndarray, Any]:
|
origin_x, origin_y, width, height = _grid_bounds(values, cell_size)
|
transform = from_origin(origin_x, origin_y, cell_size, cell_size)
|
dsm = np.full((height, width), np.nan, dtype=np.float32)
|
canopy_height = np.full((height, width), np.nan, dtype=np.float32)
|
columns = np.clip(((values[:, 0] - origin_x) / cell_size).astype(int), 0, width - 1)
|
rows = np.clip(((origin_y - values[:, 1]) / cell_size).astype(int), 0, height - 1)
|
for row, column, z, above_ground in zip(rows, columns, values[:, 2], heights, strict=True):
|
dsm[row, column] = z if not np.isfinite(dsm[row, column]) else max(dsm[row, column], z)
|
canopy_height[row, column] = above_ground if not np.isfinite(canopy_height[row, column]) else max(canopy_height[row, column], above_ground)
|
labels = np.where(np.isfinite(canopy_height) & (canopy_height >= elevated_threshold), 255, 0).astype(np.uint8)
|
return dsm, canopy_height, labels, transform
|
|
|
def _colored_preview(values: np.ndarray, labels: np.ndarray, path: Path) -> None:
|
finite = np.isfinite(values)
|
preview = np.zeros((*values.shape, 3), dtype=np.uint8)
|
if finite.any():
|
low, high = np.percentile(values[finite], [2, 98])
|
normalized = np.nan_to_num(np.clip((values - low) / max(high - low, 1e-6), 0, 1), nan=0.0)
|
preview[..., 0] = (30 + 180 * normalized).astype(np.uint8)
|
preview[..., 1] = (65 + 150 * (1 - normalized)).astype(np.uint8)
|
preview[..., 2] = (205 - 150 * normalized).astype(np.uint8)
|
preview[labels > 0] = [230, 90, 45]
|
Image.fromarray(preview).save(path)
|
|
|
def _fallback_vectors(labels: np.ndarray, transform: Any) -> gpd.GeoDataFrame:
|
records: list[dict[str, Any]] = []
|
for geometry, value in shapes(labels, mask=labels.astype(bool), transform=transform):
|
if int(value) != 255:
|
continue
|
polygon = shape(geometry)
|
if polygon.area <= 0:
|
continue
|
records.append({"class_key": "elevated_surface", "geometry": polygon})
|
return gpd.GeoDataFrame(records, geometry="geometry")
|
|
|
def vectorize(labels: np.ndarray, transform: Any, raster_path: Path, vector_path: Path) -> tuple[gpd.GeoDataFrame, str]:
|
fallback = _fallback_vectors(labels, transform)
|
vectorizer = "rasterio.features.shapes fallback"
|
try:
|
from geoai import masks_to_vector
|
|
frame = masks_to_vector(str(raster_path), min_object_area=1, simplify_tolerance=0.0)
|
if not frame.empty:
|
frame = frame[["geometry"]].copy()
|
frame["class_key"] = "elevated_surface"
|
geoai_area = float(frame.geometry.area.sum())
|
fallback_area = float(fallback.geometry.area.sum())
|
if len(frame) == len(fallback) and fallback_area and 0.98 <= geoai_area / fallback_area <= 1.02:
|
fallback = frame
|
vectorizer = "geoai.masks_to_vector"
|
else:
|
vectorizer = "geoai.masks_to_vector + rasterio completeness repair"
|
except Exception:
|
pass
|
fallback = fallback.reset_index(drop=True)
|
fallback["feature_id"] = np.arange(1, len(fallback) + 1)
|
fallback["area_local_units2"] = fallback.geometry.area.round(3)
|
fallback.to_file(vector_path, driver="GeoJSON")
|
return fallback, vectorizer
|
|
|
def reconstruct_mesh(cloud: o3d.geometry.PointCloud, output_path: Path, voxel_size: float) -> int:
|
working = cloud.voxel_down_sample(max(voxel_size, 0.02))
|
if len(working.points) < 50:
|
return 0
|
working.estimate_normals(o3d.geometry.KDTreeSearchParamHybrid(radius=max(voxel_size * 4, 0.2), max_nn=30))
|
try:
|
mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(working, max(voxel_size * 4, 0.4))
|
except RuntimeError:
|
return 0
|
if len(mesh.triangles) and o3d.io.write_triangle_mesh(str(output_path), mesh, write_ascii=False):
|
return len(mesh.triangles)
|
return 0
|
|
|
def process_point_cloud(
|
input_path: Path,
|
output_dir: Path,
|
voxel_size: float = 0.2,
|
ground_distance: float = 0.15,
|
elevated_threshold: float = 0.75,
|
ground_up_axis: str = "z",
|
) -> dict[str, Any]:
|
if voxel_size <= 0 or ground_distance <= 0 or elevated_threshold <= 0:
|
raise ValueError("voxel size, ground distance, and elevated threshold must be positive.")
|
started = time.perf_counter()
|
source_data = read_point_cloud_data(input_path)
|
original_points = len(source_data.points)
|
downsampled_data = voxel_downsample_data(source_data, voxel_size)
|
if len(downsampled_data.points) < 50:
|
raise ValueError("Voxel downsampling retained fewer than 50 points; use a smaller voxel size.")
|
downsampled = o3d.geometry.PointCloud()
|
downsampled.points = o3d.utility.Vector3dVector(downsampled_data.points)
|
downsampled.colors = o3d.utility.Vector3dVector(downsampled_data.colors)
|
downsampled, alignment = orient_ground_up(downsampled, ground_distance, ground_up_axis)
|
values = np.asarray(downsampled.points)
|
# Open3D exposes a mutable view here. Keep an owned RGB copy before the
|
# classified point cloud overwrites its colours with display labels.
|
observed_colors = np.asarray(downsampled.colors).copy()
|
ground_plane_available = True
|
try:
|
plane, inliers, _ = fit_ground(downsampled, ground_distance)
|
except (RuntimeError, ValueError):
|
ground_plane_available = False
|
plane = np.array([0.0, 0.0, 1.0, -float(np.median(values[:, 2]))])
|
inliers = np.array([], dtype=int)
|
local_ground, _, _ = estimate_local_ground(values, max(1.0, voxel_size * 5))
|
heights = values[:, 2] - local_ground
|
inliers = np.flatnonzero(heights <= max(ground_distance, voxel_size * 1.5))
|
dsm, height_raster, labels, transform = rasterize(values, heights, voxel_size, elevated_threshold)
|
output_dir.mkdir(parents=True, exist_ok=True)
|
stem = input_path.stem
|
classified = output_dir / f"{stem}.classified.ply"
|
semantic_preview_cloud = output_dir / f"{stem}.semantic-preview.ply"
|
annotation_source_cloud = output_dir / f"{stem}.semantic-annotation-source.ply"
|
semantic_las = output_dir / f"{stem}.semantic-classified.las"
|
semantic_raster_path = output_dir / f"{stem}.semantic-classes.tif"
|
semantic_preview_path = output_dir / f"{stem}.semantic-classes.preview.png"
|
semantic_vector_path = output_dir / f"{stem}.semantic-footprints.geojson"
|
semantic_summary_path = output_dir / f"{stem}.semantic-summary.csv"
|
mesh = output_dir / f"{stem}.reconstruction.ply"
|
dsm_path = output_dir / f"{stem}.dsm.tif"
|
height_path = output_dir / f"{stem}.height-above-ground.tif"
|
label_path = output_dir / f"{stem}.elevated-surface.tif"
|
preview_path = output_dir / f"{stem}.dsm.preview.png"
|
vector_path = output_dir / f"{stem}.elevated-footprints.geojson"
|
summary_path = output_dir / f"{stem}.summary.csv"
|
semantic_classes = classify_semantic_points(values, observed_colors, heights, voxel_size)
|
semantic_colors = np.asarray([np.asarray(SEMANTIC_CLASSES[int(code)]["color"], dtype=float) / 255.0 for code in semantic_classes])
|
downsampled.colors = o3d.utility.Vector3dVector(semantic_colors)
|
o3d.io.write_point_cloud(str(classified), downsampled, write_ascii=False)
|
# The console now renders the entire processed point set, so sparse
|
# conductors and all other classes are retained without preview sampling.
|
preview_indices = np.arange(len(values), dtype=np.int64)
|
preview_cloud = downsampled.select_by_index(preview_indices.tolist())
|
o3d.io.write_point_cloud(str(semantic_preview_cloud), preview_cloud, write_ascii=False)
|
# Annotation must retain observed RGB. The semantic preview uses rule
|
# colours for review only, and must never become a feature leak for a later
|
# supervised model.
|
annotation_cloud = o3d.geometry.PointCloud()
|
annotation_cloud.points = o3d.utility.Vector3dVector(values[preview_indices])
|
annotation_cloud.colors = o3d.utility.Vector3dVector(observed_colors[preview_indices])
|
o3d.io.write_point_cloud(str(annotation_source_cloud), annotation_cloud, write_ascii=False)
|
write_semantic_las(values, observed_colors, semantic_classes, semantic_las)
|
semantic_labels, semantic_transform = semantic_raster(values, semantic_classes, voxel_size)
|
with rasterio.open(semantic_raster_path, "w", driver="GTiff", height=semantic_labels.shape[0], width=semantic_labels.shape[1], count=1, dtype="uint8", transform=semantic_transform, nodata=0) as destination:
|
destination.write(semantic_labels, 1)
|
write_semantic_preview(semantic_labels, semantic_preview_path)
|
semantic_footprints = write_semantic_vectors(semantic_labels, semantic_transform, semantic_vector_path)
|
semantic_counts = pd.DataFrame([
|
{
|
"las_class_code": code,
|
"class_key": details["key"],
|
"class_label": details["label"],
|
"point_count": int((semantic_classes == code).sum()),
|
"point_ratio": round(float((semantic_classes == code).mean()), 6),
|
}
|
for code, details in SEMANTIC_CLASSES.items()
|
])
|
semantic_counts.to_csv(semantic_summary_path, index=False, encoding="utf-8-sig")
|
# Alpha-shape meshing is useful for small geometry demos but can dominate a
|
# semantic-classification run without improving its class labels.
|
triangles = reconstruct_mesh(downsampled, mesh, voxel_size) if len(values) <= 250_000 else 0
|
for path, raster, dtype, nodata in ((dsm_path, dsm, "float32", -9999.0), (height_path, height_raster, "float32", -9999.0), (label_path, labels, "uint8", 0)):
|
write_values = np.where(np.isfinite(raster), raster, nodata).astype(dtype)
|
with rasterio.open(path, "w", driver="GTiff", height=raster.shape[0], width=raster.shape[1], count=1, dtype=dtype, transform=transform, nodata=nodata) as dst:
|
dst.write(write_values, 1)
|
_colored_preview(dsm, labels, preview_path)
|
footprints, vectorizer = vectorize(labels, transform, label_path, vector_path)
|
elevated_count = int((heights >= elevated_threshold).sum())
|
summary = pd.DataFrame([{
|
"input_file": input_path.name,
|
"original_points": original_points,
|
"downsampled_points": len(values),
|
"ground_inliers": len(inliers),
|
"elevated_points": elevated_count,
|
"elevated_footprints": len(footprints),
|
"mesh_triangles": triangles,
|
"mesh_skipped_for_large_cloud": bool(len(values) > 250_000),
|
"raster_width": dsm.shape[1],
|
"raster_height": dsm.shape[0],
|
"cell_size_local_units": voxel_size,
|
}])
|
summary.to_csv(summary_path, index=False, encoding="utf-8-sig")
|
return {
|
"file": input_path.name,
|
"original_points": original_points,
|
"downsampled_points": len(values),
|
"ground_inliers": len(inliers),
|
"elevated_points": elevated_count,
|
"elevated_point_ratio": round(elevated_count / len(values), 5),
|
"elevated_footprint_count": len(footprints),
|
"mesh_triangles": triangles,
|
"raster_width": dsm.shape[1],
|
"raster_height": dsm.shape[0],
|
"coordinate_basis": "local_point_cloud_coordinates_ground_aligned" if alignment["applied"] else "local_point_cloud_coordinates",
|
"crs": None,
|
"classified_point_cloud": classified.name,
|
"semantic_preview_point_cloud": semantic_preview_cloud.name,
|
"semantic_annotation_source_point_cloud": annotation_source_cloud.name,
|
"semantic_annotation_source_kind": "complete processed RGB/XYZ point set after voxel sampling; no semantic display colours",
|
"semantic_preview_points": int(len(preview_indices)),
|
"semantic_classified_las": semantic_las.name,
|
"semantic_raster_file": semantic_raster_path.name,
|
"semantic_preview_file": semantic_preview_path.name,
|
"semantic_vector_file": semantic_vector_path.name,
|
"semantic_summary_file": semantic_summary_path.name,
|
"semantic_footprint_count": semantic_footprints,
|
"semantic_class_counts": {details["key"]: int((semantic_classes == code).sum()) for code, details in SEMANTIC_CLASSES.items()},
|
"semantic_source_has_rgb": source_data.source_has_rgb,
|
"semantic_method": "CPU RGB and local-geometry rules; narrow line and pole/tower candidates require manual review",
|
"mesh_file": mesh.name if triangles else None,
|
"dsm_file": dsm_path.name,
|
"height_file": height_path.name,
|
"label_raster_file": label_path.name,
|
"preview_file": preview_path.name,
|
"vector_file": vector_path.name,
|
"summary_file": summary_path.name,
|
"vectorizer": vectorizer,
|
"ground_plane": [round(float(value), 8) for value in plane],
|
"ground_plane_available": ground_plane_available,
|
"coordinate_alignment": alignment,
|
"elapsed_seconds": round(time.perf_counter() - started, 3),
|
}
|
|
|
def collect_inputs(input_path: Path) -> list[Path]:
|
candidates = [input_path] if input_path.is_file() else sorted(path for path in input_path.iterdir() if path.is_file()) if input_path.is_dir() else []
|
inputs = [path for path in candidates if path.suffix.lower() in SUPPORTED_SUFFIXES]
|
if not inputs:
|
raise ValueError("No supported PLY/PCD/XYZ/LAS/LAZ point-cloud inputs were found.")
|
return inputs
|
|
|
def main() -> int:
|
parser = argparse.ArgumentParser(description="Understand an existing point-cloud reconstruction on CPU.")
|
parser.add_argument("--input", type=Path, required=True)
|
parser.add_argument("--output", type=Path, required=True)
|
parser.add_argument("--voxel-size", type=float, default=0.2)
|
parser.add_argument("--ground-distance", type=float, default=0.15)
|
parser.add_argument("--elevated-threshold", type=float, default=0.75)
|
parser.add_argument("--ground-up-axis", choices={"z", "auto"}, default="z", help="Use source Z as up, or rotate a processing copy so its dominant plane is horizontal.")
|
args = parser.parse_args()
|
if args.output.exists() and any(args.output.iterdir()):
|
raise SystemExit("Output directory is not empty; use a new run directory.")
|
args.output.mkdir(parents=True, exist_ok=True)
|
started = time.perf_counter()
|
images = [process_point_cloud(path, args.output, args.voxel_size, args.ground_distance, args.elevated_threshold, args.ground_up_axis) for path in collect_inputs(args.input)]
|
metadata = {
|
"capability": "05-3d-pointcloud",
|
"classification": "B",
|
"created_at": datetime.now(UTC).isoformat(),
|
"versions": {"geoai-py": importlib_metadata.version("geoai-py"), "open3d": o3d.__version__, "laspy": importlib_metadata.version("laspy")},
|
"method": "Open3D voxel downsampling, local-ground/RGB geometry rules, alpha-shape mesh reconstruction, GeoAI elevated-footprint vectorization, and Rasterio semantic-footprint vectorization",
|
"model": "none (explainable CPU geometry/RGB classification baseline)",
|
"device": "CPU",
|
"thresholds": {"voxel_size_local_units": args.voxel_size, "ground_plane_distance": args.ground_distance, "elevated_height": args.elevated_threshold, "ground_up_axis": args.ground_up_axis},
|
"input_count": len(images),
|
"processed_point_clouds": len(images),
|
"elapsed_seconds": round(time.perf_counter() - started, 3),
|
"point_clouds": images,
|
"limitations": [
|
"This Demo consumes an existing point cloud; it does not reconstruct a point cloud from photographs or provide photogrammetric camera calibration.",
|
"Semantic labels are explainable CPU rules, not a trained point-cloud model. Building/structure, vegetation, power-line and pole/tower outputs are review candidates, not verified asset inventory or inspection conclusions.",
|
"Power-line and pole/tower rules require thin/linear or tall/vertical local evidence. Wires hidden by vegetation, bundled conductors, tree trunks, roof edges and lattice structures can be missed or falsely labelled.",
|
"Inputs without declared CRS remain in local point-cloud coordinates. The GeoTIFF and GeoJSON do not represent latitude/longitude or surveyed map coordinates.",
|
"RANSAC assumes a dominant approximately horizontal ground plane; slopes, cliffs, dense vegetation, water, or large vertical structures can cause misses and false positives.",
|
"Mesh triangles are an inspectable alpha-shape approximation, not a watertight or accuracy-validated reconstruction.",
|
] + (["At least one source had no plausible dominant horizontal RANSAC plane. Its local low-point surface is a fallback only, so ground/elevated and semantic classes need heightened manual review."] if any(not item["ground_plane_available"] for item in images) else []),
|
}
|
(args.output / "run_metadata.json").write_text(json.dumps(metadata, ensure_ascii=False, indent=2), encoding="utf-8")
|
print(json.dumps(metadata, ensure_ascii=False, indent=2))
|
return 0
|
|
|
if __name__ == "__main__":
|
raise SystemExit(main())
|
