520 lines
14 KiB
Python
520 lines
14 KiB
Python
"""Enhanced building export with point cloud roof extraction.
|
|
|
|
Combines CityGML footprints with BDOM20RGBI/LPO point clouds for:
|
|
1. RANSAC-based roof plane extraction (where BDOM available)
|
|
2. LPO height refinement
|
|
3. DOP20 facade texturing
|
|
4. Fallback to standard CityGML pipeline
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import csv
|
|
import glob
|
|
import os
|
|
from dataclasses import dataclass
|
|
from typing import List, Optional, Tuple
|
|
|
|
import numpy as np
|
|
from osgeo import gdal
|
|
|
|
from .config import Config
|
|
from .gdal_utils import open_dataset
|
|
from .pointcloud import (
|
|
PointCloud,
|
|
find_laz_file,
|
|
find_pointcloud_file,
|
|
has_bdom_data,
|
|
has_lpo_data,
|
|
read_laz_file,
|
|
read_pointcloud_file,
|
|
sample_rgb_average,
|
|
)
|
|
|
|
try:
|
|
from shapely.geometry import Polygon
|
|
from shapely.ops import unary_union
|
|
HAS_SHAPELY = True
|
|
except ImportError:
|
|
HAS_SHAPELY = False
|
|
|
|
|
|
@dataclass
|
|
class RoofPlane:
|
|
"""A detected roof plane from RANSAC."""
|
|
normal: np.ndarray # Unit normal vector
|
|
d: float # Plane distance from origin
|
|
inliers: np.ndarray # Point indices
|
|
centroid: np.ndarray # Plane centroid
|
|
|
|
|
|
def ransac_plane_fit(
|
|
points: np.ndarray,
|
|
max_iterations: int = 1000,
|
|
distance_threshold: float = 0.15,
|
|
min_inliers: int = 50,
|
|
) -> Optional[RoofPlane]:
|
|
"""Fit a single plane using RANSAC.
|
|
|
|
Args:
|
|
points: Nx3 array of points
|
|
max_iterations: Maximum RANSAC iterations
|
|
distance_threshold: Inlier distance threshold in meters
|
|
min_inliers: Minimum inliers for valid plane
|
|
|
|
Returns:
|
|
RoofPlane or None if no plane found
|
|
"""
|
|
if len(points) < 3:
|
|
return None
|
|
|
|
best_inliers = None
|
|
best_normal = None
|
|
best_d = None
|
|
|
|
n_points = len(points)
|
|
|
|
for _ in range(max_iterations):
|
|
# Sample 3 random points
|
|
idx = np.random.choice(n_points, 3, replace=False)
|
|
p1, p2, p3 = points[idx]
|
|
|
|
# Compute plane normal
|
|
v1 = p2 - p1
|
|
v2 = p3 - p1
|
|
normal = np.cross(v1, v2)
|
|
|
|
norm_len = np.linalg.norm(normal)
|
|
if norm_len < 1e-10:
|
|
continue
|
|
|
|
normal = normal / norm_len
|
|
|
|
# Plane equation: normal . (x - p1) = 0
|
|
# => normal . x = normal . p1 = d
|
|
d = np.dot(normal, p1)
|
|
|
|
# Compute distances to plane
|
|
distances = np.abs(np.dot(points, normal) - d)
|
|
|
|
# Find inliers
|
|
inlier_mask = distances < distance_threshold
|
|
inlier_count = np.sum(inlier_mask)
|
|
|
|
if inlier_count >= min_inliers:
|
|
if best_inliers is None or inlier_count > len(best_inliers):
|
|
best_inliers = np.where(inlier_mask)[0]
|
|
best_normal = normal
|
|
best_d = d
|
|
|
|
if best_inliers is None:
|
|
return None
|
|
|
|
# Refit plane using all inliers (least squares)
|
|
inlier_points = points[best_inliers]
|
|
centroid = np.mean(inlier_points, axis=0)
|
|
|
|
# SVD to find best-fit plane
|
|
centered = inlier_points - centroid
|
|
_, _, vh = np.linalg.svd(centered)
|
|
refined_normal = vh[2] # Smallest singular value direction
|
|
|
|
# Ensure normal points upward
|
|
if refined_normal[2] < 0:
|
|
refined_normal = -refined_normal
|
|
|
|
refined_d = np.dot(refined_normal, centroid)
|
|
|
|
return RoofPlane(
|
|
normal=refined_normal,
|
|
d=refined_d,
|
|
inliers=best_inliers,
|
|
centroid=centroid,
|
|
)
|
|
|
|
|
|
def extract_roof_planes(
|
|
points: PointCloud,
|
|
max_planes: int = 6,
|
|
ransac_threshold: float = 0.15,
|
|
ransac_iterations: int = 1000,
|
|
min_inliers: int = 50,
|
|
) -> List[RoofPlane]:
|
|
"""Extract multiple roof planes using iterative RANSAC.
|
|
|
|
Args:
|
|
points: Point cloud of roof points
|
|
max_planes: Maximum number of planes to extract
|
|
ransac_threshold: Distance threshold for inliers
|
|
ransac_iterations: RANSAC iterations per plane
|
|
min_inliers: Minimum points for valid plane
|
|
|
|
Returns:
|
|
List of extracted RoofPlane objects
|
|
"""
|
|
if len(points) < min_inliers:
|
|
return []
|
|
|
|
xyz = points.xyz.copy()
|
|
planes = []
|
|
remaining_mask = np.ones(len(xyz), dtype=bool)
|
|
|
|
for _ in range(max_planes):
|
|
remaining_points = xyz[remaining_mask]
|
|
|
|
if len(remaining_points) < min_inliers:
|
|
break
|
|
|
|
plane = ransac_plane_fit(
|
|
remaining_points,
|
|
max_iterations=ransac_iterations,
|
|
distance_threshold=ransac_threshold,
|
|
min_inliers=min_inliers,
|
|
)
|
|
|
|
if plane is None:
|
|
break
|
|
|
|
# Map inlier indices back to original array
|
|
remaining_indices = np.where(remaining_mask)[0]
|
|
original_inliers = remaining_indices[plane.inliers]
|
|
plane.inliers = original_inliers
|
|
|
|
planes.append(plane)
|
|
|
|
# Remove inliers from remaining points
|
|
remaining_mask[original_inliers] = False
|
|
|
|
return planes
|
|
|
|
|
|
def classify_roof_type(planes: List[RoofPlane]) -> str:
|
|
"""Classify roof type from extracted planes.
|
|
|
|
Args:
|
|
planes: List of roof planes
|
|
|
|
Returns:
|
|
Roof type string: 'flat', 'gable', 'hip', 'complex'
|
|
"""
|
|
if not planes:
|
|
return "unknown"
|
|
|
|
n_planes = len(planes)
|
|
|
|
if n_planes == 1:
|
|
# Check if plane is horizontal (flat roof)
|
|
normal = planes[0].normal
|
|
if abs(normal[2]) > 0.95: # Nearly vertical normal
|
|
return "flat"
|
|
return "shed"
|
|
|
|
if n_planes == 2:
|
|
# Check for gable roof (two planes meeting at ridge)
|
|
n1, n2 = planes[0].normal, planes[1].normal
|
|
dot = abs(np.dot(n1, n2))
|
|
if dot < 0.9: # Planes at an angle
|
|
return "gable"
|
|
|
|
if n_planes == 4:
|
|
return "hip"
|
|
|
|
return "complex"
|
|
|
|
|
|
def refine_building_with_lpo(
|
|
tile_id: str,
|
|
footprint_bounds: Tuple[float, float, float, float],
|
|
citygml_height: float,
|
|
cfg: Config,
|
|
) -> Optional[float]:
|
|
"""Refine building height using LPO points.
|
|
|
|
Args:
|
|
tile_id: Tile identifier
|
|
footprint_bounds: (xmin, ymin, xmax, ymax) of building
|
|
citygml_height: Original height from CityGML
|
|
cfg: Configuration
|
|
|
|
Returns:
|
|
Refined height or None if no LPO data
|
|
"""
|
|
lpo_file = find_pointcloud_file(cfg.pointcloud.lpo_dir, tile_id)
|
|
if not lpo_file:
|
|
return None
|
|
|
|
xmin, ymin, xmax, ymax = footprint_bounds
|
|
|
|
# Load LPO points in building bounds
|
|
lpo = read_pointcloud_file(lpo_file, bounds=(xmin, ymin, xmax, ymax))
|
|
if len(lpo) < 10:
|
|
return None
|
|
|
|
# Get 95th percentile height
|
|
max_z = float(np.percentile(lpo.z, 95))
|
|
|
|
# Load DGM for ground elevation
|
|
dgm_pattern = os.path.join(cfg.raw.dgm1_dir, f"*{tile_id}*.tif")
|
|
dgm_files = glob.glob(dgm_pattern)
|
|
if not dgm_files:
|
|
return None
|
|
|
|
dgm_ds = open_dataset(dgm_files[0], required=False)
|
|
if dgm_ds is None:
|
|
return None
|
|
|
|
dgm = dgm_ds.GetRasterBand(1).ReadAsArray()
|
|
gt = dgm_ds.GetGeoTransform()
|
|
|
|
# Sample ground elevation at building center
|
|
center_x = (xmin + xmax) / 2
|
|
center_y = (ymin + ymax) / 2
|
|
col = int((center_x - gt[0]) / gt[1])
|
|
row = int((center_y - gt[3]) / gt[5])
|
|
|
|
if 0 <= row < dgm.shape[0] and 0 <= col < dgm.shape[1]:
|
|
ground_z = float(dgm[row, col])
|
|
refined_height = max_z - ground_z
|
|
if refined_height > 0:
|
|
return refined_height
|
|
|
|
return None
|
|
|
|
|
|
def sample_facade_color(
|
|
tile_id: str,
|
|
world_x: float,
|
|
world_y: float,
|
|
cfg: Config,
|
|
) -> Tuple[int, int, int]:
|
|
"""Sample facade color from orthophoto.
|
|
|
|
Args:
|
|
tile_id: Tile identifier
|
|
world_x, world_y: World coordinates
|
|
cfg: Configuration
|
|
|
|
Returns:
|
|
RGB tuple (0-255)
|
|
"""
|
|
ortho_path = os.path.join(cfg.export.ortho_dir, f"{tile_id}.jpg")
|
|
if not os.path.exists(ortho_path):
|
|
return (180, 170, 160) # Default neutral
|
|
|
|
ds = open_dataset(ortho_path, required=False)
|
|
if ds is None:
|
|
return (180, 170, 160)
|
|
|
|
gt = ds.GetGeoTransform()
|
|
col = int((world_x - gt[0]) / gt[1])
|
|
row = int((world_y - gt[3]) / gt[5])
|
|
|
|
if 0 <= row < ds.RasterYSize and 0 <= col < ds.RasterXSize:
|
|
try:
|
|
r = ds.GetRasterBand(1).ReadAsArray(col, row, 1, 1)[0, 0]
|
|
g = ds.GetRasterBand(2).ReadAsArray(col, row, 1, 1)[0, 0]
|
|
b = ds.GetRasterBand(3).ReadAsArray(col, row, 1, 1)[0, 0]
|
|
return (int(r), int(g), int(b))
|
|
except Exception:
|
|
pass
|
|
|
|
return (180, 170, 160)
|
|
|
|
|
|
def process_building_with_bdom(
|
|
tile_id: str,
|
|
footprint_bounds: Tuple[float, float, float, float],
|
|
cfg: Config,
|
|
) -> Optional[dict]:
|
|
"""Process a building using BDOM point cloud for roof geometry.
|
|
|
|
Args:
|
|
tile_id: Tile identifier
|
|
footprint_bounds: (xmin, ymin, xmax, ymax) of building
|
|
cfg: Configuration
|
|
|
|
Returns:
|
|
Dict with roof planes and metadata, or None if processing fails
|
|
"""
|
|
eb_cfg = cfg.buildings_enhanced
|
|
|
|
bdom_file = find_laz_file(cfg.pointcloud.bdom_dir, tile_id)
|
|
if not bdom_file:
|
|
return None
|
|
|
|
xmin, ymin, xmax, ymax = footprint_bounds
|
|
|
|
# Load BDOM points in building bounds (with small buffer)
|
|
buffer = 0.5
|
|
bdom = read_laz_file(bdom_file, bounds=(
|
|
xmin - buffer, ymin - buffer,
|
|
xmax + buffer, ymax + buffer
|
|
))
|
|
|
|
if len(bdom) < eb_cfg.min_plane_inliers:
|
|
return None
|
|
|
|
# Load DGM for ground reference
|
|
dgm_pattern = os.path.join(cfg.raw.dgm1_dir, f"*{tile_id}*.tif")
|
|
dgm_files = glob.glob(dgm_pattern)
|
|
if not dgm_files:
|
|
return None
|
|
|
|
dgm_ds = open_dataset(dgm_files[0], required=False)
|
|
if dgm_ds is None:
|
|
return None
|
|
|
|
dgm = dgm_ds.GetRasterBand(1).ReadAsArray()
|
|
gt = dgm_ds.GetGeoTransform()
|
|
|
|
# Get ground elevation at building center
|
|
center_x = (xmin + xmax) / 2
|
|
center_y = (ymin + ymax) / 2
|
|
col = int((center_x - gt[0]) / gt[1])
|
|
row = int((center_y - gt[3]) / gt[5])
|
|
|
|
if 0 <= row < dgm.shape[0] and 0 <= col < dgm.shape[1]:
|
|
ground_z = float(dgm[row, col])
|
|
else:
|
|
ground_z = float(np.min(bdom.z))
|
|
|
|
# Filter to points above ground (roof points)
|
|
roof_threshold = ground_z + 2.0 # At least 2m above ground
|
|
roof_mask = bdom.z > roof_threshold
|
|
roof_points = PointCloud(
|
|
x=bdom.x[roof_mask],
|
|
y=bdom.y[roof_mask],
|
|
z=bdom.z[roof_mask],
|
|
r=bdom.r[roof_mask] if bdom.r is not None else None,
|
|
g=bdom.g[roof_mask] if bdom.g is not None else None,
|
|
b=bdom.b[roof_mask] if bdom.b is not None else None,
|
|
)
|
|
|
|
if len(roof_points) < eb_cfg.min_plane_inliers:
|
|
return None
|
|
|
|
# Extract roof planes
|
|
planes = extract_roof_planes(
|
|
roof_points,
|
|
max_planes=eb_cfg.max_roof_planes,
|
|
ransac_threshold=eb_cfg.ransac_threshold_m,
|
|
ransac_iterations=eb_cfg.ransac_iterations,
|
|
min_inliers=eb_cfg.min_plane_inliers,
|
|
)
|
|
|
|
if not planes:
|
|
return None
|
|
|
|
# Get roof RGB if available
|
|
roof_rgb = None
|
|
if roof_points.has_rgb:
|
|
r_avg = int(np.mean(roof_points.r))
|
|
g_avg = int(np.mean(roof_points.g))
|
|
b_avg = int(np.mean(roof_points.b))
|
|
roof_rgb = (r_avg, g_avg, b_avg)
|
|
|
|
# Compute building height
|
|
max_z = float(np.max(roof_points.z))
|
|
height = max_z - ground_z
|
|
|
|
return {
|
|
"planes": planes,
|
|
"roof_type": classify_roof_type(planes),
|
|
"height": height,
|
|
"ground_z": ground_z,
|
|
"roof_rgb": roof_rgb,
|
|
"source": "bdom",
|
|
}
|
|
|
|
|
|
def export_buildings_enhanced(cfg: Config) -> int:
|
|
"""Export enhanced buildings with point cloud fusion.
|
|
|
|
For tiles with BDOM20RGBI:
|
|
- Extract roof planes using RANSAC
|
|
- Use RGB from point cloud for facade colors
|
|
|
|
For tiles with only LPO:
|
|
- Refine building heights from LPO
|
|
- Use DOP20 for facade colors
|
|
|
|
For tiles without point cloud:
|
|
- Fall back to standard CityGML pipeline
|
|
|
|
Args:
|
|
cfg: Configuration
|
|
|
|
Returns:
|
|
Exit code (0 = success)
|
|
"""
|
|
eb_cfg = cfg.buildings_enhanced
|
|
|
|
# Check for manifest
|
|
if not os.path.exists(cfg.export.manifest_path):
|
|
print(f"[buildings_enhanced] ERROR: Manifest not found: {cfg.export.manifest_path}")
|
|
print("[buildings_enhanced] Run heightmap export first.")
|
|
return 1
|
|
|
|
# Load tiles
|
|
tiles = []
|
|
with open(cfg.export.manifest_path, "r") as f:
|
|
reader = csv.DictReader(f)
|
|
for row in reader:
|
|
tiles.append({
|
|
"tile_id": row["tile_id"],
|
|
"xmin": float(row["xmin"]),
|
|
"ymin": float(row["ymin"]),
|
|
"xmax": float(row["xmax"]),
|
|
"ymax": float(row["ymax"]),
|
|
})
|
|
|
|
if not tiles:
|
|
print("[buildings_enhanced] No tiles in manifest.")
|
|
return 1
|
|
|
|
os.makedirs(eb_cfg.out_dir, exist_ok=True)
|
|
|
|
stats = {"bdom": 0, "lpo": 0, "citygml": 0, "failed": 0}
|
|
|
|
for tile in tiles:
|
|
tile_id = tile["tile_id"]
|
|
print(f"[buildings_enhanced] Processing {tile_id}...")
|
|
|
|
# Check data availability
|
|
has_bdom = has_bdom_data(tile_id, cfg.pointcloud.bdom_dir)
|
|
has_lpo = has_lpo_data(tile_id, cfg.pointcloud.lpo_dir)
|
|
|
|
if has_bdom and eb_cfg.use_bdom_roof:
|
|
print(f"[buildings_enhanced] {tile_id}: Using BDOM for roof extraction")
|
|
stats["bdom"] += 1
|
|
# BDOM processing would create enhanced GLB here
|
|
# For now, we log that BDOM is available
|
|
|
|
elif has_lpo and eb_cfg.use_lpo_refinement:
|
|
print(f"[buildings_enhanced] {tile_id}: Using LPO for height refinement")
|
|
stats["lpo"] += 1
|
|
# LPO refinement would adjust CityGML heights here
|
|
|
|
elif eb_cfg.fallback_to_citygml:
|
|
print(f"[buildings_enhanced] {tile_id}: Falling back to CityGML")
|
|
stats["citygml"] += 1
|
|
# Standard CityGML export
|
|
|
|
else:
|
|
print(f"[buildings_enhanced] {tile_id}: No data available, skipping")
|
|
stats["failed"] += 1
|
|
|
|
print(f"[buildings_enhanced] DONE. "
|
|
f"BDOM={stats['bdom']}, LPO={stats['lpo']}, "
|
|
f"CityGML={stats['citygml']}, Failed={stats['failed']}")
|
|
|
|
# Note: Full GLB export would integrate with existing buildings.py
|
|
# This module provides the data fusion logic; actual mesh export
|
|
# requires integration with trimesh and GLB generation from buildings.py
|
|
|
|
if stats["bdom"] > 0 or stats["lpo"] > 0:
|
|
print("[buildings_enhanced] Note: Full GLB export requires integration with buildings.py")
|
|
print("[buildings_enhanced] Use --export buildings for standard CityGML export")
|
|
|
|
return 0
|