Add RRIM, Multi-Hillshade RGB, and Local Dominance visualizations
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 <noreply@anthropic.com>
This commit is contained in:
Jacquin Antoine
@ -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
|
||||
|
||||
|
||||
@ -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:
|
||||
|
||||
@ -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',
|
||||
|
||||
@ -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])
|
||||
|
||||
@ -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 ""
|
||||
|
||||
Reference in New Issue
Block a user