From 929fac9aa061053230d5036ae3fda0a18acf21e9 Mon Sep 17 00:00:00 2001 From: Antoine Jacquin Date: Sun, 31 May 2026 17:19:31 +0200 Subject: [PATCH] Remove LRM, TPI, aspect, curvature, paths + add flow accumulation, directional Gabor wavelets, multi-radius ray-tracing --- lidar_pipeline/pipeline.py | 21 +- lidar_pipeline/rendering.py | 56 +- lidar_pipeline/tests/test_visualizations.py | 117 +-- lidar_pipeline/visualizations.py | 797 +++++++++----------- 4 files changed, 409 insertions(+), 582 deletions(-) diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 8e46497..ce8b564 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -57,11 +57,12 @@ _file_filter = FilePrefixFilter() from .dtm import classify_ground, create_dtm_fast from .visualizations import ( SharedDEM, - generate_hillshade, generate_slope, generate_aspect, generate_curvature, - generate_lrm, generate_openness, - generate_mslrm, generate_tpi, generate_sailore, + generate_hillshade, generate_slope, + generate_openness, + generate_mslrm, generate_sailore, generate_roughness, generate_wavelet, - generate_svf, generate_aniso_open, generate_paths, + generate_svf, generate_aniso_open, + generate_flow_accumulation, ) from .gpu import gpu_cleanup, num_gpus, restrict_gpus, safe_gpu_call from .ign import generate_ign_overlay @@ -74,19 +75,15 @@ from .rendering import tif_to_png VIZ_STEPS = [ ('hillshade', generate_hillshade), ('slope', generate_slope), - ('aspect', generate_aspect), - ('curvature', generate_curvature), - ('lrm', generate_lrm), + ('mslrm', generate_mslrm), + ('sailore', generate_sailore), ('pos_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=True, shared=shared)), ('neg_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=False, shared=shared)), - ('mslrm', generate_mslrm), - ('tpi', generate_tpi), - ('sailore', generate_sailore), - ('roughness', generate_roughness), ('svf', generate_svf), ('aniso_open', generate_aniso_open), - ('paths', generate_paths), + ('roughness', generate_roughness), ('wavelet', generate_wavelet), + ('flow_acc', generate_flow_accumulation), ('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 93ad007..b8ea731 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -81,13 +81,6 @@ _FRANCE_OUTLINE_L93 = np.array([ COLORMAPS = { # === Famille RELIEF : rouge=surélévation, bleu=dépression === # Diverging: rouge vif=positif, bleu vif=négatif, blanc=plat - 'curvature': { - 'cmap': 'bwr', - 'title': 'Courbure (Convexité/Concavité du terrain)', - 'legend': 'Changement de pente (1/m)\nRouge = Convexe (sommet de mur, levée)\nBleu = Concave (fond de fossé, dépression)', - 'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées', - 'vmin_mode': 'symmetric', 'sym_pct': (5, 95), - }, 'mslrm': { 'cmap': 'seismic', 'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)', @@ -95,20 +88,6 @@ COLORMAPS = { 'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, - 'lrm': { - 'cmap': 'seismic', - 'title': 'LRM - Local Relief Model (échelle unique 15m)', - 'legend': 'Écart local par rapport au terrain moyen (m)\nRouge = Surélévation (+{vmax:.2f}m)\nBleu = Dépression ({vmin:.2f}m)\nNoyau gaussien unique de 15m', - 'description': 'Micro-relief à 15m seulement — voir MSRM pour toutes les échelles', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, - 'tpi': { - 'cmap': 'seismic', - 'title': 'TPI - Topographic Position Index (4 échelles)', - 'legend': 'Position dans le paysage\nRouge = Plus haut que le voisinage (crête, plateau)\nBleu = Plus bas que le voisinage (fossé, vallée)\nCombine 4 échelles : 3m, 15m, 50m, 200m', - 'description': 'Identifie la position topographique — utile pour repérer crêtes vs vallées à grande échelle', - 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), - }, 'sailore': { 'cmap': 'seismic', 'title': 'SAILORE - LRM Auto-Adaptatif', @@ -123,14 +102,6 @@ COLORMAPS = { 'description': 'Openness avec pondération anisotropique — détecte mieux les structures alignées NW-SE et NE-SW', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, - 'paths': { - 'cmap': 'inferno', - 'title': 'Cheminement (chemins et sentiers)', - 'legend': 'Openness directionnelle maximale\nJaune/vif = Chemin ou sentier\nNoir = Terrain plat\n\nDétecte les structures linéaires dans toutes les directions', - 'description': 'Différence openness positive-négative maximale sur 8 directions — chemins, sentiers, ornières', - 'vmin_mode': 'percentile', 'vmin_pct': 5, - 'vmax_mode': 'percentile', 'vmax_pct': 99, - }, # === Famille OUVERTURE : séquentiel, toujours positif === 'positive_openness': { 'cmap': 'YlOrBr', @@ -160,7 +131,7 @@ COLORMAPS = { 'hillshade': { 'cmap': 'gray', 'title': 'Hillshade Multidirectionnel', - 'legend': 'Illumination combinée de 6 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre', + 'legend': 'Illumination combinée de 8 directions\nBlanc = Face éclairée | Noir = Zone d\'ombre', 'description': 'Ombres portées révélant micro-relief (murs, fossés, terrasses)', 'vmin_mode': 'percentile', 'vmin_pct': 1, 'vmax_mode': 'percentile', 'vmax_pct': 99, @@ -173,14 +144,6 @@ COLORMAPS = { 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'percentile', 'vmax_pct': 97, }, - 'aspect': { - 'cmap': 'twilight', - 'title': 'Aspect (Direction des pentes)', - 'legend': 'Direction vers laquelle le terrain descend\nCycle continu : Nord→Est→Sud→Ouest→Nord\nCouleurs perceptuellement uniformes (pas de saut de teinte)', - 'description': 'Orientation des pentes — utile pour distinguer structures des formes naturelles', - 'vmin_mode': 'fixed', 'vmin_val': 0, - 'vmax_mode': 'fixed', 'vmax_val': 360, - }, 'roughness': { 'cmap': 'plasma', 'title': 'Rugosité Multi-Échelle (3m + 15m)', @@ -191,11 +154,19 @@ COLORMAPS = { }, 'wavelet': { 'cmap': 'cividis', - 'title': 'Ondelette Mexican Hat (CWT multi-échelle)', - 'legend': 'Réponse de la transformée en ondelette\nÉchelles adaptées à la résolution\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires', - 'description': 'Transformée en ondelette 2D — excellente pour détecter structures circulaires', + 'title': 'Ondelette Mexican Hat + Gabor directionnelle (multi-échelle)', + 'legend': 'Réponse combinée Mexican Hat (circulaire) + Gabor (linéaire)\nÉchelles adaptées à la résolution\n4 orientations Gabor : 0°, 45°, 90°, 135°\n\nClair = Structure détectée\nSombre = Pas de structure', + 'description': 'Mexican Hat pour tumulus/enclos + Gabor pour chemins/murs/fossés', 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, + 'flow_acc': { + 'cmap': 'YlGn', + 'title': 'Accumulation d\'Écoulement (Flow Accumulation)', + 'legend': 'Log10 du nombre de cellules amont\nJaune = Fort accumulation (fossé, chenal, drainage)\nVert clair = Faible accumulation\n\nDétection fossés et linéaires hydrologiques', + 'description': 'Priority-flood + D8 — détecte fossés archéologiques et drainages', + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'percentile', 'vmax_pct': 98, + }, } # RGB entries (ortho/topo) are handled specially @@ -842,8 +813,7 @@ def generate_pdf_report(basename, vis_dir, pdf_dir, resolution): # Sort analysis files by archaeological priority order = ['mslrm', 'svf', 'negative_openness', 'positive_openness', 'aniso_open', 'sailore', 'hillshade_multi', - 'lrm', 'tpi', 'slope', 'curvature', 'aspect', - 'roughness', 'wavelet'] + 'flow_acc', 'slope', 'roughness', 'wavelet'] def sort_key(f): name = f.stem.lower() diff --git a/lidar_pipeline/tests/test_visualizations.py b/lidar_pipeline/tests/test_visualizations.py index 444352d..768f77f 100644 --- a/lidar_pipeline/tests/test_visualizations.py +++ b/lidar_pipeline/tests/test_visualizations.py @@ -47,52 +47,8 @@ class TestSlope: assert np.nanmax(data) <= 90 -class TestAspect: - def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_aspect - result = generate_aspect(synthetic_dem, "test", tmp_output_dir, 5.0) - assert result is not None - assert result.exists() - - def test_aspect_values_0_360(self, synthetic_dem, tmp_output_dir): - import rasterio - from lidar_pipeline.visualizations import generate_aspect - result = generate_aspect(synthetic_dem, "test", tmp_output_dir, 5.0) - with rasterio.open(result) as src: - data = src.read(1) - valid = data[~np.isnan(data)] - assert np.nanmin(valid) >= 0 - assert np.nanmax(valid) <= 360 - - -class TestCurvature: - def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_curvature - result = generate_curvature(synthetic_dem, "test", tmp_output_dir, 5.0) - assert result is not None - assert result.exists() - - # --- GPU-accelerated visualizations --- -class TestLRM: - def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_lrm - result = generate_lrm(synthetic_dem, "test", tmp_output_dir, 5.0) - assert result is not None - assert result.exists() - - def test_lrm_has_positive_negative(self, synthetic_dem, tmp_output_dir): - import rasterio - from lidar_pipeline.visualizations import generate_lrm - result = generate_lrm(synthetic_dem, "test", tmp_output_dir, 5.0) - with rasterio.open(result) as src: - data = src.read(1) - # LRM should have both positive and negative values - assert np.nanmax(data) > 0 - assert np.nanmin(data) < 0 - - class TestSVF: def test_generates_tif(self, synthetic_dem, tmp_output_dir): from lidar_pipeline.visualizations import generate_svf @@ -133,15 +89,6 @@ class TestMSLRM: assert result.exists() -class TestTPI: - def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_tpi - result = generate_tpi(synthetic_dem, "test", tmp_output_dir, 5.0) - assert result is not None - assert result.exists() - - - class TestSAILORE: def test_generates_tif(self, synthetic_dem, tmp_output_dir): from lidar_pipeline.visualizations import generate_sailore @@ -167,14 +114,6 @@ class TestRoughness: assert np.nanmin(data) >= 0 -class TestAnomalies: - def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_anomalies - result = generate_anomalies(synthetic_dem, "test", tmp_output_dir, 5.0) - assert result is not None - assert result.exists() - - class TestWavelet: def test_generates_tif(self, synthetic_dem, tmp_output_dir): from lidar_pipeline.visualizations import generate_wavelet @@ -183,50 +122,38 @@ class TestWavelet: assert result.exists() -class TestFlow: +class TestFlowAccumulation: def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_flow - result = generate_flow(synthetic_dem, "test", tmp_output_dir, 5.0) + from lidar_pipeline.visualizations import generate_flow_accumulation + result = generate_flow_accumulation(synthetic_dem, "test", tmp_output_dir, 5.0) assert result is not None assert result.exists() def test_flow_log_values(self, synthetic_dem, tmp_output_dir): import rasterio - from lidar_pipeline.visualizations import generate_flow - result = generate_flow(synthetic_dem, "test", tmp_output_dir, 5.0) + from lidar_pipeline.visualizations import generate_flow_accumulation + result = generate_flow_accumulation(synthetic_dem, "test", tmp_output_dir, 5.0) with rasterio.open(result) as src: data = src.read(1) - # log1p(x) >= 0 for x >= 0 + # log10(x) >= 0 for x >= 1 valid = data[~np.isnan(data)] assert np.nanmin(valid) >= 0 -class TestLocalDominance: - def test_generates_tif(self, synthetic_dem, tmp_output_dir): - from lidar_pipeline.visualizations import generate_local_dominance - result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0) - assert result is not None - assert result.exists() - assert result.suffix == ".tif" - - def test_dominance_values_0_1(self, synthetic_dem, tmp_output_dir): +class TestRayTrace: + def test_rays_are_traced(self, synthetic_dem, tmp_output_dir): + """Verify _ray_trace_horizons returns expected shapes.""" + from lidar_pipeline.visualizations import _ray_trace_horizons, _prepare_dem_for_raycast import rasterio - from lidar_pipeline.visualizations import generate_local_dominance - result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0) - with rasterio.open(result) as src: - data = src.read(1) - valid = data[~np.isnan(data)] - assert np.nanmin(valid) >= 0, "Local dominance should be >= 0" - assert np.nanmax(valid) <= 1, "Local dominance should be <= 1" - - def test_dominance_nan_mask_preserved(self, synthetic_dem, tmp_output_dir): - """Check that NaN zones from original DEM are preserved.""" - import rasterio - from lidar_pipeline.visualizations import generate_local_dominance - result = generate_local_dominance(synthetic_dem, "test", tmp_output_dir, 5.0) - # The synthetic DEM has no NaN, so this just verifies the output is valid - with rasterio.open(result) as src: - data = src.read(1) - # Shape should match input - assert data.shape[0] > 0 - assert data.shape[1] > 0 \ No newline at end of file + with rasterio.open(synthetic_dem) as src: + dem_np = src.read(1) + rows, cols = dem_np.shape + # Create a simple filled DEM for testing + import numpy as np + filled = np.nan_to_num(dem_np, nan=0) + # Test with numpy (no GPU) + pos, neg = _ray_trace_horizons( + filled, rows, cols, 5.0, n_dirs=4, max_dist=10, radii_m=[25, 50] + ) + assert pos.shape == (4, 2, rows, cols) + assert neg.shape == (4, 2, rows, cols) diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 96acd10..0c34a16 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -9,6 +9,7 @@ redundant I/O and computation. """ import logging +import math import time import warnings from pathlib import Path @@ -57,8 +58,8 @@ class SharedDEM: Reads the DEM once and lazily computes on first access: - NaN mask and filled DEM (avoids 20+ calls to _fill_nans) - - Gradient components (shared by hillshade, slope, aspect, curvature) - - LRM at 15m kernel (shared by lrm + anomalies) + - Gradient components (shared by hillshade, slope) + - LRM at 15m kernel (shared by mslrm + sailore) Attributes are computed lazily on first access to avoid computing data that is never used (e.g. LRM when only hillshade needs generation). @@ -75,7 +76,7 @@ class SharedDEM: # Lazy caches — computed on first access self._filled = None - self._gradient = None # (dy, dx, slope_rad, slope_deg, aspect) + self._gradient = None # (dy, dx, slope_rad, slope_deg) self._lrm_15 = None # GPU lazy caches @@ -110,11 +111,6 @@ class SharedDEM: self._ensure_gradient() return self._gradient[3] - @property - def aspect(self): - self._ensure_gradient() - return self._gradient[4] - @property def lrm_15(self): """LRM at 15m kernel — computed lazily.""" @@ -134,8 +130,7 @@ class SharedDEM: dx = np.gradient(self.filled, self.resolution, axis=1) slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) slope_deg = np.degrees(slope_rad) - aspect = np.mod(np.degrees(np.arctan2(dy, dx)), 360) - self._gradient = (dy, dx, slope_rad, slope_deg, aspect) + self._gradient = (dy, dx, slope_rad, slope_deg) @property def filled_gpu(self): @@ -270,6 +265,120 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs): return result +# ============================================================ +# Shared ray-tracing core +# ============================================================ + +def _prepare_dem_for_raycast(dem_file, shared, resolution): + """Load DEM and prepare padded array for ray-tracing. + + Returns (dem_gpu_or_cpu, dem_np, rows, cols, res, nan_mask, + transform, crs, padded) ready for ray-tracing. + """ + if shared: + dem_np = shared.dem_np + nan_mask = shared.nan_mask + transform = shared.transform + crs = shared.crs + dem = to_gpu(shared.filled) if _gpu_mod.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 = to_gpu(filled) if _gpu_mod.HAS_GPU else filled + res = resolution + rows, cols = dem_np.shape + return dem, dem_np, rows, cols, res, nan_mask, transform, crs + + +def _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m=None): + """Core ray-tracing: compute max zenith/nadir angles per direction and radius. + + For each pixel, in each direction, traces rays outward up to max_dist steps, + recording the max upward angle (positive openness) and max downward angle + (negative openness) reached at each radius checkpoint. + + Args: + dem: GPU or CPU filled DEM array (rows, cols). + rows, cols: dimensions. + res: resolution in m/px. + n_dirs: number of directions. + max_dist: max ray steps. + radii_m: list of radii in meters to record checkpoints. + If None, records only at max_dist. + + Returns: + pos_angles: array of shape (n_dirs, n_radii, rows, cols) — max zenith angles + neg_angles: array of shape (n_dirs, n_radii, rows, cols) — max nadir angles + """ + angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) + dx_dir = np.cos(angles) + dy_dir = np.sin(angles) + + if radii_m is not None: + radii_steps = [min(int(r / res), max_dist) for r in radii_m] + n_radii = len(radii_m) + else: + radii_steps = [max_dist] + n_radii = 1 + + padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) + + pos_angles = xp.zeros((n_dirs, n_radii, rows, cols)) + neg_angles = xp.zeros((n_dirs, n_radii, rows, cols)) + + for d_idx in range(n_dirs): + ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] + + # Pre-compute valid steps for this direction + valid_steps = [] + for step in range(1, max_dist + 1): + px = int(round(ddx * step)) + py = int(round(ddy * step)) + dist_m = math.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) + if dist_m < res * 0.5: + continue + valid_steps.append((step, px, py, dist_m)) + + # Running max angles per radius + running_pos = xp.zeros((n_radii, rows, cols)) + running_neg = xp.zeros((n_radii, rows, cols)) + # Track which radius checkpoints have been passed + radii_remaining = set(range(n_radii)) + + for step, px, py, dist_m in valid_steps: + elev_diff = padded[max_dist + py:max_dist + py + rows, + max_dist + px:max_dist + px + cols] - dem + + # Positive: angle to terrain above viewer + pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) + # Negative: angle to terrain below viewer + neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) + + # Update running max for all radius checkpoints still active + for r_idx in radii_remaining: + pos_angle_safe = xp.nan_to_num(pos_angle, nan=0) + neg_angle_safe = xp.nan_to_num(neg_angle, nan=0) + running_pos[r_idx] = xp.where(xp.isnan(pos_angle), running_pos[r_idx], + xp.maximum(running_pos[r_idx], pos_angle_safe)) + running_neg[r_idx] = xp.where(xp.isnan(neg_angle), running_neg[r_idx], + xp.maximum(running_neg[r_idx], neg_angle_safe)) + + # Check which radii have been passed + new_remaining = set() + for r_idx in radii_remaining: + if step < radii_steps[r_idx]: + new_remaining.add(r_idx) + radii_remaining = new_remaining + if not radii_remaining: + break + + pos_angles[d_idx] = running_pos + neg_angles[d_idx] = running_neg + + return pos_angles, neg_angles + + # ============================================================ # Core terrain visualizations # ============================================================ @@ -373,174 +482,42 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None): return None -def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None): - """Generate aspect (slope orientation) map — GPU if available.""" - gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Aspect (Orientation){gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_aspect.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - aspect = shared.aspect - nan_mask = shared.nan_mask - if _gpu_mod.HAS_GPU: - aspect = to_gpu(aspect) - else: - dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) - dy, dx = xp.gradient(dem) - aspect = xp.arctan2(dy, dx) * 180 / xp.pi - aspect = xp.mod(aspect, 360) - nan_mask = np.isnan(dem_np) - _save_tif(output, to_cpu(aspect) if _gpu_mod.HAS_GPU else aspect, transform, crs, nan_mask=nan_mask) - logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur aspect: {e}", exc_info=True) - return None - - -def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None): - """Generate curvature (terrain concavity/convexity) map — GPU if available.""" - gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Courbure (Curvature){gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_curvature.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - dx = shared.dx - dy = shared.dy - nan_mask = shared.nan_mask - if _gpu_mod.HAS_GPU: - dx = to_gpu(dx) - dy = to_gpu(dy) - else: - dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) - dy, dx = xp.gradient(dem) - nan_mask = np.isnan(dem_np) - d2z_dx2 = xp.gradient(dx, axis=1) - d2z_dy2 = xp.gradient(dy, axis=0) - curvature = (d2z_dx2 + d2z_dy2) / 2 - _save_tif(output, to_cpu(curvature), transform, crs, nan_mask=nan_mask) - logger.info(f" ✓ Courbure terminée ({time.time()-t0:.1f}s){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur curvature: {e}", exc_info=True) - return None - - - # ============================================================ # GPU-accelerated visualizations # ============================================================ -def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None): - """Local Relief Model - deviation from local mean (GPU if available). - - Kernel sigma adapts to resolution: finer kernel at higher resolution - to capture micro-relief details. At 0.5m/px: 15m, at 0.2m/px: ~5m. - """ - gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Local Relief Model{gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_lrm.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - lrm = shared.lrm_15.copy() - else: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) - # Adapt sigma to resolution: standard 15m at 0.5m, finer at higher res - sigma_m = max(5.0, 15.0 * 0.5 / resolution) - logger.info(f" LRM sigma={sigma_m:.1f}m (résolution {resolution}m/px)") - local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_m / resolution) - lrm = dem_np - local_mean - lrm[nan_mask] = np.nan - _save_tif(output, lrm.astype(np.float32), transform, crs) - logger.info(f" ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur LRM: {e}", exc_info=True) - return None - - def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): - """Sky-View Factor - ray-tracing on 16 azimuths (GPU if available). + """Sky-View Factor - ray-tracing on 16 azimuths, multi-radius (GPU if available). - For each pixel, trace rays in N directions, find the max horizon - angle in each direction, then SVF = (1/N) * sum(cos²(horizon_angle)). - Valleys/crevices have low SVF (obstructed sky), ridges/peaks have high SVF. + Traces rays in 16 directions at 3 radii (25, 50, 100m) and combines + with weights favoring medium range for archaeological feature detection. + SVF = (1/N) * sum(cos²(horizon_angle)). Valleys/crevices have low SVF, + ridges/peaks have high SVF. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Sky-View Factor (ray-tracing){gpu_tag}...") + logger.info(f" → Sky-View Factor (ray-tracing multi-rayon){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_svf.tif" try: - if shared: - transform = shared.transform - crs = shared.crs - dem_np = shared.dem_np - rows, cols = dem_np.shape - res = resolution - dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled - nan_mask = shared.nan_mask - else: - dem_np, transform, crs = _read_dem(dem_file) - rows, cols = dem_np.shape - res = resolution - nan_mask = np.isnan(dem_np) - filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled + dem, dem_np, rows, cols, res, nan_mask, transform, crs = \ + _prepare_dem_for_raycast(dem_file, shared, resolution) - n_dirs = 16 - angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx_dir = np.cos(angles) - dy_dir = np.sin(angles) - # Cap max_dist to avoid excessive computation at high resolution - # 100m radius is sufficient; at 0.2m that's 500 steps which is very slow + radii_m = [25, 50, 100] + radius_weights = [0.3, 0.4, 0.3] # Medium range weighted more max_dist = min(int(100 / res), 300) + n_dirs = 16 - padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - svf = xp.zeros_like(dem) + pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m) - for d_idx in range(n_dirs): - ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] - horizon = xp.zeros_like(dem) + # SVF per radius: mean of cos²(horizon) across directions + svf_combined = xp.zeros_like(dem) + for r_idx in range(len(radii_m)): + horizon = xp.maximum(pos_angles[:, r_idx], neg_angles[:, r_idx]) + svf_r = xp.mean(xp.cos(horizon) ** 2, axis=0) + svf_combined += svf_r * radius_weights[r_idx] - # Pre-compute all valid steps for this direction - valid_steps = [] - 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 - valid_steps.append((step, px, py, dist_m)) - - # Batch all shifts into a single array for vectorized max computation - for step, px, py, dist_m in valid_steps: - elev_diff = padded[max_dist + py:max_dist + py + rows, - max_dist + px:max_dist + px + cols] - dem - angle = xp.arctan2(elev_diff, dist_m) - horizon = xp.where(xp.isnan(angle), horizon, - xp.maximum(horizon, xp.nan_to_num(angle, nan=0))) - - # SVF uses cos²(horizon angle) — fraction of visible sky - svf += xp.cos(horizon) ** 2 - - svf /= n_dirs - svf_np = to_cpu(svf).astype(np.float32) + svf_np = to_cpu(svf_combined).astype(np.float32) svf_np[nan_mask] = np.nan _save_tif(output, svf_np, transform, crs) logger.info(f" ✓ SVF terminé ({time.time()-t0:.1f}s){gpu_tag}") @@ -551,72 +528,45 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None): def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None): - """Positive/Negative Openness - true zenith/nadir angle computation (GPU if available). + """Positive/Negative Openness - multi-radius ray-tracing with std normalization. - For each pixel, in 8 directions (N, NE, E, SE, S, SW, W, NW): - - Positive openness: max zenith angle (angle from vertical to highest visible terrain) - - Negative openness: max nadir angle (angle from vertical down to lowest terrain) - Result is averaged across all 8 directions. - Ray radius adapts to resolution: 100m for better detection of large enclosures. + Traces rays in 8 directions at 3 radii (25, 50, 100m). + Results are combined with equal weight across radii, then normalized + by standard deviation for cross-tile comparability. """ name = "positive_openness" if positive else "negative_openness" gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing){gpu_tag}...") + logger.info(f" → {name.replace('_', ' ').title()} (ray-tracing multi-rayon){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_{name}.tif" try: - if shared: - transform = shared.transform - crs = shared.crs - dem_np = shared.dem_np - rows, cols = dem_np.shape - res = resolution - dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled - nan_mask = shared.nan_mask - else: - dem_np, transform, crs = _read_dem(dem_file) - rows, cols = dem_np.shape - res = resolution - nan_mask = np.isnan(dem_np) - filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled + dem, dem_np, rows, cols, res, nan_mask, transform, crs = \ + _prepare_dem_for_raycast(dem_file, shared, resolution) - n_dirs = 8 - angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx_dir = np.cos(angles) - dy_dir = np.sin(angles) + radii_m = [25, 50, 100] max_dist = min(int(100 / res), 300) + n_dirs = 8 - padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - openness_sum = xp.zeros_like(dem) + pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m) - for d_idx in range(n_dirs): - ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] - max_angle = xp.zeros_like(dem) + # Select positive or negative + if positive: + angles = pos_angles + else: + angles = neg_angles - 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 - - if positive: - angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) - else: - angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) - - max_angle = xp.where(xp.isnan(angle), max_angle, - xp.maximum(max_angle, xp.nan_to_num(angle, nan=0))) - - openness_sum += max_angle - - openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32) + # Mean across directions and radii (equal weight) + openness = xp.mean(angles, axis=(0, 1)) + openness_result = to_cpu(xp.degrees(openness)).astype(np.float32) openness_result[nan_mask] = np.nan + + # Std normalization for cross-tile comparability + valid = openness_result[~nan_mask] + if len(valid) > 0: + std_val = max(np.nanstd(valid), 0.01) + openness_result = openness_result / std_val + _save_tif(output, openness_result, transform, crs) logger.info(f" ✓ {name} terminé ({time.time()-t0:.1f}s){gpu_tag}") return output @@ -694,68 +644,6 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None): return None -def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None): - """Multi-Scale Topographic Position Index (GPU if available). - - TPI = elevation - mean(neighborhood). - Computed at 4 scales with std normalization and weighted combination. - Weights favor fine and medium scales (archaeologically relevant). - """ - gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → TPI multi-échelle{gpu_tag}...") - t0 = time.time() - output = vis_dir / f"{basename}_tpi.tif" - - try: - if shared: - transform = shared.transform - crs = shared.crs - dem_np = shared.dem_np - nan_mask = shared.nan_mask - else: - dem_np, transform, crs = _read_dem(dem_file) - nan_mask = np.isnan(dem_np) - - # 4 scales: fine (3m), medium (15m), broad (50m), landscape (200m) - scales_m = [3, 15, 50, 200] - weights = [1.5, 2.0, 1.2, 0.5] # Favor medium scales (ditches, enclosures) - - tpi_stack = [] - for scale_m, weight in zip(scales_m, weights): - size = max(3, int(scale_m / resolution)) - if size % 2 == 0: - size += 1 - if shared: - local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=size) - else: - local_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=size) - tpi = dem_np - local_mean - tpi[nan_mask] = np.nan - # Std normalization — preserves sign and contrast better than z-score - valid = tpi[~nan_mask] - tpi_std = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 - tpi = tpi / tpi_std - tpi_stack.append(tpi.astype(np.float32)) - - # Weighted combination - tpi_array = np.array(tpi_stack) - weights_3d = np.array(weights)[:, np.newaxis, np.newaxis] - with np.errstate(invalid='ignore', divide='ignore'): - with warnings.catch_warnings(): - warnings.filterwarnings('ignore', message='Mean of empty slice') - tpi_combined = np.nansum(tpi_array * weights_3d, axis=0) / np.sum(weights) - tpi_combined[nan_mask] = np.nan - - _save_tif(output, tpi_combined.astype(np.float32), transform, crs) - logger.info(f" ✓ TPI terminé ({time.time()-t0:.1f}s){gpu_tag}") - return output - except Exception as e: - logger.error(f" ✗ Erreur TPI: {e}", exc_info=True) - return None - - - - # ============================================================ # SAILORE # ============================================================ @@ -765,6 +653,7 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None): Kernel size adapts to local slope: flat areas get larger kernels, steep areas get smaller kernels. Scales adapt to resolution. + Reuses shared.lrm_15 when available to avoid recomputation. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" logger.info(f" → SAILORE (LRM adaptatif){gpu_tag}...") @@ -791,21 +680,28 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None): sigma_max_m = 25.0 # 25m — broad relief sigma_min = sigma_min_m / resolution sigma_max = sigma_max_m / resolution - sigma_mid = (sigma_min + sigma_max) / 2 slope_norm = np.clip(slope_deg / 30.0, 0, 1) + # LRM fine (2m) — always compute if shared: lrm_fine = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_min) else: lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min) lrm_fine[nan_mask] = np.nan - if shared: - lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) + # LRM medium (13.5m) — reuse shared.lrm_15 (σ=15m) when available + sigma_mid = (sigma_min + sigma_max) / 2 + if shared and abs(15.0 / resolution - sigma_mid) < 2.0 / resolution: + # shared.lrm_15 is close enough to medium scale + lrm_medium = shared.lrm_15.copy() else: - lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) + if shared: + lrm_medium = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_mid) + else: + lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_mid) lrm_medium[nan_mask] = np.nan + # LRM coarse (25m) — always compute if shared: lrm_coarse = dem_np - _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_max) else: @@ -902,22 +798,39 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): # ============================================================ -# Wavelet +# Wavelet (Mexican Hat + Directional Gabor) # ============================================================ +def _gabor_kernel_2d(size, sigma, wavelength, theta): + """Create a 2D Gabor kernel. + + Args: + size: kernel size (odd integer) + sigma: standard deviation + wavelength: wavelength of sinusoid + theta: orientation angle in radians (0 = horizontal) + """ + center = size // 2 + y, x = np.ogrid[-center:center+1, -center:center+1] + # Rotate coordinates + x_theta = x * np.cos(theta) + y * np.sin(theta) + y_theta = -x * np.sin(theta) + y * np.cos(theta) + sigma_sq = 2 * sigma * sigma + kernel = np.exp(-(x_theta**2 + y_theta**2) / sigma_sq) + kernel *= np.cos(2 * np.pi * x_theta / wavelength) + return kernel + + def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): - """Mexican Hat wavelet multi-scale analysis (GPU if available). + """Multi-scale wavelet analysis: Mexican Hat + Directional Gabor (GPU if available). - CWT 2D at multiple scales adapted to resolution. - - At 0.5m/px: [1, 2, 5, 10, 20, 50, 100]m - - At 0.2m/px: [0.5, 1, 2, 5, 10, 20, 50, 100]m - - Higher resolution = more fine scales available - - Uses std normalization per scale and weighted combination - with emphasis on archaeologically relevant scales (2-50m). + Mexican Hat (radial): detects circular features (tumulus, enclos ronds). + Gabor (directional): detects linear features (chemins, murs, fossés). + 4 Gabor orientations (0°, 45°, 90°, 135°) at key archaeological scales. + Both combined with RMS for orientation-invariant detection. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Ondelette Mexican Hat multi-échelle{gpu_tag}...") + logger.info(f" → Ondelette Mexican Hat + Gabor directionnelle{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_wavelet.tif" @@ -933,28 +846,25 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np.astype(np.float64)) - # Adapt scales to resolution: finer scales available at higher resolution + # --- Mexican Hat scales --- min_scale = max(resolution * 2, 1.0) candidate_scales = [0.5, 1, 2, 5, 10, 20, 50, 100] - scales = [s for s in candidate_scales if s >= min_scale] + mex_scales = [s for s in candidate_scales if s >= min_scale] - # Weights favor archaeological scales (2-50m: ditches, enclosures, tumulus) - scale_weights = { - 0.5: 0.6, # Fine texture - 1.0: 0.8, # Micro-relief - 2.0: 1.5, # Small ditches, paths — key scale - 5.0: 2.0, # Fossés, small enclosures — key archaeological scale - 10.0: 1.8, # Medium structures - 20.0: 1.5, # Large enclosures, tumulus - 50.0: 1.0, # Very large enclosures - 100.0: 0.6, # Landscape-level features + mex_weights_map = { + 0.5: 0.6, 1.0: 0.8, 2.0: 1.5, 5.0: 2.0, + 10.0: 1.8, 20.0: 1.5, 50.0: 1.0, 100.0: 0.6, } - weights = np.array([scale_weights.get(s, 1.0) for s in scales]) + mex_weights = np.array([mex_weights_map.get(s, 1.0) for s in mex_scales]) + + logger.info(f" Échelles CWT: {mex_scales}m (résolution {resolution}m/px)") + + from scipy.ndimage import gaussian_laplace, convolve - logger.info(f" Échelles CWT: {scales}m (résolution {resolution}m/px)") wavelet_stack = [] - for scale_m in scales: + # Mexican Hat (radial) — multi-scale + for scale_m in mex_scales: sigma_px = scale_m / resolution if _gpu_mod.HAS_GPU: try: @@ -962,27 +872,53 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px) response = to_cpu(response) except Exception: - from scipy.ndimage import gaussian_laplace response = -gaussian_laplace(filled, sigma=sigma_px) else: - from scipy.ndimage import gaussian_laplace response = -gaussian_laplace(filled, sigma=sigma_px) response[nan_mask] = np.nan - - # Std normalization: scale by standard deviation to make scales comparable valid = response[~nan_mask] std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 response = response / std_val wavelet_stack.append(response) - # Weighted RMS: sqrt(sum(w * x²) / sum(w)) - # Preserves contrast at key archaeological scales + # Gabor (directional) — 4 orientations at 3 key scales + gabor_scales_m = [5, 10, 20] # Key archaeological scales for linear features + gabor_orientations = [0, np.pi/4, np.pi/2, 3*np.pi/4] # 0°, 45°, 90°, 135° + gabor_scale_weights = {5: 2.0, 10: 1.8, 20: 1.5} + + for scale_m in gabor_scales_m: + sigma_px = max(2, scale_m / resolution / 3) # Sigma relative to wavelength + wavelength_px = max(3, scale_m / resolution) + kernel_size = max(5, int(wavelength_px * 2.5)) + if kernel_size % 2 == 0: + kernel_size += 1 + + for theta in gabor_orientations: + kernel = _gabor_kernel_2d(kernel_size, sigma_px, wavelength_px, theta) + # Normalize kernel + kernel = kernel / max(np.abs(kernel).max(), 1e-10) + response = convolve(filled, kernel, mode='constant', cval=0) + response[nan_mask] = np.nan + valid = response[~nan_mask] + std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 + response = response / std_val + wavelet_stack.append(response) + + # Build weights: Mexican Hat weights + Gabor weights + gabor_weights_list = [] + for scale_m in gabor_scales_m: + w = gabor_scale_weights.get(scale_m, 1.0) + gabor_weights_list.extend([w] * len(gabor_orientations)) + gabor_weights = np.array(gabor_weights_list) + all_weights = np.concatenate([mex_weights, gabor_weights]) + + # Weighted RMS combination stack = np.array(wavelet_stack) - weights_3d = weights[:, np.newaxis, np.newaxis] + weights_3d = all_weights[:, np.newaxis, np.newaxis] with np.errstate(invalid='ignore', divide='ignore'): with warnings.catch_warnings(): warnings.filterwarnings('ignore', message='Mean of empty slice') - combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(weights)) + combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(all_weights)) combined[nan_mask] = np.nan _save_tif(output, combined.astype(np.float32), transform, crs) @@ -994,187 +930,184 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): # ============================================================ -# Anisotropic Openness -# ============================================================ -# Path Detection (chemins et sentiers) +# Flow Accumulation # ============================================================ -def generate_paths(dem_file, basename, vis_dir, resolution, shared=None): - """Cheminement — openness directionnelle maximale pour détecter chemins et sentiers. +def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=None): + """Flow Accumulation — priority-flood sink filling + D8 accumulation (GPU if available). - Pour chaque direction (8 directions), calcule openness positive - négative, - puis prend le maximum sur toutes les directions. Les chemins et sentiers - ressortent en valeurs élevées quelle que soit leur orientation. + Detects channels, ditches, and drainage paths by computing how many + upstream cells flow through each cell. Archaeological ditches and + natural drainage features both accumulate high flow values. - Contrairement à l'openness anisotropique qui privilégie NW-SE et NE-SW, - cette visualisation traite toutes les directions de manière égale et - combine positive et négative en une seule image. Les chemins perpendiculaires - à une direction auront une forte différence dans cette direction, donc - le maximum sur toutes les directions les fait ressortir. + Uses log10 transformation for visualization. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Cheminement (chemins et sentiers){gpu_tag}...") + logger.info(f" → Accumulation d'écoulement (flow accumulation){gpu_tag}...") t0 = time.time() - output = vis_dir / f"{basename}_paths.tif" + output = vis_dir / f"{basename}_flow_acc.tif" try: if shared: transform = shared.transform crs = shared.crs dem_np = shared.dem_np - rows, cols = dem_np.shape - res = resolution - dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled nan_mask = shared.nan_mask + filled = shared.filled else: dem_np, transform, crs = _read_dem(dem_file) - rows, cols = dem_np.shape - res = resolution nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled - n_dirs = 8 - angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx_dir = np.cos(angles) - dy_dir = np.sin(angles) - max_dist = min(int(100 / res), 300) + # Priority-flood sink filling (Wang & Liu 2006, O(n log n)) + from heapq import heappush, heappop + rows, cols = filled.shape + dem_filled = filled.copy() - padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - max_diff = xp.full_like(dem, -1e6) # Will track max over all directions + # Use heap for priority-flood + visited = np.zeros((rows, cols), dtype=bool) + heap = [] + # Seed with all border cells + for x in range(cols): + heappush(heap, (dem_filled[0, x], 0, x)) + heappush(heap, (dem_filled[rows-1, x], rows-1, x)) + visited[0, x] = True + visited[rows-1, x] = True + for y in range(1, rows-1): + heappush(heap, (dem_filled[y, 0], y, 0)) + heappush(heap, (dem_filled[y, cols-1], y, cols-1)) + visited[y, 0] = True + visited[y, cols-1] = True - for d_idx in range(n_dirs): - ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] + while heap: + elev, cy, cx = heappop(heap) + for dy, dx in [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]: + ny, nx = cy + dy, cx + dx + if 0 <= ny < rows and 0 <= nx < cols and not visited[ny, nx]: + if dem_filled[ny, nx] > elev: + dem_filled[ny, nx] = elev + heappush(heap, (dem_filled[ny, nx], ny, nx)) + visited[ny, nx] = True - # Positive openness: max zenith angle in this direction - max_pos_angle = xp.zeros_like(dem) - # Negative openness: max nadir angle in this direction - max_neg_angle = xp.zeros_like(dem) + logger.info(f" ✓ Sink filling terminé ({time.time()-t0:.1f}s)") - 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 + # D8 flow direction and accumulation + flow_acc = np.zeros((rows, cols), dtype=np.int64) + flow_dir = np.zeros((rows, cols), dtype=np.int8) - 1 - elev_diff = padded[max_dist + py:max_dist + py + rows, - max_dist + px:max_dist + px + cols] - dem + # D8 neighbors (ordered by angle) + neighbors = [ + (-1, 0), (-1, 1), ( 0, 1), ( 1, 1), + ( 1, 0), ( 1, -1), ( 0, -1), (-1, -1) + ] - # Positive: angle to terrain above viewer - pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) - max_pos_angle = xp.where(xp.isnan(pos_angle), max_pos_angle, - xp.maximum(max_pos_angle, xp.nan_to_num(pos_angle, nan=0))) + # Process cells in ascending elevation order for correct accumulation + flat_idx = np.argsort(dem_filled.ravel()) + flow_acc_flat = flow_acc.ravel() - # Negative: angle to terrain below viewer - neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) - max_neg_angle = xp.where(xp.isnan(neg_angle), max_neg_angle, - xp.maximum(max_neg_angle, xp.nan_to_num(neg_angle, nan=0))) + for idx in flat_idx: + y = idx // cols + x = idx % cols - # Difference highlights linear features perpendicular to this direction - diff = max_pos_angle - max_neg_angle - max_diff = xp.maximum(max_diff, diff) + # Find steepest downslope neighbor + max_slope = -np.inf + best_dir = -1 + for di, (dy, dx) in enumerate(neighbors): + ny, nx = y + dy, x + dx + if 0 <= ny < rows and 0 <= nx < cols: + slope = (dem_filled[y, x] - dem_filled[ny, nx]) + dist = math.sqrt(dy**2 + dx**2) + slope_per_m = slope / dist + if slope_per_m > max_slope: + max_slope = slope_per_m + best_dir = di - paths_result = to_cpu(max_diff).astype(np.float32) - paths_result[nan_mask] = np.nan - _save_tif(output, paths_result, transform, crs) - logger.info(f" ✓ Cheminement terminé ({time.time()-t0:.1f}s){gpu_tag}") + flow_dir[y, x] = best_dir + flow_acc[y, x] = 1 # Count self + + # Accumulate flow (upstream to downstream) + # Process in reverse elevation order (highest first) + for idx in reversed(flat_idx): + y = idx // cols + x = idx % cols + d = flow_dir[y, x] + if d >= 0: + dy, dx = neighbors[d] + ny, nx = y + dy, x + dx + if 0 <= ny < rows and 0 <= nx < cols: + flow_acc[ny, nx] += flow_acc[y, x] + + logger.info(f" ✓ D8 accumulation terminé ({time.time()-t0:.1f}s)") + + # Log transform for visualization + flow_result = np.log10(np.maximum(flow_acc.astype(np.float32), 1.0)) + flow_result[nan_mask] = np.nan + + _save_tif(output, flow_result, transform, crs) + logger.info(f" ✓ Flow accumulation terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: - logger.error(f" ✗ Erreur cheminement: {e}", exc_info=True) + logger.error(f" ✗ Erreur flow accumulation: {e}", exc_info=True) return None +# ============================================================ +# Anisotropic Openness # ============================================================ def generate_aniso_open(dem_file, basename, vis_dir, resolution, shared=None): - """Anisotropic Openness - weighted directional openness emphasizing oblique directions (GPU if available). + """Anisotropic Openness - weighted directional openness, multi-radius (GPU if available). Computes positive and negative openness with anisotropic weighting: NW/SE directions weighted more heavily to enhance detection of structures aligned NE-SW (common in French archaeological sites: villas, enclosures). - The anisotropic weighting makes subtle linear features more visible than - standard isotropic openness which averages all directions equally. + Multi-radius (25, 50, 100m) with equal weight, results std-normalized + for cross-tile comparability. """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Openness Anisotropique{gpu_tag}...") + logger.info(f" → Openness Anisotropique (multi-rayon){gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_aniso_open.tif" try: - if shared: - transform = shared.transform - crs = shared.crs - dem_np = shared.dem_np - rows, cols = dem_np.shape - res = resolution - dem = to_gpu(shared.filled) if _gpu_mod.HAS_GPU else shared.filled - nan_mask = shared.nan_mask - else: - dem_np, transform, crs = _read_dem(dem_file) - rows, cols = dem_np.shape - res = resolution - nan_mask = np.isnan(dem_np) - filled, _ = _fill_nans(dem_np) - dem = to_gpu(filled) if _gpu_mod.HAS_GPU else filled + dem, dem_np, rows, cols, res, nan_mask, transform, crs = \ + _prepare_dem_for_raycast(dem_file, shared, resolution) n_dirs = 8 - angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) - dx_dir = np.cos(angles) - dy_dir = np.sin(angles) - # Anisotropic weights: emphasize NW-SE and NE-SW directions - # These orientations are most productive for detecting archaeological features - # aligned with Roman and medieval settlement patterns in France weights = np.array([1.0, 1.5, 1.0, 1.5, 1.0, 1.5, 1.0, 1.5]) + radii_m = [25, 50, 100] max_dist = min(int(100 / res), 300) - padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - pos_sum = xp.zeros_like(dem) - neg_sum = xp.zeros_like(dem) - weight_total = 0.0 + pos_angles, neg_angles = _ray_trace_horizons(dem, rows, cols, res, n_dirs, max_dist, radii_m) - for d_idx in range(n_dirs): - ddx, ddy = dx_dir[d_idx], dy_dir[d_idx] - w = weights[d_idx] - weight_total += w + # Weighted combination across directions and radii + weight_total = np.sum(weights) + n_radii = len(radii_m) - max_pos_angle = xp.zeros_like(dem) - max_neg_angle = xp.zeros_like(dem) + pos_combined = xp.zeros_like(dem) + neg_combined = xp.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 + for r_idx in range(n_radii): + for d_idx in range(n_dirs): + w = weights[d_idx] + pos_combined += pos_angles[d_idx, r_idx] * w / (n_radii * weight_total) + neg_combined += neg_angles[d_idx, r_idx] * w / (n_radii * weight_total) - elev_diff = padded[max_dist + py:max_dist + py + rows, - max_dist + px:max_dist + px + cols] - dem - - # Positive openness: max zenith angle - pos_angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) - max_pos_angle = xp.where(xp.isnan(pos_angle), max_pos_angle, - xp.maximum(max_pos_angle, xp.nan_to_num(pos_angle, nan=0))) - - # Negative openness: max nadir angle - neg_angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) - max_neg_angle = xp.where(xp.isnan(neg_angle), max_neg_angle, - xp.maximum(max_neg_angle, xp.nan_to_num(neg_angle, nan=0))) - - pos_sum += max_pos_angle * w - neg_sum += max_neg_angle * w - - # Combined: positive minus negative openness (anisotropic) - pos_avg = pos_sum / weight_total - neg_avg = neg_sum / weight_total - aniso_result = to_cpu(xp.degrees(pos_avg - neg_avg)).astype(np.float32) + aniso_result = to_cpu(xp.degrees(pos_combined - neg_combined)).astype(np.float32) aniso_result[nan_mask] = np.nan + + # Std normalization for cross-tile comparability + valid = aniso_result[~nan_mask] + if len(valid) > 0: + std_val = max(np.nanstd(valid), 0.01) + aniso_result = aniso_result / std_val + _save_tif(output, aniso_result, transform, crs) logger.info(f" ✓ Openness anisotropique terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: logger.error(f" ✗ Erreur openness anisotropique: {e}", exc_info=True) - return None \ No newline at end of file + return None