diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py index 406b601..e8863ef 100644 --- a/lidar_pipeline/visualizations.py +++ b/lidar_pipeline/visualizations.py @@ -58,6 +58,59 @@ def _read_dem(dem_file): return src.read(1), src.transform, src.crs +def _fill_nans(arr): + """Fill NaN values using nearest-neighbor interpolation. + + Returns (filled_array, nan_mask) so the caller can restore NaN after filtering. + """ + from scipy.interpolate import NearestNDInterpolator + nan_mask = np.isnan(arr) + if not np.any(nan_mask): + return arr, nan_mask + valid = ~nan_mask + y_coords, x_coords = np.where(valid) + if len(y_coords) == 0: + return np.zeros_like(arr), nan_mask + z_values = arr[valid] + interp = NearestNDInterpolator( + np.column_stack((y_coords, x_coords)), z_values + ) + y_missing, x_missing = np.where(nan_mask) + filled = arr.copy() + filled[y_missing, x_missing] = interp(y_missing, x_missing) + return filled, nan_mask + + +def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs): + """Apply a filter to an array while preserving NaN zones. + + 1. Fill NaN with nearest-neighbor interpolation + 2. Apply the filter + 3. Restore original NaN mask on the result + + Args: + arr: Input array (numpy or cupy). + filter_func: Function that takes (array, *args, **kwargs) and returns filtered array. + use_gpu: If True, apply filter on GPU (send filled array to GPU first). + Returns: + Filtered array with original NaN positions preserved. + """ + is_gpu_arr = HAS_GPU and isinstance(arr, cp.ndarray) + arr_np = to_cpu(arr) if is_gpu_arr else arr + + filled, nan_mask = _fill_nans(arr_np) + + if use_gpu and HAS_GPU: + filled_gpu = to_gpu(filled) + result_gpu = filter_func(filled_gpu, *args, **kwargs) + result = to_cpu(result_gpu) + else: + result = filter_func(filled, *args, **kwargs) + + result[nan_mask] = np.nan + return result + + # ============================================================ # Core terrain visualizations # ============================================================ @@ -179,11 +232,13 @@ def generate_lrm(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) - local_mean = xp_gaussian_filter(dem, sigma=15/resolution) - lrm = dem - local_mean - lrm_np = to_cpu(lrm).astype(np.float32) - _save_tif(output, lrm_np, transform, crs) + nan_mask = np.isnan(dem_np) + local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15/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 logger.info(f" ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -320,21 +375,25 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) + nan_mask = np.isnan(dem_np) sigmas = [5, 10, 25, 50, 100] lrm_stack = [] for sigma in sigmas: sigma_px = sigma / resolution - local_mean = xp_gaussian_filter(dem, sigma=sigma_px) - lrm = dem - local_mean - lrm_norm = lrm / max(float(xp.nanstd(lrm)), 0.01) - lrm_stack.append(lrm_norm) + local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px) + lrm = dem_np - local_mean + lrm[nan_mask] = np.nan + valid_lrm = lrm[~nan_mask] + lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 + lrm_norm = lrm / lrm_std + lrm_stack.append(lrm_norm.astype(np.float32)) - mslrm = xp.sqrt(xp.mean(xp.array(lrm_stack) ** 2, axis=0)) - mslrm_np = to_cpu(mslrm).astype(np.float32) - _save_tif(output, mslrm_np, transform, crs) + lrm_array = np.array(lrm_stack) + mslrm = np.sqrt(np.nanmean(lrm_array ** 2, axis=0)) + mslrm[nan_mask] = np.nan + _save_tif(output, mslrm.astype(np.float32), transform, crs) logger.info(f" ✓ MSRM terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -355,24 +414,28 @@ def generate_tpi(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) + nan_mask = np.isnan(dem_np) fine_size = int(5 / resolution) if fine_size % 2 == 0: fine_size += 1 - tpi_fine = dem - xp_uniform_filter(dem, size=fine_size) + fine_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=fine_size) + tpi_fine = dem_np - fine_mean + tpi_fine[nan_mask] = np.nan broad_size = int(100 / resolution) if broad_size % 2 == 0: broad_size += 1 - tpi_broad = dem - xp_uniform_filter(dem, size=broad_size) + broad_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=broad_size) + tpi_broad = dem_np - broad_mean + tpi_broad[nan_mask] = np.nan - fine_std = max(float(xp.nanstd(tpi_fine)), 0.01) - broad_std = max(float(xp.nanstd(tpi_broad)), 0.01) + fine_std = max(np.nanstd(tpi_fine[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01 + broad_std = max(np.nanstd(tpi_broad[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01 tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std) - tpi_np = to_cpu(tpi_combined).astype(np.float32) + tpi_combined[nan_mask] = np.nan - _save_tif(output, tpi_np, transform, crs) + _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: @@ -448,31 +511,34 @@ def generate_sailore(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) + nan_mask = np.isnan(dem_np) - gy, gx = xp.gradient(dem, resolution) - slope = xp.arctan(xp.sqrt(gx**2 + gy**2)) - slope_deg = xp.degrees(slope) + gy, gx = np.gradient(dem_np, resolution) + slope = np.arctan(np.sqrt(gx**2 + gy**2)) + slope_deg = np.degrees(slope) + slope_deg[nan_mask] = np.nan sigma_min = 2.0 / resolution sigma_max = 25.0 / resolution - slope_norm = xp.clip(slope_deg / 30.0, 0, 1) - adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min) + slope_norm = np.clip(slope_deg / 30.0, 0, 1) - lrm_fine = dem - xp_gaussian_filter(dem, sigma=sigma_min) - lrm_medium = dem - xp_gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2) - lrm_coarse = dem - xp_gaussian_filter(dem, sigma=sigma_max) + lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min) + lrm_fine[nan_mask] = np.nan + lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2) + lrm_medium[nan_mask] = np.nan + lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max) + lrm_coarse[nan_mask] = np.nan w_fine = slope_norm - w_medium = 1 - 2 * xp.abs(slope_norm - 0.5) + w_medium = 1 - 2 * np.abs(slope_norm - 0.5) w_coarse = 1 - slope_norm w_total = w_fine + w_medium + w_coarse w_total[w_total == 0] = 1 sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total - sailore_np = to_cpu(sailore).astype(np.float32) + sailore[nan_mask] = np.nan - _save_tif(output, sailore_np, transform, crs) + _save_tif(output, sailore.astype(np.float32), transform, crs) logger.info(f" ✓ SAILORE terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -493,15 +559,16 @@ def generate_roughness(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np.astype(np.float64)) + nan_mask = np.isnan(dem_np) window_size = int(5 / resolution) if window_size % 2 == 0: window_size += 1 - # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (GPU-accelerated) - local_mean = xp_uniform_filter(dem, size=window_size) - local_mean_sq = xp_uniform_filter(dem * dem, size=window_size) - roughness = xp.sqrt(local_mean_sq - local_mean * local_mean) + # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (NaN-aware) + local_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=window_size) + local_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=window_size) + roughness = np.sqrt(np.maximum(local_mean_sq - local_mean * local_mean, 0)) + roughness[nan_mask] = np.nan roughness = to_cpu(roughness) _save_tif(output, roughness, transform, crs) @@ -525,24 +592,28 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) + nan_mask = np.isnan(dem_np) - lrm = dem - xp_gaussian_filter(dem, sigma=15 / resolution) - lrm_mean = xp.nanmean(lrm) - lrm_std = max(float(xp.nanstd(lrm)), 0.01) + lrm_mean_val = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15 / resolution) + lrm = dem_np - lrm_mean_val + lrm[nan_mask] = np.nan + valid_lrm = lrm[~nan_mask] + lrm_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0 + lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01 z_score = (lrm - lrm_mean) / lrm_std window = int(10 / resolution) if window % 2 == 0: window += 1 - local_mean = xp_uniform_filter(z_score, size=window) - z_mean = xp.nanmean(z_score) - z_std = max(float(xp.nanstd(z_score)), 0.01) + local_mean = _filter_nanaware(z_score, xp_uniform_filter, size=window) + z_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0 + z_std = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01 morans_i = z_score * (local_mean - z_mean) / z_std - anomaly_score = xp.abs(z_score) * xp.sign(morans_i) + anomaly_score = np.abs(z_score) * np.sign(morans_i) + anomaly_score[nan_mask] = np.nan - _save_tif(output, to_cpu(anomaly_score), transform, crs) + _save_tif(output, anomaly_score.astype(np.float32), transform, crs) logger.info(f" ✓ Anomalies terminé ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: @@ -566,26 +637,31 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution): try: dem_np, transform, crs = _read_dem(dem_file) - dem = to_gpu(dem_np) + nan_mask = np.isnan(dem_np) scales = [2, 5, 10, 20, 50] wavelet_stack = [] for scale_m in scales: sigma_px = scale_m / resolution + filled, _ = _fill_nans(dem_np.astype(np.float64)) if HAS_GPU: from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace - response = -gpu_gaussian_laplace(dem, sigma=sigma_px) + response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px) + response = to_cpu(response) else: from scipy.ndimage import gaussian_laplace - response = to_gpu(-gaussian_laplace(to_cpu(dem).astype(np.float64), sigma=sigma_px)) - response /= max(float(xp.nanstd(response)), 0.01) + response = -gaussian_laplace(filled, sigma=sigma_px) + 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) - combined = xp.sqrt(xp.mean(xp.array(wavelet_stack) ** 2, axis=0)) - combined_np = to_cpu(combined).astype(np.float32) + combined = np.sqrt(np.nanmean(np.array(wavelet_stack) ** 2, axis=0)) + combined[nan_mask] = np.nan - _save_tif(output, combined_np, transform, crs) + _save_tif(output, combined.astype(np.float32), transform, crs) logger.info(f" ✓ Ondelette terminée ({time.time()-t0:.1f}s){gpu_tag}") return output except Exception as e: