From 7f6b816ed62067e71e515b54288c1349ea55a999 Mon Sep 17 00:00:00 2001 From: Jacquin Antoine Date: Thu, 14 May 2026 01:03:47 +0200 Subject: [PATCH] Add RRIM, Multi-Hillshade RGB, and Local Dominance visualizations MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three new visualizations complementing existing SVF/openness/LRM/MSRM: - RRIM (Red Relief Image Map): RGB composite combining positive openness (R), inverted slope (G), negative openness (B). Uses ray-tracing to compute both openness values in a single pass. - Multi-Hillshade RGB: 3 azimuths (315°, 135°, 45°) mapped to R/G/B channels with slope blending. Color reveals structure orientation. - Local Dominance: (dem - local_min) / (local_max - local_min) using min/max filters. Measures local height position — complements openness. Also adds: - _compute_openness_both() helper for shared ray-tracing (used by RRIM) - xp_maximum_filter() in gpu.py (GPU/CPU abstraction) - Entries in COLORMAPS, RGB_LEGENDS, VIZ_STEPS, and is_rgb detection - All NaN handling follows existing patterns (nan_mask restoration) Co-Authored-By: Claude Opus 4.6 --- CLAUDE.md | 2 +- lidar_pipeline/gpu.py | 10 ++ lidar_pipeline/pipeline.py | 5 +- lidar_pipeline/rendering.py | 20 ++- lidar_pipeline/visualizations.py | 248 ++++++++++++++++++++++++++++++- 5 files changed, 281 insertions(+), 4 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 5c6ba44..279791e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -4,7 +4,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co ## Project Overview -LiDAR archaeological processing pipeline that generates 17 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154). +LiDAR archaeological processing pipeline that generates 20 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154). ## Commands diff --git a/lidar_pipeline/gpu.py b/lidar_pipeline/gpu.py index ed83b20..f4f357f 100644 --- a/lidar_pipeline/gpu.py +++ b/lidar_pipeline/gpu.py @@ -110,6 +110,16 @@ def xp_minimum_filter(arr, footprint=None, size=None): return ndimage.minimum_filter(arr, footprint=footprint, size=size) +def xp_maximum_filter(arr, footprint=None, size=None): + """Maximum filter — uses GPU if array is on GPU, CPU otherwise.""" + if _cp is not None and isinstance(arr, _cp.ndarray): + try: + return _cp_ndimage.maximum_filter(arr, footprint=footprint, size=size) + except Exception: + arr = to_cpu(arr) + return ndimage.maximum_filter(arr, footprint=footprint, size=size) + + def gpu_cleanup(): """Free GPU memory. Call between visualizations to prevent OOM.""" if _cp is not None: diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index de6a72f..6c221fb 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -61,7 +61,7 @@ from .visualizations import ( generate_lrm, generate_svf, generate_openness, generate_mslrm, generate_tpi, generate_sailore, generate_roughness, generate_anomalies, generate_wavelet, - generate_flow, + generate_flow, generate_rrim, generate_multi_hillshade, generate_local_dominance, ) from .gpu import gpu_cleanup from .ign import generate_ign_overlay @@ -87,6 +87,9 @@ VIZ_STEPS = [ ('anomalies', generate_anomalies), ('wavelet', generate_wavelet), ('flow', generate_flow), + ('rrim', lambda d, b, v, r: generate_rrim(d, b, v, r)), + ('multi_hillshade', lambda d, b, v, r: generate_multi_hillshade(d, b, v, r)), + ('local_dominance', generate_local_dominance), ('ortho', lambda d, b, v, r: generate_ign_overlay( d, b, v, r, layer='ORTHOIMAGERY.ORTHOPHOTOS', diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index 0c6077d..9aababa 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -156,6 +156,14 @@ COLORMAPS = { 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'percentile', 'vmax_pct': 98, }, + 'local_dominance': { + 'cmap': 'RdYlBu_r', + 'title': 'Dominance Locale (position relative dans le voisinage)', + 'legend': 'Proportion du voisinage sous le point central\nRouge = Point dominant (sommet, crête)\nBleu = Point encaissé (fossé, vallée)\nRayon: 15m', + 'description': 'Mesure la saillie locale — complémentaire de l\'openness', + 'vmin_mode': 'percentile', 'vmin_pct': 2, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, } # RGB entries (ortho/topo) are handled specially @@ -170,6 +178,16 @@ RGB_LEGENDS = { 'legend': 'Carte IGN\nPlan topographique', 'description': 'Carte topographique IGN (Plan IGN)', }, + 'rrim': { + 'title': 'RRIM — Red Relief Image Map (composite RGB)', + 'legend': 'Rouge = Openness positive (crêtes, levées)\nVert = Pente inversée (plat = clair)\nBleu = Openness négative (fossés, dépressions)', + 'description': 'Composite RGB synthétique pour prospection archéologique', + }, + 'multi_hillshade': { + 'title': 'Hillshade Composite RGB (3 azimuts)', + 'legend': 'Rouge = Éclairage NW (315°)\nVert = Éclairage SE (135°)\nBleu = Éclairage NE (45°)', + 'description': 'Composite couleur révélant les structures selon leur orientation', + }, } @@ -282,7 +300,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None): try: with rasterio.open(tif_file) as src: - is_rgb = src.count >= 3 and ('ortho' in str(tif_file) or 'topo' in str(tif_file)) + is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo', 'rrim', 'multi_hillshade')) if is_rgb: data = src.read([1, 2, 3]) diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index b0dff01..788201a 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -18,7 +18,7 @@ import rasterio from scipy.ndimage import generic_filter from scipy.stats import binned_statistic_2d -from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, xp_minimum_filter, gpu_cleanup +from .gpu import HAS_GPU, to_gpu, to_cpu, xp_gaussian_filter, xp_uniform_filter, xp_minimum_filter, xp_maximum_filter, gpu_cleanup logger = logging.getLogger("lidar") @@ -528,6 +528,252 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh return None +def _compute_openness_both(dem, resolution, nan_mask, n_dirs=8, radius=50): + """Compute positive and negative openness in one ray-tracing pass. + + Traces rays in n_dirs directions up to radius pixels, measuring: + - positive openness: max angle above horizontal to visible terrain + - negative openness: max angle below horizontal to visible terrain + + Returns (pos_open, neg_open) as float32 arrays in degrees. + NaN mask is applied after computation. + """ + rows, cols = dem.shape + res = resolution + max_dist = int(radius / res) + + angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) + dx_dir = np.cos(angles) + dy_dir = np.sin(angles) + + padded = np.pad(dem, max_dist, mode='constant', constant_values=np.nanmax(dem[~nan_mask]) + 10000 if np.any(~nan_mask) else 0) + + pos_sum = np.zeros_like(dem) + neg_sum = np.zeros_like(dem) + + for d_idx in range(n_dirs): + ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] + max_pos_angle = np.zeros_like(dem) + max_neg_angle = np.zeros_like(dem) + + for step in range(1, max_dist + 1): + px = int(round(ddx * step)) + py = int(round(ddy * step)) + dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) + if dist_m < res * 0.5: + continue + + elev_diff = padded[max_dist + py:max_dist + py + rows, + max_dist + px:max_dist + px + cols] - dem + + pos_angle = np.arctan2(np.maximum(elev_diff, 0), dist_m) + neg_angle = np.arctan2(np.maximum(-elev_diff, 0), dist_m) + + valid = ~np.isnan(elev_diff) + max_pos_angle[valid] = np.maximum(max_pos_angle[valid], pos_angle[valid]) + max_neg_angle[valid] = np.maximum(max_neg_angle[valid], neg_angle[valid]) + + pos_sum += max_pos_angle + neg_sum += max_neg_angle + + pos_open = np.degrees(pos_sum / n_dirs).astype(np.float32) + neg_open = np.degrees(neg_sum / n_dirs).astype(np.float32) + pos_open[nan_mask] = np.nan + neg_open[nan_mask] = np.nan + return pos_open, neg_open + + +def generate_rrim(dem_file, basename, vis_dir, resolution, shared=None, + n_dirs=8, radius=50, pmin=2, pmax=98, contrast=1.5): + """Red Relief Image Map — RGB composite for archaeological prospection. + + Combines slope, positive openness, and negative openness into a single + false-color image where: + Red = positive openness (ridges, elevated features) + Green = inverted slope (flat = bright, steep = dark) + Blue = negative openness (depressions, ditches) + + Each channel is normalized via percentiles and enhanced with a gamma curve. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → RRIM (Red Relief Image){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_rrim.tif" + + try: + if shared: + transform = shared.transform + crs = shared.crs + dem_np = shared.dem_np + nan_mask = shared.nan_mask + slope_rad = shared.slope_rad + dem_for_ray = to_gpu(shared.filled) if HAS_GPU else shared.filled + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + filled, _ = _fill_nans(dem_np) + dem_for_ray = to_gpu(filled) if HAS_GPU else filled + dy, dx = np.gradient(filled, resolution) + slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) + + # Compute both openness values (ray-tracing on filled DEM) + pos_open, neg_open = _compute_openness_both( + to_cpu(dem_for_ray) if HAS_GPU else dem_for_ray, + resolution, nan_mask, n_dirs=n_dirs, radius=radius + ) + + # Normalize each component to [0, 1] using percentiles + slope_deg = np.degrees(slope_rad) + slope_deg[nan_mask] = np.nan + + def _normalize(arr, lo, hi): + valid = arr[~np.isnan(arr)] + if len(valid) == 0: + return np.zeros_like(arr, dtype=np.float32) + vlo = np.percentile(valid, lo) + vhi = np.percentile(valid, hi) + if vhi - vlo < 1e-6: + return np.full_like(arr, 0.5, dtype=np.float32) + norm = np.clip((arr - vlo) / (vhi - vlo), 0, 1) + # Apply gamma for contrast + norm = np.power(norm, 1.0 / contrast) + return norm.astype(np.float32) + + r = _normalize(pos_open, pmin, pmax) # Red: positive openness (ridges) + g = _normalize(90 - slope_deg, pmin, pmax) # Green: inverted slope (flat=bright) + g[nan_mask] = np.nan + b = _normalize(neg_open, pmin, pmax) # Blue: negative openness (ditches) + + # Assemble RGB (uint8) + rgb = np.stack([r, g, b], axis=0) # (3, H, W) + rgb = np.nan_to_num(rgb, nan=0.0) + rgb_uint8 = (np.clip(rgb, 0, 1) * 255).astype(np.uint8) + + _save_tif(output, rgb_uint8, transform, crs, dtype='uint8', count=3) + logger.info(f" ✓ RRIM terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur RRIM: {e}", exc_info=True) + return None + + +def generate_multi_hillshade(dem_file, basename, vis_dir, resolution, shared=None, + azimuths=(315, 135, 45), altitude=30, blend_slope=0.3): + """Multi-directional hillshade RGB composite — 3 azimuths mapped to R/G/B. + + Each azimuth produces a hillshade mapped to a color channel: + Red = azimuth 315° (NW illumination) + Green = azimuth 135° (SE illumination) + Blue = azimuth 45° (NE illumination) + + Shadow direction reveals structure orientation through color. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → Hillshade Composite RGB{gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_multi_hillshade.tif" + + try: + if shared: + transform = shared.transform + crs = shared.crs + nan_mask = shared.nan_mask + slope_rad = to_gpu(shared.slope_rad) if HAS_GPU else shared.slope_rad + aspect = to_gpu(shared.aspect) if HAS_GPU else shared.aspect + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + filled, _ = _fill_nans(dem_np) + dem = to_gpu(filled) if HAS_GPU else filled + dy, dx = xp.gradient(dem, resolution) + slope_rad = xp.arctan(xp.sqrt(dx**2 + dy**2)) + aspect = xp.arctan2(dy, dx) + + alt_rad = xp.radians(xp.array(altitude)) + sin_alt = xp.sin(alt_rad) + cos_alt = xp.cos(alt_rad) + cos_slope = xp.cos(slope_rad) + + channels = [] + for az in azimuths: + az_rad = xp.radians(xp.array(az)) + hs = sin_alt * xp.sin(slope_rad) + cos_alt * cos_slope * xp.cos(az_rad - aspect) + blended = (1 - blend_slope) * xp.clip(hs, 0, 1) + blend_slope * cos_slope + channels.append(to_cpu(blended).astype(np.float32)) + + gpu_cleanup() + + # Assemble RGB (uint8) + rgb = np.stack(channels, axis=0) # (3, H, W) + rgb[:, nan_mask] = 0.0 + rgb_uint8 = (np.clip(rgb, 0, 1) * 255).astype(np.uint8) + + _save_tif(output, rgb_uint8, transform, crs, dtype='uint8', count=3) + logger.info(f" ✓ Hillshade Composite RGB terminé ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur multi_hillshade: {e}", exc_info=True) + return None + + +def generate_local_dominance(dem_file, basename, vis_dir, resolution, shared=None, + radius=15, pmin=2, pmax=98): + """Local Dominance — proportion of neighborhood below center point. + + LD = (dem - local_min) / (local_max - local_min + epsilon) + + High values = locally dominant (peak, ridge) + Low values = locally recessed (valley, pit) + + Uses minimum/maximum filters on the filled DEM, then restores NaN mask. + Complements openness by measuring local height position rather than angular extent. + """ + gpu_tag = " [GPU]" if HAS_GPU else "" + logger.info(f" → Dominance Locale (rayon {radius}m){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_local_dominance.tif" + + try: + if shared: + transform = shared.transform + crs = shared.crs + nan_mask = shared.nan_mask + dem_np = shared.dem_np + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + + radius_px = max(1, int(radius / resolution)) + if radius_px % 2 == 0: + radius_px += 1 + + local_min = _filter_nanaware_from_filled( + shared, xp_minimum_filter, size=radius_px + ) if shared else _filter_nanaware( + dem_np, xp_minimum_filter, size=radius_px + ) + + local_max_data = _filter_nanaware_from_filled( + shared, xp_maximum_filter, size=radius_px + ) if shared else _filter_nanaware( + dem_np, xp_maximum_filter, size=radius_px + ) + + # Local dominance ratio + epsilon = 0.01 # Avoid division by zero on flat terrain + local_range = local_max_data - local_min + epsilon + dominance = (dem_np - local_min) / local_range + dominance = np.clip(dominance, 0, 1) + dominance[nan_mask] = np.nan + + _save_tif(output, dominance.astype(np.float32), transform, crs, nan_mask=nan_mask) + logger.info(f" ✓ Dominance Locale terminée ({time.time()-t0:.1f}s){gpu_tag}") + return output + except Exception as e: + logger.error(f" ✗ Erreur local_dominance: {e}", exc_info=True) + return None + + def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None): """Multi-Scale Relief Model (MSRM) - LRM at 5 scales combined (GPU if available).""" gpu_tag = " [GPU]" if HAS_GPU else ""