From 54800cb51673fe73dc7b582e5788b2ab291077ac Mon Sep 17 00:00:00 2001 From: Jacquin Antoine Date: Sat, 9 May 2026 22:22:28 +0200 Subject: [PATCH] =?UTF-8?q?Pipeline=20LiDAR:=20acc=C3=A9l=C3=A9ration=20GP?= =?UTF-8?q?U=20(CuPy),=20sortie=20WebP,=20script=20run.sh?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Accélération GPU via CuPy pour SVF, Openness, LRM, MSRM, SAILORE, TPI, wavelet - Fallback automatique vers numpy si GPU non disponible - Sortie WebP sans perte (remplace PNG, fichiers plus petits) - Script run.sh avec options -g (GPU), -w (workers), -r (résolution) - Docker basé sur nvidia/cuda:12.4.0-devel pour support CuPy - Docker tourne en uid/gid 1000:1000 - Légendes explicites différenciant LRM vs MSRM vs SAILORE - Correction bug ordre elif (mslrm avant lrm) - Retrait de geomorphons et VAT (demande utilisateur) Co-Authored-By: Claude Opus 4.6 --- Dockerfile | 9 +- README.md | 13 +- process_lidar.py | 326 +++++++++++++++++++++++++++-------------------- run.sh | 72 +++++++++++ 4 files changed, 278 insertions(+), 142 deletions(-) create mode 100755 run.sh diff --git a/Dockerfile b/Dockerfile index 2406fbc..6ed79fc 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,9 +1,9 @@ -FROM ubuntu:22.04 +FROM nvidia/cuda:12.4.0-devel-ubuntu22.04 ENV DEBIAN_FRONTEND=noninteractive ENV TZ=Europe/Paris -# Install PDAL and Python from Ubuntu packages +# Install PDAL and system packages RUN apt-get update && apt-get install -y --no-install-recommends \ pdal \ gdal-bin \ @@ -29,6 +29,9 @@ RUN pip3 install --no-cache-dir \ scipy \ tqdm +# Install CuPy for GPU acceleration (optional - will fallback to numpy if not available) +RUN pip3 install --no-cache-dir cupy-cuda12x || echo "CuPy not available - GPU acceleration disabled" + # Copy scripts COPY process_lidar.py /usr/local/bin/ RUN chmod +x /usr/local/bin/process_lidar.py @@ -43,4 +46,4 @@ USER lidar VOLUME ["/data"] -CMD ["python3", "/usr/local/bin/process_lidar.py", "/data/input", "-o", "/data/output"] +CMD ["python3", "/usr/local/bin/process_lidar.py", "/data/input", "-o", "/data/output"] \ No newline at end of file diff --git a/README.md b/README.md index 27622b4..6dd1db2 100644 --- a/README.md +++ b/README.md @@ -46,13 +46,22 @@ mkdir -p input # Copiez vos fichiers .laz dans input/ cp /chemin/vos/fichiers/*.laz input/ -# Build l'image Docker +# Build l'image Docker (avec support GPU NVIDIA) docker build -t lidar-archeo . ``` ## Utilisation -### Traitement standard (1 CPU) +### Traitement avec accélération GPU (recommandé) +```bash +# Nécessite une carte NVIDIA + nvidia-container-toolkit +docker run --rm --gpus all \ + -v $(pwd)/input:/data/input:ro \ + -v $(pwd)/output:/data/output \ + lidar-archeo +``` + +### Traitement standard (CPU seul, sans GPU) ```bash docker run --rm \ -v $(pwd)/input:/data/input:ro \ diff --git a/process_lidar.py b/process_lidar.py index cbb59a6..8da10f7 100755 --- a/process_lidar.py +++ b/process_lidar.py @@ -1,7 +1,8 @@ #!/usr/bin/env python3 """ -Pipeline LiDAR pour détection archéologique - Version Python pur -Visualisations générées avec numpy/rasterio/matplotlib +Pipeline LiDAR pour détection archéologique +Visualisations générées avec numpy/cupy + rasterio/matplotlib +Support GPU via CuPy si disponible, fallback numpy sinon """ import os @@ -33,6 +34,55 @@ try: except ImportError: HAS_WARP = False +# GPU acceleration via CuPy +try: + import cupy as cp + import cupyx.scipy.ndimage as cp_ndimage + HAS_GPU = True + _xp = cp # Default array module (GPU) + _gpu_info = cp.cuda.runtime.getDeviceProperties(0) + _gpu_name = _gpu_info['name'].decode() if isinstance(_gpu_info['name'], bytes) else str(_gpu_info['name']) + print(f"✓ GPU détectée: {_gpu_name}") + print(f" Mémoire GPU: {_gpu_info['totalGlobalMem'] // (1024**3)} Go") +except (ImportError, Exception): + HAS_GPU = False + _xp = np # Fallback CPU + + +def to_gpu(arr): + """Send array to GPU if available.""" + if HAS_GPU: + return cp.asarray(arr.astype(np.float64)) + return arr.astype(np.float64) + + +def to_cpu(arr): + """Bring array back to CPU (numpy).""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp.asnumpy(arr) + return arr + + +def xp_gaussian_filter(arr, sigma): + """Gaussian filter using GPU if available.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp_ndimage.gaussian_filter(arr, sigma) + return ndimage.gaussian_filter(arr, sigma) + + +def xp_uniform_filter(arr, size): + """Uniform filter using GPU if available.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp_ndimage.uniform_filter(arr, size) + return ndimage.uniform_filter(arr, size) + + +def xp_minimum_filter(arr, footprint=None, size=None): + """Minimum filter using GPU if available.""" + if HAS_GPU and isinstance(arr, cp.ndarray): + return cp_ndimage.minimum_filter(arr, footprint=footprint, size=size) + return ndimage.minimum_filter(arr, footprint=footprint, size=size) + rcParams['figure.dpi'] = 150 rcParams['savefig.dpi'] = 300 rcParams['font.size'] = 10 @@ -586,37 +636,36 @@ class LidarArchaeoPipeline: def generate_lrm(self, dem_file, basename): """Local Relief Model - deviation from local mean""" - print(f" → Local Relief Model...") + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → Local Relief Model{gpu_tag}...") output = self.vis_dir / f"{basename}_lrm.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - # Calculate local mean with gaussian filter - from scipy.ndimage import gaussian_filter - local_mean = gaussian_filter(dem, sigma=15/self.resolution) - - # Deviation from mean + dem = to_gpu(dem_np) + local_mean = xp_gaussian_filter(dem, sigma=15/self.resolution) lrm = dem - local_mean + lrm_np = to_cpu(lrm).astype(np.float32) # Save with rasterio.open( output, 'w', driver='GTiff', - height=lrm.shape[0], - width=lrm.shape[1], + height=lrm_np.shape[0], + width=lrm_np.shape[1], count=1, dtype='float32', crs=crs, transform=transform, compress='lzw' ) as dst: - dst.write(lrm.astype('float32'), 1) + dst.write(lrm_np, 1) return output except Exception as e: @@ -627,85 +676,82 @@ class LidarArchaeoPipeline: """Sky-View Factor - ray-tracing on 16 azimuths. 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.""" - print(f" → Sky-View Factor (ray-tracing)...") + Valleys/crevices have low SVF (obstructed sky), ridges/peaks have high SVF. + Uses GPU (CuPy) if available for acceleration.""" + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → Sky-View Factor (ray-tracing){gpu_tag}...") output = self.vis_dir / f"{basename}_svf.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - rows, cols = dem.shape - res = self.resolution # meters per pixel + rows, cols = dem_np.shape + res = self.resolution - # 16 azimuth directions (0°, 22.5°, 45°, ... 337.5°) + # Move to GPU if available + dem = to_gpu(dem_np) + + # 16 azimuth directions n_dirs = 16 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx = np.cos(angles) dy = np.sin(angles) - # Maximum search distance (in pixels) max_dist = int(50 / res) # 50m search radius - # Pad DEM with NaN to handle boundaries - padded = np.pad(dem.astype(np.float64), max_dist, mode='constant', - constant_values=np.nan) + # Pad DEM with NaN for boundary handling + xp = cp if HAS_GPU else np + padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - svf = np.zeros_like(dem, dtype=np.float64) + svf = xp.zeros_like(dem) for d_idx in range(n_dirs): - # Direction vector ddx, ddy = dx[d_idx], dy[d_idx] - - # For each step along the ray, compute the horizon angle - horizon = np.zeros_like(dem, dtype=np.float64) + horizon = xp.zeros_like(dem) for step in range(1, max_dist + 1): - # Offset in pixels (rounded to nearest integer) px = int(round(ddx * step)) py = int(round(ddy * step)) - - # Distance in meters dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2) if dist_m < res * 0.5: continue - # Elevation difference relative to center pixel - # Source is at (max_dist + row, max_dist + col) in padded array - # Target is at (max_dist + row + py, max_dist + col + px) elev_diff = padded[max_dist + py:max_dist + py + rows, max_dist + px:max_dist + px + cols] - dem - # Horizon angle = arctan(elev_diff / dist) - angle = np.arctan2(elev_diff, dist_m) + angle = xp.arctan2(elev_diff, dist_m) + horizon = xp.where(xp.isnan(angle), horizon, + xp.maximum(horizon, xp.nan_to_num(angle, nan=0))) - # Keep maximum horizon angle - horizon = np.where(np.isnan(angle), horizon, - np.maximum(horizon, np.nan_to_num(angle, nan=0))) + svf += xp.cos(xp.pi / 2 - horizon) ** 2 - # SVF contribution: cos²(zenith) where zenith = pi/2 - horizon - # Higher horizon = less sky visible - svf += np.cos(np.pi / 2 - horizon) ** 2 + svf /= n_dirs - svf /= n_dirs # Average over all directions + # Bring back to CPU for saving + svf_np = to_cpu(svf).astype(np.float32) - # Save with rasterio.open( output, 'w', driver='GTiff', - height=svf.shape[0], - width=svf.shape[1], + height=svf_np.shape[0], + width=svf_np.shape[1], count=1, dtype='float32', crs=crs, transform=transform, compress='lzw' ) as dst: - dst.write(svf.astype('float32'), 1) + dst.write(svf_np, 1) + + return output + except Exception as e: + print(f" ✗ Erreur SVF: {e}") + return None return output except Exception as e: @@ -717,46 +763,44 @@ class LidarArchaeoPipeline: 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.""" + Result is averaged across all 8 directions. + Uses GPU (CuPy) if available for acceleration.""" name = "positive_openness" if positive else "negative_openness" - print(f" → {name.replace('_', ' ').title()} (ray-tracing)...") + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → {name.replace('_', ' ').title()} (ray-tracing){gpu_tag}...") output = self.vis_dir / f"{basename}_{name}.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - rows, cols = dem.shape + rows, cols = dem_np.shape res = self.resolution - # 8 cardinal directions + dem = to_gpu(dem_np) + xp = cp if HAS_GPU else np + n_dirs = 8 angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False) dx = np.cos(angles) dy = np.sin(angles) - max_dist = int(50 / res) # 50m search radius + max_dist = int(50 / res) - # Pad DEM with NaN to handle boundaries - padded = np.pad(dem.astype(np.float64), max_dist, mode='constant', - constant_values=np.nan) + padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan) - openness_sum = np.zeros_like(dem, dtype=np.float64) + openness_sum = xp.zeros_like(dem) for d_idx in range(n_dirs): ddx, ddy = dx[d_idx], dy[d_idx] - - # For positive openness: find max upward angle (zenith) - # For negative openness: find max downward angle (nadir) - max_angle = np.zeros_like(dem, dtype=np.float64) + max_angle = 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 @@ -765,23 +809,17 @@ class LidarArchaeoPipeline: max_dist + px:max_dist + px + cols] - dem if positive: - # Positive openness: angle from vertical to terrain above - # Only consider points higher than center (elev_diff > 0) - angle = np.arctan2(np.maximum(elev_diff, 0), dist_m) + angle = xp.arctan2(xp.maximum(elev_diff, 0), dist_m) else: - # Negative openness: angle from vertical to terrain below - # Only consider points lower than center (elev_diff < 0) - angle = np.arctan2(np.maximum(-elev_diff, 0), dist_m) + angle = xp.arctan2(xp.maximum(-elev_diff, 0), dist_m) - max_angle = np.where(np.isnan(angle), max_angle, - np.maximum(max_angle, np.nan_to_num(angle, nan=0))) + max_angle = xp.where(xp.isnan(angle), max_angle, + xp.maximum(max_angle, xp.nan_to_num(angle, nan=0))) openness_sum += max_angle # Average across all directions (convert to degrees) - openness_result = np.degrees(openness_sum / n_dirs) - - # Save + openness_result = to_cpu(xp.degrees(openness_sum / n_dirs)).astype(np.float32) with rasterio.open( output, 'w', @@ -804,18 +842,21 @@ class LidarArchaeoPipeline: def generate_mslrm(self, dem_file, basename): """Multi-Scale Relief Model (MSRM) - LRM computed at multiple scales (sigma=5,10,25,50,100) and combined. Reveals features at all scales: - small (walls, ditches), medium (enclosures, tumulus), large (landscapes).""" - print(f" → Multi-Scale Relief Model (MSRM)...") + small (walls, ditches), medium (enclosures, tumulus), large (landscapes). + Uses GPU if available for acceleration.""" + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → Multi-Scale Relief Model (MSRM){gpu_tag}...") output = self.vis_dir / f"{basename}_mslrm.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - from scipy.ndimage import gaussian_filter + dem = to_gpu(dem_np) + xp = cp if HAS_GPU else np # Compute LRM at multiple scales sigmas = [5, 10, 25, 50, 100] @@ -823,28 +864,28 @@ class LidarArchaeoPipeline: for sigma in sigmas: sigma_px = sigma / self.resolution - local_mean = gaussian_filter(dem, sigma=sigma_px) + local_mean = xp_gaussian_filter(dem, sigma=sigma_px) lrm = dem - local_mean - # Normalize each scale - lrm_norm = lrm / max(np.nanstd(lrm), 0.01) + lrm_norm = lrm / max(xp.nanstd(lrm), 0.01) lrm_stack.append(lrm_norm) # Combine: RMS of normalized LRM at all scales - mslrm = np.sqrt(np.mean(np.array(lrm_stack) ** 2, axis=0)) + mslrm = xp.sqrt(xp.mean(xp.array(lrm_stack) ** 2, axis=0)) + mslrm_np = to_cpu(mslrm).astype(np.float32) with rasterio.open( output, 'w', driver='GTiff', - height=mslrm.shape[0], - width=mslrm.shape[1], + height=mslrm_np.shape[0], + width=mslrm_np.shape[1], count=1, dtype='float32', crs=crs, transform=transform, compress='lzw' ) as dst: - dst.write(mslrm.astype('float32'), 1) + dst.write(mslrm_np, 1) return output except Exception as e: @@ -855,50 +896,52 @@ class LidarArchaeoPipeline: """Multi-Scale Topographic Position Index. TPI = elevation - mean(neighborhood). Computed at fine (5m) and broad (100m) scales to identify - ridges, valleys, and intermediate landforms.""" - print(f" → TPI multi-echelle...") + ridges, valleys, and intermediate landforms. + Uses GPU if available for acceleration.""" + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → TPI multi-echelle{gpu_tag}...") output = self.vis_dir / f"{basename}_tpi.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - from scipy.ndimage import uniform_filter + dem = to_gpu(dem_np) - # Fine-scale TPI (5m radius) fine_size = int(5 / self.resolution) if fine_size % 2 == 0: fine_size += 1 - tpi_fine = dem - uniform_filter(dem, size=fine_size) + tpi_fine = dem - xp_uniform_filter(dem, size=fine_size) - # Broad-scale TPI (100m radius) broad_size = int(100 / self.resolution) if broad_size % 2 == 0: broad_size += 1 - tpi_broad = dem - uniform_filter(dem, size=broad_size) + tpi_broad = dem - xp_uniform_filter(dem, size=broad_size) # Combine: fine-scale weighted more for archaeological features # Normalize each scale then weight: 0.6 fine + 0.4 broad - fine_std = max(np.nanstd(tpi_fine), 0.01) - broad_std = max(np.nanstd(tpi_broad), 0.01) + xp = cp if HAS_GPU else np + fine_std = max(float(xp.nanstd(tpi_fine)), 0.01) + broad_std = max(float(xp.nanstd(tpi_broad)), 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) with rasterio.open( output, 'w', driver='GTiff', - height=tpi_combined.shape[0], - width=tpi_combined.shape[1], + height=tpi_np.shape[0], + width=tpi_np.shape[1], count=1, dtype='float32', crs=crs, transform=transform, compress='lzw' ) as dst: - dst.write(tpi_combined.astype('float32'), 1) + dst.write(tpi_np, 1) return output except Exception as e: @@ -1080,62 +1123,57 @@ class LidarArchaeoPipeline: def generate_sailore(self, dem_file, basename): """SAILORE - Self-Adaptive Improved Local Relief Model. Kernel size adapts to local slope: flat areas get larger kernels, - steep areas get smaller kernels. Better than fixed LRM for heterogeneous terrain.""" - print(f" → SAILORE (LRM adaptatif)...") + steep areas get smaller kernels. Better than fixed LRM for heterogeneous terrain. + Uses GPU if available for acceleration.""" + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → SAILORE (LRM adaptatif){gpu_tag}...") output = self.vis_dir / f"{basename}_sailore.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - from scipy.ndimage import gaussian_filter, uniform_filter + dem = to_gpu(dem_np) + xp = cp if HAS_GPU else np - # Compute slope for adaptive kernel sizing - gy, gx = np.gradient(dem, self.resolution) - slope = np.arctan(np.sqrt(gx**2 + gy**2)) - slope_deg = np.degrees(slope) + gy, gx = xp.gradient(dem, self.resolution) + slope = xp.arctan(xp.sqrt(gx**2 + gy**2)) + slope_deg = xp.degrees(slope) - # Adaptive sigma: flat terrain (low slope) = large kernel, steep = small - # slope 0° -> sigma_max (25m), slope 30° -> sigma_min (2m) sigma_min = 2.0 / self.resolution sigma_max = 25.0 / self.resolution - # Normalize slope to 0-1 range - slope_norm = np.clip(slope_deg / 30.0, 0, 1) - # Invert: flat areas get high sigma, steep areas get low sigma + slope_norm = xp.clip(slope_deg / 30.0, 0, 1) adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min) - # Compute local mean with adaptive Gaussian (approximated by blending) - # For efficiency, compute at 3 fixed scales and blend based on slope - lrm_fine = dem - gaussian_filter(dem, sigma=sigma_min) - lrm_medium = dem - gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2) - lrm_coarse = dem - gaussian_filter(dem, sigma=sigma_max) + 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) - # Blend based on slope: steep -> fine, flat -> coarse w_fine = slope_norm - w_medium = 1 - 2 * np.abs(slope_norm - 0.5) + w_medium = 1 - 2 * xp.abs(slope_norm - 0.5) w_coarse = 1 - slope_norm - # Normalize weights 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) with rasterio.open( output, 'w', driver='GTiff', - height=sailore.shape[0], - width=sailore.shape[1], + height=sailore_np.shape[0], + width=sailore_np.shape[1], count=1, dtype='float32', crs=crs, transform=transform, compress='lzw' ) as dst: - dst.write(sailore.astype('float32'), 1) + dst.write(sailore_np, 1) return output except Exception as e: @@ -1361,34 +1399,40 @@ class LidarArchaeoPipeline: def generate_wavelet(self, dem_file, basename): """Mexican Hat wavelet (Ricker wavelet) multi-scale analysis. Continuous Wavelet Transform at multiple scales to detect - circular/circular features like tumulus, ring ditches, enclosures.""" - print(f" → Ondelette Mexican Hat multi-echelle...") + circular/circular features like tumulus, ring ditches, enclosures. + Uses GPU if available for acceleration.""" + gpu_tag = " [GPU]" if HAS_GPU else "" + print(f" → Ondelette Mexican Hat multi-echelle{gpu_tag}...") output = self.vis_dir / f"{basename}_wavelet.tif" try: with rasterio.open(dem_file) as src: - dem = src.read(1) + dem_np = src.read(1) transform = src.transform crs = src.crs - # Mexican Hat (Ricker) wavelet at multiple scales - # Uses scipy.ndimage.gaussian_laplace as 2D Mexican Hat approximation + dem = to_gpu(dem_np) + xp = cp if HAS_GPU else np + scales = [2, 5, 10, 20, 50] # meters wavelet_stack = [] for scale_m in scales: - # Create 1D Ricker wavelet and apply as 2D separable filter sigma_px = scale_m / self.resolution # 2D Mexican Hat = Laplacian of Gaussian - from scipy.ndimage import gaussian_laplace - response = -gaussian_laplace(dem.astype(np.float64), sigma=sigma_px) - # Normalize by scale - response /= max(np.nanstd(response), 0.01) + if HAS_GPU: + from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace + response = -gpu_gaussian_laplace(dem, sigma=sigma_px) + else: + from scipy.ndimage import gaussian_laplace + response = -gaussian_laplace(to_cpu(dem).astype(np.float64), sigma=sigma_px) + response = to_gpu(response) + response /= max(float(xp.nanstd(response)), 0.01) wavelet_stack.append(response) - # Combine: RMS of all scales - combined = np.sqrt(np.mean(np.array(wavelet_stack) ** 2, axis=0)) + combined = xp.sqrt(xp.mean(xp.array(wavelet_stack) ** 2, axis=0)) + combined_np = to_cpu(combined).astype(np.float32) with rasterio.open( output, @@ -1402,7 +1446,7 @@ class LidarArchaeoPipeline: transform=transform, compress='lzw' ) as dst: - dst.write(combined.astype('float32'), 1) + dst.write(combined_np, 1) return output except Exception as e: @@ -1935,7 +1979,7 @@ class LidarArchaeoPipeline: if not tif_file or not tif_file.exists(): return None - jpg_file = self.vis_dir / f"{tif_file.stem}.png" + jpg_file = self.vis_dir / f"{tif_file.stem}.webp" try: with rasterio.open(tif_file) as src: @@ -2331,10 +2375,18 @@ class LidarArchaeoPipeline: fig.patch.set_facecolor('white') - plt.savefig(jpg_file, dpi=150, bbox_inches='tight', pad_inches=0.15, + # Save as PNG first (matplotlib doesn't support WebP well) + png_temp = self.vis_dir / f"{tif_file.stem}_temp.png" + plt.savefig(png_temp, dpi=150, bbox_inches='tight', pad_inches=0.15, facecolor='white', format='png') plt.close() + # Convert PNG to lossless WebP using PIL + from PIL import Image as PILImage + img = PILImage.open(str(png_temp)) + img.save(str(jpg_file), format='WEBP', lossless=True) + png_temp.unlink() # Delete temp PNG + # Delete the source TIFF file to save space tif_file.unlink() @@ -2425,9 +2477,9 @@ class LidarArchaeoPipeline: # Look for PNGs in per-file subdirectory first, then fallback to main dir file_vis_dir = self.vis_dir / basename if file_vis_dir.exists(): - png_files = sorted(file_vis_dir.glob("*.png")) + png_files = sorted(file_vis_dir.glob("*.webp")) else: - png_files = sorted(self.vis_dir.glob(f"{basename}_*.png")) + png_files = sorted(self.vis_dir.glob(f"{basename}_*.webp")) if not png_files: print(f" ✗ Aucune image PNG trouvée pour {basename}") return None diff --git a/run.sh b/run.sh new file mode 100755 index 0000000..aceec6e --- /dev/null +++ b/run.sh @@ -0,0 +1,72 @@ +#!/bin/bash +# Pipeline LiDAR Archéologique +# Utilisation: ./run.sh [options] +# Options: +# -r RESOLUTION Résolution en m/px (défaut: 0.5) +# -w WORKERS Nombre de workers parallèles (défaut: 1) +# -g Activer l'accélération GPU +# -h Afficher l'aide + +set -e + +SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)" +INPUT_DIR="${SCRIPT_DIR}/input" +OUTPUT_DIR="${SCRIPT_DIR}/output" +IMAGE_NAME="lidar-lidar" +RESOLUTION=0.5 +WORKERS=1 +GPU_FLAG="" + +while getopts "r:w:gh" opt; do + case $opt in + r) RESOLUTION="$OPTARG" ;; + w) WORKERS="$OPTARG" ;; + g) GPU_FLAG="--gpus all" ;; + h) + echo "Pipeline LiDAR Archéologique" + echo "" + echo "Usage: $0 [-r RESOLUTION] [-w WORKERS] [-g]" + echo "" + echo " -r RESOLUTION Résolution en m/px (défaut: 0.5)" + echo " -w WORKERS Nombre de workers CPU parallèles (défaut: 1)" + echo " -g Activer l'accélération GPU NVIDIA" + echo " -h Afficher cette aide" + echo "" + echo "Exemples:" + echo " $0 # Traitement standard" + echo " $0 -g # Avec accélération GPU" + echo " $0 -g -w 4 # GPU + 4 workers parallèles" + echo " $0 -r 0.2 -g # Haute résolution + GPU" + exit 0 + ;; + *) echo "Option invalide. Utilisez -h pour l'aide." >&2; exit 1 ;; + esac +done + +# Build l'image si elle n'existe pas +if ! docker image inspect "$IMAGE_NAME" >/dev/null 2>&1; then + echo "Build de l'image Docker..." + docker build -t "$IMAGE_NAME" "$SCRIPT_DIR" +fi + +# Créer les répertoires s'ils n'existent pas +mkdir -p "$INPUT_DIR" "$OUTPUT_DIR" + +# Lancer le pipeline +echo "============================================" +echo " Pipeline LiDAR Archéologique" +echo "============================================" +echo " Résolution : ${RESOLUTION}m/px" +echo " Workers : ${WORKERS}" +echo " GPU : $([ -n "$GPU_FLAG" ] && echo 'OUI' || echo 'non')" +echo "============================================" + +docker run --rm $GPU_FLAG \ + --user 1000:1000 \ + -v "${INPUT_DIR}:/data/input:ro" \ + -v "${OUTPUT_DIR}:/data/output" \ + "$IMAGE_NAME" \ + python3 /usr/local/bin/process_lidar.py /data/input \ + -o /data/output \ + -r "$RESOLUTION" \ + -w "$WORKERS" \ No newline at end of file