diff --git a/lidar_pipeline/gpu.py b/lidar_pipeline/gpu.py index a7e97d3..f2f5ddd 100644 --- a/lidar_pipeline/gpu.py +++ b/lidar_pipeline/gpu.py @@ -24,6 +24,12 @@ _gpu_mem_gb = 0 _best_gpu_id: int | None = None _gpu_reason = None +# GPU restriction from -g flag (host-level indices) +_restricted_gpu_ids: list[int] | None = None + +# Discovered GPU candidates (populated by _pick_gpu) +_gpu_candidates: list[dict] = [] + def _pick_gpu() -> list: """List all GPUs from the system, sorted by compute capability (highest first).""" @@ -75,41 +81,104 @@ _cp_ndimage = None _gpu_initialized = False +def _filter_candidates(gpus: list) -> list: + """Filter GPU candidates by CUDA_VISIBLE_DEVICES and _restricted_gpu_ids. + + nvidia-smi lists ALL GPUs even when CUDA_VISIBLE_DEVICES is set + (driver 580.x behavior), so we must filter manually. + """ + # Filter by CUDA_VISIBLE_DEVICES if set + cuda_visible = os.environ.get('CUDA_VISIBLE_DEVICES') + if cuda_visible is not None: + try: + visible = {int(i.strip()) for i in cuda_visible.split(',')} + gpus = [g for g in gpus if g[0] in visible] + except ValueError: + pass + + # Filter by programmatic restriction (-g 0, -g 0,2) + if _restricted_gpu_ids is not None: + allowed = set(_restricted_gpu_ids) + gpus = [g for g in gpus if g[0] in allowed] + + return gpus + + def _init_gpu(): """Lazily initialize CuPy on first GPU use. - Uses a subprocess to test each GPU (highest compute capability first). - The subprocess sets CUDA_VISIBLE_DEVICES before importing CuPy and - runs a warm-up kernel. First GPU that passes wins. - - This is necessary because CUDA_VISIBLE_DEVICES must be set in the - process environment BEFORE CuPy imports, not via os.environ in Python. + 1. Filters candidates by CUDA_VISIBLE_DEVICES and _restricted_gpu_ids + 2. If CUDA_VISIBLE_DEVICES is already set (e.g. run.sh -g 0): + import CuPy directly (no subprocess test needed) + 3. Otherwise: test each GPU in a subprocess, pick the first that works """ global _xp, _cp, _cp_ndimage, _gpu_initialized, HAS_GPU, _best_gpu_id, _gpu_name, _gpu_mem_gb if _gpu_initialized: return _gpu_initialized = True - if not _candidate_gpus: + candidates = _filter_candidates(_candidate_gpus) + if not candidates: + logger.info("Pas de GPU utilisable — mode CPU uniquement") _xp = np _cp = None _cp_ndimage = None HAS_GPU = False return - # Test each GPU in a subprocess + cuda_visible = os.environ.get('CUDA_VISIBLE_DEVICES') + + if cuda_visible is not None: + # CUDA_VISIBLE_DEVICES already set (e.g. by run.sh -g 0). + # Pick the best visible GPU and import CuPy directly. + idx, name, cap_str, mem_mi, score, major = candidates[0] + try: + import cupy as _real_cupy + import cupyx.scipy.ndimage as _real_cupy_ndimage + + # Warm-up kernel to verify GPU works + x = _real_cupy.array([1.0, 2.0], dtype=_real_cupy.float32) + s = _real_cupy.sum(x).get() + if s != 3.0: + raise RuntimeError("GPU warm-up failed") + + _best_gpu_id = idx + _gpu_name = name + _gpu_mem_gb = mem_mi // 1024 + HAS_GPU = True + _xp = _real_cupy + _cp = _real_cupy + _cp_ndimage = _real_cupy_ndimage + return + except Exception as e: + logger.warning(f"GPU indisponible (CUDA_VISIBLE_DEVICES={cuda_visible}): {e}") + _xp = np + _cp = None + _cp_ndimage = None + HAS_GPU = False + return + + # No CUDA_VISIBLE_DEVICES set — test each GPU in subprocess import subprocess _working_gpu = None - for idx, name, cap_str, mem_mi, score, major in _candidate_gpus: + for idx, name, cap_str, mem_mi, score, major in candidates: result = subprocess.run( ['python3', '-c', - 'import cupy; a=cupy.array([1.0,2.0],dtype=cupy.float32); print(cupy.sum(a).get())'], + 'import cupy; a=cupy.array([1.0,2.0],dtype=cupy.float32); ' + 'dev=cupy.cuda.runtime.getDevice(); ' + 'print(f"OK:{dev}:{cupy.sum(a).get()}")'], capture_output=True, text=True, timeout=120, env={**os.environ, 'CUDA_VISIBLE_DEVICES': str(idx)}, ) - if result.returncode == 0 and '3.0' in result.stdout: - _working_gpu = (idx, name, mem_mi) - break + stdout = result.stdout.strip() + if result.returncode == 0 and stdout.startswith('OK:') and '3.0' in stdout: + # Verify the device actually used is device 0 (the GPU we targeted) + parts = stdout.split(':') + if len(parts) >= 2 and parts[1] == '0': + _working_gpu = (idx, name, mem_mi) + break + else: + logger.warning(f"GPU {idx} ({name}, sm_{cap_str}) faux positif CPU fallback") logger.warning(f"GPU {idx} ({name}, sm_{cap_str}) non compatible: " f"{result.stderr.strip().splitlines()[-1] if result.stderr else 'inconnue'}") @@ -141,27 +210,38 @@ def _init_gpu(): # --------------------------------------------------------------------------- def num_gpus(): - """Return 1 if GPU is active, 0 otherwise.""" - return 1 if HAS_GPU else 0 + """Return the number of available GPUs (after restrict_gpus filtering).""" + return len(_gpu_candidates) + + +def available_gpu_ids(): + """Return list of host-level GPU indices available for processing. + + Respects any prior restrict_gpus() call. + """ + return [c['id'] for c in _gpu_candidates] def restrict_gpus(gpu_ids: list[int], set_env_var: bool = False): - """No-op — GPU is auto-selected at import time.""" - pass + """Restrict GPU selection to specific host-level indices. + + Stores the restriction to be applied during _init_gpu(). + """ + global _restricted_gpu_ids + _restricted_gpu_ids = gpu_ids def set_active_gpu(gpu_id): - """No-op — GPU is auto-selected at import time.""" - pass + """Restrict to a single GPU by host-level index.""" + global _restricted_gpu_ids + _restricted_gpu_ids = [gpu_id] def _gpu_available(): """Check if GPU is usable right now.""" - if not HAS_GPU: - return False try: _init_gpu() - return _cp is not None + return HAS_GPU and _cp is not None except Exception: return False @@ -212,7 +292,7 @@ def to_gpu(arr): try: return _cp.asarray(arr.astype(np.float32)) except Exception: - pass + disable_gpu() return arr.astype(np.float32) @@ -302,7 +382,7 @@ def safe_gpu_call(func, *args, **kwargs): return func(*args, **kwargs) except Exception as e: err_msg = str(e) - if _cp is not None and ('CUDA' in err_msg or 'cuda' in err_msg or 'GPU' in err_msg): + if _cp is not None and ('CUDA' in err_msg or 'cuda' in err_msg or 'GPU' in err_msg or 'Out of memory' in err_msg): logger.warning(f"Erreur GPU ({e.__class__.__name__}), retry en CPU...") disable_gpu() return func(*args, **kwargs) diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py index 7231e80..fe343f8 100644 --- a/lidar_pipeline/pipeline.py +++ b/lidar_pipeline/pipeline.py @@ -57,14 +57,15 @@ _file_filter = FilePrefixFilter() from .dtm import classify_ground, create_dtm_fast from .visualizations import ( SharedDEM, - generate_hillshade, generate_slope, + generate_hillshade, generate_slope, generate_aspect, generate_openness, generate_mslrm, generate_sailore, generate_roughness, generate_wavelet, - generate_svf, generate_aniso_open, - generate_flow_accumulation, - generate_anomaly_mask, - ) + generate_solar, + generate_svf, generate_aniso_open, + generate_flow_accumulation, + generate_anomaly_mask, +) from .gpu import gpu_cleanup, num_gpus, restrict_gpus, safe_gpu_call from .ign import generate_ign_overlay from .rendering import tif_to_png @@ -76,6 +77,7 @@ from .rendering import tif_to_png VIZ_STEPS = [ ('hillshade', generate_hillshade), ('slope', generate_slope), + ('aspect', generate_aspect), ('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)), @@ -85,6 +87,7 @@ VIZ_STEPS = [ ('roughness', generate_roughness), ('wavelet', generate_wavelet), ('flow_acc', generate_flow_accumulation), + ('solar', generate_solar), ('anomaly', generate_anomaly_mask), ('ortho', lambda d, b, v, r: generate_ign_overlay( d, b, v, r, @@ -463,11 +466,12 @@ class LidarArchaeoPipeline: logger.info(f"Fichiers: {len(files)}") with ProcessPoolExecutor(max_workers=self.workers) as executor: - # Pass resolutions as comma-separated string for multiprocessing serialization + # Round-robin assign each file to a real GPU host index + active_ids = self.gpu_ids if self.gpu_ids else _gpu_mod.available_gpu_ids() resolutions_str = ','.join(str(r) for r in self.resolutions) future_to_file = { - executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), resolutions_str, self.force, self.ground_method, self.force_classify, self.keep_tif, self.quality, self.only_viz, self.skip_viz, self.output_format, gpu_id % n_gpus, gpu_ids=self.gpu_ids): laz_file - for gpu_id, laz_file in enumerate(files) + executor.submit(_process_file_standalone, str(laz_file), str(self.input_dir), str(self.output_dir), resolutions_str, self.force, self.ground_method, self.force_classify, self.keep_tif, self.quality, self.only_viz, self.skip_viz, self.output_format, active_ids[file_idx % len(active_ids)] if active_ids else None): laz_file + for file_idx, laz_file in enumerate(files) } done = 0 try: @@ -535,18 +539,13 @@ class LidarArchaeoPipeline: logger.warning(f" Note: Impossible de supprimer les fichiers temporaires: {e}") -def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_id=None, gpu_ids=None): +def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, force=False, ground_method='auto', force_classify=False, keep_tif=False, quality=98, only_viz=None, skip_viz=None, output_format='avif', gpu_id=None): """Standalone function for multiprocessing — creates its own pipeline instance. Each worker gets its own temp directory to avoid file conflicts. When multiple GPUs are available, each worker is assigned a GPU via CUDA_VISIBLE_DEVICES to balance load across GPUs. """ - # Restrict visible GPUs first, then pick one for this worker - if gpu_ids is not None: - from .gpu import restrict_gpus - restrict_gpus(gpu_ids) - if gpu_id is not None and gpu_id >= 0: from .gpu import set_active_gpu set_active_gpu(gpu_id) diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py index c5e4509..962878a 100644 --- a/lidar_pipeline/rendering.py +++ b/lidar_pipeline/rendering.py @@ -144,6 +144,14 @@ 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)', @@ -154,11 +162,10 @@ COLORMAPS = { }, 'wavelet': { 'cmap': 'cividis', - 'title': 'Ondelette Mexican Hat + Gabor directionnelle (multi-échelle)', - 'legend': 'Réponse RMS 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': 'fixed', 'vmin_val': 0, - 'vmax_mode': 'percentile', 'vmax_pct': 98, + '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', + 'vmin_mode': 'symmetric', 'sym_pct': (2, 98), }, 'flow_acc': { 'cmap': 'YlGn', @@ -176,6 +183,14 @@ COLORMAPS = { 'vmin_mode': 'fixed', 'vmin_val': 0, 'vmax_mode': 'fixed', 'vmax_val': 1, }, + 'solar': { + 'cmap': 'gray', + 'title': 'Éclairage Solaire', + 'legend': "Illumination solaire (azimut 90°, altitude 30°)\nClair = Face éclairée | Sombre = Zone d'ombre", + 'description': "Simulation de l'éclairage solaire matinal", + 'vmin_mode': 'fixed', 'vmin_val': 0, + 'vmax_mode': 'fixed', 'vmax_val': 1, + }, } # RGB entries (ortho/topo) are handled specially @@ -822,7 +837,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', - 'flow_acc', 'slope', 'roughness', 'wavelet'] + 'flow_acc', 'solar', 'slope', 'roughness', 'wavelet'] def sort_key(f): name = f.stem.lower() diff --git a/lidar_pipeline/tests/test_pipeline.py b/lidar_pipeline/tests/test_pipeline.py index bce852f..06deb2d 100644 --- a/lidar_pipeline/tests/test_pipeline.py +++ b/lidar_pipeline/tests/test_pipeline.py @@ -19,16 +19,10 @@ class TestVizSteps: names = [name for name, _ in VIZ_STEPS] assert len(names) == len(set(names)), "VIZ_STEPS has duplicate names" - def test_no_solar_in_viz_steps(self): - """Solar visualization was removed.""" - from lidar_pipeline.pipeline import VIZ_STEPS - names = [name for name, _ in VIZ_STEPS] - assert "solar" not in names - def test_expected_visualization_count(self): - """Should have 13 visualizations (11 terrain + ortho + topo).""" + """Should have 16 visualizations (14 terrain + ortho + topo).""" from lidar_pipeline.pipeline import VIZ_STEPS - assert len(VIZ_STEPS) == 13 + assert len(VIZ_STEPS) == 16 def test_ortho_and_topo_present(self): from lidar_pipeline.pipeline import VIZ_STEPS diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 041dffd..d6416ce 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -513,6 +513,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. + + 0° = North, 90° = East, 180° = South, 270° = West. + Direction toward which the terrain descends. + """ + gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" + logger.info(f" → Aspect (Orientation des pentes){gpu_tag}...") + t0 = time.time() + output = vis_dir / f"{basename}_aspect.tif" + + try: + if shared: + transform = shared.transform + crs = shared.crs + dy = shared.dy + dx = shared.dx + nan_mask = shared.nan_mask + if _gpu_mod.HAS_GPU: + dy = to_gpu(dy) + dx = to_gpu(dx) + 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) + aspect = xp.arctan2(dy, dx) * 180 / xp.pi + aspect = xp.mod(aspect, 360) + _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]' if _gpu_mod.HAS_GPU else ''}") + return output + except Exception as e: + logger.error(f" ✗ Erreur aspect: {e}", exc_info=True) + return None + + # ============================================================ # GPU-accelerated visualizations # ============================================================ @@ -829,39 +865,76 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None): # ============================================================ -# Wavelet (Mexican Hat + Directional Gabor) +# Exposition des surfaces (Éclairage Solaire) # ============================================================ -def _gabor_kernel_2d(size, sigma, wavelength, theta): - """Create a 2D Gabor kernel. +def generate_solar(dem_file, basename, vis_dir, resolution, shared=None): + """Generate solar irradiance simulation. - Args: - size: kernel size (odd integer) - sigma: standard deviation - wavelength: wavelength of sinusoid - theta: orientation angle in radians (0 = horizontal) + Simulates morning sunlight (azimuth 90°, altitude 30°) to reveal + subtle topographic features through shadow effects. """ - 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 + logger.info(" → Exposition des surfaces (Éclairage Solaire)...") + t0 = time.time() + output = vis_dir / f"{basename}_solar.tif" + try: + if shared: + transform = shared.transform + crs = shared.crs + nan_mask = shared.nan_mask + dx = shared.dx + dy = shared.dy + else: + dem_np, transform, crs = _read_dem(dem_file) + nan_mask = np.isnan(dem_np) + dem_filled, _ = _fill_nans(dem_np) + dy, dx = np.gradient(dem_filled, resolution, resolution) + + # Solar parameters: morning sun (azimuth 90° = east, altitude 30°) + sun_azimuth = np.radians(90) + sun_altitude = np.radians(30) + + # Aspect from gradient + aspect = np.degrees(np.arctan2(-dx, -dy)) + aspect[aspect < 0] += 360 + + # Slope in radians + slope_rad = np.arctan(np.sqrt(dx**2 + dy**2)) + + # Solar irradiance calculation + irradiance = (np.sin(sun_altitude) * np.sin(slope_rad) + + np.cos(sun_altitude) * np.cos(slope_rad) * + np.cos(np.radians(aspect) - sun_azimuth)) + + # Clip to valid range [0, 1] + irradiance = np.clip(irradiance, 0, 1) + irradiance[nan_mask] = np.nan + + _save_tif(output, irradiance.astype(np.float32), transform, crs) + logger.info(f" ✓ Exposition des surfaces terminée ({time.time()-t0:.1f}s)") + return output + except Exception as e: + logger.error(f" ✗ Erreur exposition: {e}", exc_info=True) + return None + + +# ============================================================ +# Wavelet (Mexican Hat) +# ============================================================ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): - """Multi-scale wavelet analysis: Mexican Hat + Directional Gabor (GPU if available). + """Mexican Hat wavelet multi-scale analysis (GPU if available). - 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. + 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 + + Uses std normalization per scale and weighted RMS combination + with emphasis on archaeologically relevant scales (2-50m). """ gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Ondelette Mexican Hat + Gabor directionnelle{gpu_tag}...") + logger.info(f" → Ondelette Mexican Hat multi-échelle{gpu_tag}...") t0 = time.time() output = vis_dir / f"{basename}_wavelet.tif" @@ -877,25 +950,23 @@ 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)) - # --- Mexican Hat scales --- min_scale = max(resolution * 2, 1.0) candidate_scales = [0.5, 1, 2, 5, 10, 20, 50, 100] - mex_scales = [s for s in candidate_scales if s >= min_scale] + scales = [s for s in candidate_scales if s >= min_scale] - mex_weights_map = { + scale_weights = { 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, } - mex_weights = np.array([mex_weights_map.get(s, 1.0) for s in mex_scales]) + weights = np.array([scale_weights.get(s, 1.0) for s in scales]) - logger.info(f" Échelles CWT: {mex_scales}m (résolution {resolution}m/px)") + logger.info(f" Échelles CWT: {scales}m (résolution {resolution}m/px)") - from scipy.ndimage import gaussian_laplace, convolve + from scipy.ndimage import gaussian_laplace wavelet_stack = [] - # Mexican Hat (radial) — multi-scale - for scale_m in mex_scales: + for scale_m in scales: sigma_px = scale_m / resolution if _gpu_mod.HAS_GPU: try: @@ -912,44 +983,12 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): response = response / std_val wavelet_stack.append(response) - # 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 = all_weights[:, np.newaxis, np.newaxis] + weights_3d = 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(all_weights)) + combined = np.sqrt(np.nansum((stack ** 2) * weights_3d, axis=0) / np.sum(weights)) combined[nan_mask] = np.nan _save_tif(output, combined.astype(np.float32), transform, crs) @@ -960,21 +999,115 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None): return None +# ============================================================ +# Flow Accumulation helpers (module-level for numba caching) +# ============================================================ + +def _priority_flood(dem, nodata_mask): + """Priority-flood algorithm for sink filling (Wang & Liu 2006). + + O(n log n) compared to 50 iterations of minimum_filter. + Fills pits so water can flow downhill. NaN cells are treated as + closed (never pushed into the heap). + """ + import heapq + + rows, cols = dem.shape + filled = dem.copy() + closed = nodata_mask.copy() + open_queue = [] + + # Initialize border cells (skip NaN cells) + for r in range(rows): + for c in [0, cols - 1]: + if not closed[r, c]: + heapq.heappush(open_queue, (filled[r, c], r, c)) + closed[r, c] = True + for c in range(1, cols - 1): + for r in [0, rows - 1]: + if not closed[r, c]: + heapq.heappush(open_queue, (filled[r, c], r, c)) + closed[r, c] = True + + dx8 = [1, 1, 0, -1, -1, -1, 0, 1] + dy8 = [0, 1, 1, 1, 0, -1, -1, -1] + + while open_queue: + elev, r, c = heapq.heappop(open_queue) + for d in range(8): + nr, nc = r + dy8[d], c + dx8[d] + if 0 <= nr < rows and 0 <= nc < cols and not closed[nr, nc]: + if filled[nr, nc] < elev: + filled[nr, nc] = elev # Fill the pit + closed[nr, nc] = True + heapq.heappush(open_queue, (filled[nr, nc], nr, nc)) + + return filled + + +def _d8_accumulate_numba(dem_filled, flow_dir, nodata_mask, rows, cols): + """JIT-compiled D8 flow accumulation (top-down via elevation sort). + + Uses numba for ~100x speedup over pure Python loop. + Falls back to pure Python if numba is unavailable. + """ + try: + from numba import njit + + @njit(cache=True) + def _accumulate(dem, fdir, nodata, rows, cols): + dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int8) + dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int8) + + flow_acc = np.ones((rows, cols), dtype=np.float32) + for r in range(rows): + for c in range(cols): + if nodata[r, c]: + flow_acc[r, c] = 0.0 + + # Sort cells by elevation descending (highest first) + n_cells = rows * cols + flat_dem = dem.ravel() + sort_idx = np.argsort(-flat_dem) + + # Accumulate top-down (highest cell first) + for i in range(n_cells): + cell = sort_idx[i] + r = cell // cols + c = cell % cols + if nodata[r, c]: + continue + d = fdir[r, c] + if d < 0: + continue + nr = r + dy8[d] + nc = c + dx8[d] + if 0 <= nr < rows and 0 <= nc < cols and not nodata[nr, nc]: + flow_acc[nr, nc] += flow_acc[r, c] + + return flow_acc + + return _accumulate(dem_filled, flow_dir, nodata_mask, rows, cols) + + except ImportError: + return None + + # ============================================================ # Flow Accumulation # ============================================================ def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=None): - """Flow Accumulation — priority-flood sink filling + D8 accumulation (GPU if available). + """Flow Accumulation — priority-flood sink filling + D8 accumulation. 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. - Uses log10 transformation for visualization. + D8 direction is computed via vectorized numpy slicing. + Accumulation uses numba JIT (cached at module level) or pure Python fallback. """ - gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else "" - logger.info(f" → Accumulation d'écoulement (flow accumulation){gpu_tag}...") + logger.info(f" → Accumulation d'écoulement (flow accumulation)...") t0 = time.time() output = vis_dir / f"{basename}_flow_acc.tif" @@ -990,92 +1123,70 @@ def generate_flow_accumulation(dem_file, basename, vis_dir, resolution, shared=N nan_mask = np.isnan(dem_np) filled, _ = _fill_nans(dem_np) - # Priority-flood sink filling (Wang & Liu 2006, O(n log n)) - from heapq import heappush, heappop - rows, cols = filled.shape - dem_filled = filled.copy() + rows, cols = dem_np.shape - # 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 - - 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 + # Sink filling — priority-flood (O(n log n), NaN-aware) + dem_filled = _priority_flood(filled, nan_mask) logger.info(f" ✓ Sink filling terminé ({time.time()-t0:.1f}s)") - # D8 flow direction and accumulation - flow_acc = np.zeros((rows, cols), dtype=np.int64) - flow_dir = np.zeros((rows, cols), dtype=np.int8) - 1 + # D8 flow direction — vectorized via numpy slicing + dx8 = np.array([1, 1, 0, -1, -1, -1, 0, 1], dtype=np.int32) + dy8 = np.array([0, 1, 1, 1, 0, -1, -1, -1], dtype=np.int32) + dist8 = np.array([1.0, np.sqrt(2), 1.0, np.sqrt(2), + 1.0, np.sqrt(2), 1.0, np.sqrt(2)]) - # D8 neighbors (ordered by angle) - neighbors = [ - (-1, 0), (-1, 1), ( 0, 1), ( 1, 1), - ( 1, 0), ( 1, -1), ( 0, -1), (-1, -1) - ] + flow_dir = np.full((rows, cols), -1, dtype=np.int8) + max_slope = np.zeros((rows, cols), dtype=np.float64) - # Process cells in ascending elevation order for correct accumulation - flat_idx = np.argsort(dem_filled.ravel()) - flow_acc_flat = flow_acc.ravel() + padded = np.pad(dem_filled, 1, mode='constant', + constant_values=np.nanmax(dem_filled[~np.isnan(dem_filled)]) + 10000) - for idx in flat_idx: - y = idx // cols - x = idx % cols + for d in range(8): + nx = 1 + dx8[d] + ny = 1 + dy8[d] + neighbor_elev = padded[ny:ny + rows, nx:nx + cols] + slope = (dem_filled - neighbor_elev) / (dist8[d] * resolution) + slope[nan_mask] = -1 + better = slope > max_slope + flow_dir[better] = d + max_slope[better] = slope[better] - # 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 + logger.info(f" ✓ Direction D8 terminée ({time.time()-t0:.1f}s)") - flow_dir[y, x] = best_dir - flow_acc[y, x] = 1 # Count self + # D8 accumulation — numba JIT (module-level cache) or pure Python fallback + result = _d8_accumulate_numba(dem_filled, flow_dir, nan_mask.astype(np.bool_), rows, cols) - # 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] + if result is not None: + flow_acc = result + logger.info(f" Accumulation D8 via numba") + else: + # Pure Python fallback + logger.info(f" Accumulation D8 via Python (installez numba pour accélérer)") + flat_dem = dem_filled[~nan_mask].flatten() + valid_indices = np.where(~nan_mask.flatten())[0] + sort_order = valid_indices[np.argsort(-flat_dem)] + + flow_acc = np.ones((rows, cols), dtype=np.float32) + flow_acc[nan_mask] = 0 + + for idx in sort_order: + r, c = divmod(idx, cols) + d = flow_dir[r, c] + if d < 0: + continue + nr, nc = r + dy8[d], c + dx8[d] + if 0 <= nr < rows and 0 <= nc < cols and not nan_mask[nr, nc]: + flow_acc[nr, nc] += flow_acc[r, c] 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)) + # Log transform + flow_result = np.log1p(flow_acc) 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]' if _gpu_mod.HAS_GPU else ''}") + logger.info(f" ✓ Flow accumulation terminé ({time.time()-t0:.1f}s)") return output except Exception as e: logger.error(f" ✗ Erreur flow accumulation: {e}", exc_info=True) @@ -1185,16 +1296,18 @@ def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, rows, cols = nan_mask.shape # Collect available visualization layers from disk - # Each is loaded, normalized to z-score, and contributes to the composite + # Each is loaded, converted to |z-score|, and contributes to a weighted sum. + # The weighted sum of |z| acts as a "vote": pixels where multiple layers + # show anomalies get higher scores than pixels where only one layer fires. layer_configs = [ # (filename_pattern, weight) - ("mslrm", 2.5), # Multi-scale relief — strongest signal - ("negative_openness", 1.8), # Fossés, dolines - ("roughness", 1.5), # Surface irregularity - ("wavelet", 1.5), # Circular + linear structures - ("svf", 1.3), # Sky-view depressions - ("flow_acc", 1.2), # Drainage channels / ditches - ("aniso_open", 1.0), # Anisotropic structures + ("mslrm", 2.5), # Multi-scale relief — strongest signal + ("negative_openness", 2.0), # Fossés, dolines + ("roughness", 1.8), # Surface irregularity + ("wavelet", 1.5), # Circular + linear structures + ("svf", 1.3), # Sky-view depressions + ("flow_acc", 1.2), # Drainage channels / ditches + ("aniso_open", 1.0), # Anisotropic structures ("positive_openness", 0.8), # Surélevations ] @@ -1208,13 +1321,13 @@ def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, try: with rasterio.open(layer_path) as src: data = src.read(1).astype(np.float64) - # Z-score normalization + # Absolute z-score (captures both positive and negative deviations) valid = data[~nan_mask] if len(valid) == 0: continue mean_val = np.nanmean(valid) std_val = max(np.nanstd(valid), 0.01) - zscore = (data - mean_val) / std_val + zscore = np.abs(data - mean_val) / std_val zscore[nan_mask] = 0.0 layers.append((zscore, weight)) except Exception as e: @@ -1230,7 +1343,7 @@ def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, sigma_px = max(5, 10.0 / resolution) local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px) if shared else \ _filter_nanaware(dem_np_safe, xp_gaussian_filter, sigma=sigma_px) - quick_mslrm = dem_np_safe - local_mean + quick_mslrm = np.abs(dem_np_safe - local_mean) quick_mslrm[nan_mask] = np.nan valid = quick_mslrm[~nan_mask] std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01 @@ -1254,25 +1367,29 @@ def generate_anomaly_mask(dem_file, basename, vis_dir, resolution, shared=None, roughness = roughness / std_val_r layers.append((roughness, 1.5)) - # Weighted RMS combination (unsigned — all deviations are suspicious) + # Weighted SUM of |z-scores| (not RMS — each layer votes independently) combined = np.zeros((rows, cols), dtype=np.float64) total_weight = 0.0 for layer_data, weight in layers: - combined += (layer_data ** 2) * weight + combined += layer_data * weight total_weight += weight if total_weight > 0: - combined = np.sqrt(combined / total_weight) + combined = combined / total_weight - # Apply n_sigma threshold: pixels below n_sigma are suppressed - combined = np.maximum(combined - n_sigma, 0.0) + # Adaptive threshold: suppress pixels below the (100 - n_sigma*10)th percentile. + # With n_sigma=2.0 → 80th percentile: keep only the top 20% of signal. + # This adapts to each tile's terrain instead of a fixed z-score cutoff. + threshold_pct = max(50, min(95, 100 - n_sigma * 10)) + threshold_val = np.percentile(combined[~nan_mask], threshold_pct) + combined = np.clip(combined - threshold_val, 0, None) - # Rescale to 0–1 for visualization (percentile-based stretch) - valid_combined = combined[~nan_mask] - if len(valid_combined) > 0 and np.nanmax(valid_combined) > 0: - p99 = np.percentile(valid_combined, 99) - if p99 > 0: - combined = np.clip(combined / p99, 0, 1) + # Rescale survivors to 0–1 + above_thresh = combined[combined > 0] + if len(above_thresh) > 0: + p95 = np.percentile(above_thresh, 95) + if p95 > 0: + combined = np.clip(combined / p95, 0, 1) combined[nan_mask] = np.nan