diff --git a/CLAUDE.md b/CLAUDE.md
index 279791e..c50821a 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -4,7 +4,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 
 ## Project Overview
 
-LiDAR archaeological processing pipeline that generates 20 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
+LiDAR archaeological processing pipeline that generates 16 terrain visualizations from LAZ/LAS point clouds. Runs in Docker with optional NVIDIA GPU acceleration (CuPy). Designed for French LiDAR HD data in Lambert 93 (EPSG:2154).
 
 ## Commands
 
@@ -14,10 +14,15 @@ All commands run inside Docker. Use `./run.sh` as the primary interface.
 ./run.sh -g                                    # Standard run with GPU
 ./run.sh -g -w 4                               # GPU + 4 parallel workers
 ./run.sh -g -r 0.2                             # High resolution (0.2m/px)
+./run.sh -g -r 0.5,0.2                         # Multi-resolution (0.5m + 0.2m)
 ./run.sh --test                                 # Run unit tests
 ./run.sh -g --file LHD_FXX_1000_6882_PTS_LAMB93_IGN69.copc  # Single file
 ./run.sh --ground-classification csf            # Force CSF ground classification (complex terrain)
 ./run.sh -g --keep-tif                          # Keep TIFF files (allows WebP regeneration without recalculating DTM)
+./run.sh -g --only hillshade svf lrm            # Only generate specific visualizations
+./run.sh -g --skip ortho topo                   # Exclude specific visualizations
+./run.sh -g --quality 90                       # WebP quality 90 (default: 85)
+./run.sh -g --lossless                          # Lossless WebP compression
 ./run.sh                                        # Print help (no args)
 ```
 
@@ -32,12 +37,12 @@ docker run --rm --gpus all -v $(pwd)/input:/data/input:ro -v $(pwd)/output:/data
 ### Module responsibilities
 
 - **`cli.py`** — argparse + logging setup. Entry point via `python -m lidar_pipeline`.
-- **`pipeline.py`** — `LidarArchaeoPipeline` orchestrator. `VIZ_STEPS` registry maps names to generate functions. `FilePrefixFilter` for parallel logging. Creates `SharedDEM` once per file and passes it to all visualizations.
-- **`dtm.py`** — PDAL ground classification (SMRF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`.
-- **`visualizations.py`** — 15 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution, shared=None)` and return a TIF path or None. `SharedDEM` class pre-computes gradient, NaN mask, LRM to avoid redundant I/O and computation.
+- **`pipeline.py`** — `LidarArchaeoPipeline` orchestrator. `VIZ_STEPS` registry maps names to generate functions. `FilePrefixFilter` for parallel logging. Creates `SharedDEM` once per file and passes it to all visualizations. Multi-resolution support: `self.resolutions` list, `_res_suffix()` for naming, `generate_all_visualizations()` accepts `vis_dir` override.
+- **`dtm.py`** — PDAL ground classification (SMRF/CSF + auto-detection) and DTM generation via scipy `binned_statistic_2d`. `create_dtm_fast()` accepts `output_suffix` for multi-resolution DTM naming.
+- **`visualizations.py`** — 13 `generate_*` functions + 2 IGN overlay lambdas. All take `(dem_file, basename, vis_dir, resolution, shared=None)` and return a TIF path or None. `SharedDEM` class pre-computes gradient, NaN mask, LRM to avoid redundant I/O and computation. Lazy evaluation: properties computed on first access.
 - **`gpu.py`** — CuPy/numpy abstraction: `HAS_GPU`, `to_gpu()`, `to_cpu()`, `xp_gaussian_filter()`, `xp_uniform_filter()`, `xp_minimum_filter()`, `gpu_cleanup()`. Falls back to CPU gracefully.
 - **`ign.py`** — IGN WMTS tile download + overlay generation for orthophoto and topographic maps.
-- **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. `generate_pdf_report()` creates A3 PDF.
+- **`rendering.py`** — `COLORMAPS` dict maps filename keywords to (cmap, title, legend, description). `tif_to_png()` converts TIF→WebP with legend/scale/north arrow. Quality parameter controls WebP compression (default 85).
 
 ### SharedDEM optimization
 
@@ -73,6 +78,10 @@ DTM small gaps (< 1m from existing data) are filled using `rasterio.fill.fillnod
 
 Uses priority-flood algorithm (Wang & Liu 2006) for sink filling, which is O(n log n) instead of iterative minimum_filter. D8 accumulation uses numba JIT; falls back to pure Python if numba unavailable.
 
+### Multi-resolution
+
+`-r 0.5,0.2` processes each tile at both 0.5m and 0.2m. Ground classification is shared (done once per tile). Each resolution gets its own DTM (`_dtm.tif` / `_dtm_r0p2.tif`) and visualization subdirectory (`basename/` / `basename_r0p2/`).
+
 ### Parallel processing
 
 Uses `ProcessPoolExecutor` with `'spawn'` start method (required for CUDA). Each worker gets its own temp directory (`temp_{basename}`). `_process_file_standalone()` configures its own logger with `_file_filter` for per-file log prefixes.
@@ -82,6 +91,6 @@ Uses `ProcessPoolExecutor` with `'spawn'` start method (required for CUDA). Each
 - **Language**: UI messages and comments in French. Code identifiers in English.
 - **Logging**: Use `logger = logging.getLogger("lidar")`. Prefix per-file logs via `_file_filter.basename`.
 - **GPU pattern**: `arr_gpu = to_gpu(arr)` → compute → `result = to_cpu(arr_gpu)` → `gpu_cleanup()` between visualizations.
-- **Output format**: Visualizations saved as WebP. TIFF intermediates deleted by default. Use `--keep-tif` to keep DTM+TIF for WebP regeneration with `--force`. No COGs or viewer.
+- **Output format**: Visualizations saved as AVIF (quality 98 by default, best quality/size ratio). Use `--format webp` for WebP output. TIFF intermediates deleted by default. Use `--keep-tif` to keep DTM+TIF for regeneration with `--force`. No PDF reports, no COGs or viewer.
 - **Compression**: TIF intermediates use `deflate` compression (faster than LZW for float32 data).
 - **Tests**: Run only inside Docker via `./run.sh --test`. Synthetic DEM fixture in `tests/conftest.py`.
\ No newline at end of file
diff --git a/Dockerfile b/Dockerfile
index 5094986..f46d1ba 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -41,7 +41,8 @@ RUN pip3 install --no-cache-dir \
     titiler.core \
     fastapi \
     uvicorn \
-    piexif
+    piexif \
+    pillow-avif-plugin
 
 # 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"
diff --git a/lidar_pipeline/cli.py b/lidar_pipeline/cli.py
index 98e3766..8e0855e 100644
--- a/lidar_pipeline/cli.py
+++ b/lidar_pipeline/cli.py
@@ -131,13 +131,19 @@ Exemples:
     parser.add_argument(
         "--quality",
         type=int,
-        default=85,
-        help="Qualité WebP (1-100, défaut: 85). Utilisez 100 pour lossless."
+        default=98,
+        help="Qualité image (1-100, défaut: 98). Utilisez 100 pour lossless."
     )
     parser.add_argument(
         "--lossless",
         action="store_true",
-        help="Forcer la compression WebP lossless (équivalent à --quality 100)"
+        help="Forcer la compression lossless (équivalent à --quality 100)"
+    )
+    parser.add_argument(
+        "--format",
+        choices=["webp", "avif"],
+        default="avif",
+        help="Format de sortie : avif (défaut, meilleure qualité) ou webp"
     )
     parser.add_argument(
         "--only",
@@ -194,6 +200,13 @@ Exemples:
 
     try:
         quality = 100 if args.lossless else args.quality
+        # Parse --only and --skip: accept comma-separated values
+        only_viz = None
+        if args.only:
+            only_viz = [v.strip() for item in args.only for v in item.split(',')]
+        skip_viz = None
+        if args.skip:
+            skip_viz = [v.strip() for item in args.skip for v in item.split(',')]
         pipeline = LidarArchaeoPipeline(
             input_dir=args.input,
             output_dir=args.output,
@@ -204,8 +217,9 @@ Exemples:
             force_classify=args.force_classification,
             keep_tif=args.keep_tif,
             quality=quality,
-            only_viz=args.only,
-            skip_viz=args.skip,
+            only_viz=only_viz,
+            skip_viz=skip_viz,
+            output_format=args.format,
         )
 
         # If --file is specified, process only matching files
@@ -270,9 +284,15 @@ def _kill_orphan_pdal(signum=None, frame=None):
     """Kill orphan PDAL processes on interrupt or exit."""
     import subprocess
     try:
-        subprocess.run(["pkill", "-f", "pdal"], capture_output=True, timeout=5)
+        subprocess.run(["pkill", "-9", "-f", "pdal"], capture_output=True, timeout=3)
     except Exception:
         pass
     if signum is not None:
-        logger.info("Interruption — nettoyage des processus PDAL")
+        logger.info("Interruption — nettoyage des processus")
+        # Force-kill all child processes immediately
+        try:
+            import os
+            os.killpg(os.getpgrp(), signal.SIGKILL)
+        except Exception:
+            pass
         sys.exit(130)
\ No newline at end of file
diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py
index ae64972..9e7f15c 100644
--- a/lidar_pipeline/dtm.py
+++ b/lidar_pipeline/dtm.py
@@ -483,7 +483,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution, force=False, output
             logger.info(f"    {remaining:,} pixels restent sans données (grands écarts)")
 
         # Save as GeoTIFF
-        output_tif = dtm_dir / f"{basename}_dtm.tif"
+        output_tif = dtm_dir / f"{basename}_dtm{output_suffix}.tif"
         transform = from_bounds(min_x, min_y, max_x, max_y, width, height)
 
         with rasterio.open(
diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py
index 85cdc49..0f4f23e 100644
--- a/lidar_pipeline/pipeline.py
+++ b/lidar_pipeline/pipeline.py
@@ -58,10 +58,10 @@ from .dtm import classify_ground, create_dtm_fast
 from .visualizations import (
     SharedDEM,
     generate_hillshade, generate_slope, generate_aspect, generate_curvature,
-    generate_lrm, generate_svf, generate_openness,
+    generate_lrm, generate_openness,
     generate_mslrm, generate_tpi, generate_sailore,
     generate_roughness, generate_anomalies, generate_wavelet,
-    generate_flow, generate_local_dominance,
+    generate_flow,
 )
 from .gpu import gpu_cleanup
 from .ign import generate_ign_overlay
@@ -76,7 +76,6 @@ VIZ_STEPS = [
     ('slope',      generate_slope),
     ('aspect',     generate_aspect),
     ('curvature',  generate_curvature),
-    ('svf',        generate_svf),
     ('lrm',        generate_lrm),
     ('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)),
@@ -87,7 +86,6 @@ VIZ_STEPS = [
     ('anomalies',  generate_anomalies),
     ('wavelet',    generate_wavelet),
     ('flow',       generate_flow),
-    ('local_dominance', generate_local_dominance),
     ('ortho',      lambda d, b, v, r: generate_ign_overlay(
                     d, b, v, r,
                     layer='ORTHOIMAGERY.ORTHOPHOTOS',
@@ -108,7 +106,7 @@ VIZ_STEPS = [
 class LidarArchaeoPipeline:
     """Orchestrates the LiDAR archaeological analysis pipeline."""
 
-    def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False, quality=85, only_viz=None, skip_viz=None):
+    def __init__(self, input_dir, output_dir, resolution=0.5, workers=1, force=False, ground_method='auto', force_classify=False, keep_tif=False, quality=98, only_viz=None, skip_viz=None, output_format='avif'):
         self.input_dir = Path(input_dir)
         self.output_dir = Path(output_dir)
         # Accept single float or comma-separated string for multi-resolution
@@ -127,6 +125,7 @@ class LidarArchaeoPipeline:
         self.quality = quality
         self.only_viz = only_viz
         self.skip_viz = skip_viz
+        self.output_format = output_format
         self.temp_dir = self.output_dir / "temp"
 
         if not self.input_dir.exists():
@@ -137,9 +136,8 @@ class LidarArchaeoPipeline:
 
         self.dtm_dir = self.output_dir / "DTM"
         self.vis_dir = self.output_dir / "visualisations"
-        self.pdf_dir = self.output_dir / "rapports"
 
-        for d in [self.dtm_dir, self.vis_dir, self.pdf_dir]:
+        for d in [self.dtm_dir, self.vis_dir]:
             d.mkdir(exist_ok=True)
 
         # Filter visualizations based on --only / --skip
@@ -169,7 +167,7 @@ class LidarArchaeoPipeline:
         logger.info(f"  Classification sol : {self.ground_method}")
         logger.info(f"  Force classif.: {'OUI' if self.force_classify else 'non'}")
         logger.info(f"  Keep TIFF    : {'OUI' if self.keep_tif else 'non'}")
-        logger.info(f"  Qualité WebP: {self.quality if self.quality < 100 else 'lossless'}")
+        logger.info(f"  Qualité {self.output_format.upper()}: {self.quality if self.quality < 100 else 'lossless'}")
         if only_viz:
             logger.info(f"  Visualisations: uniquement {', '.join(only_viz)}")
         elif skip_viz:
@@ -197,18 +195,19 @@ class LidarArchaeoPipeline:
         return True
 
     @staticmethod
-    def _expected_webp_path(name, basename, file_vis_dir):
-        """Return the expected WebP filename for a visualization step."""
+    def _expected_output_path(name, basename, file_vis_dir, output_format='avif'):
+        """Return the expected output filename for a visualization step."""
+        ext = 'avif' if output_format == 'avif' else 'webp'
         if name == 'pos_open':
-            return file_vis_dir / f"{basename}_positive_openness.webp"
+            return file_vis_dir / f"{basename}_positive_openness.{ext}"
         elif name == 'neg_open':
-            return file_vis_dir / f"{basename}_negative_openness.webp"
+            return file_vis_dir / f"{basename}_negative_openness.{ext}"
         elif name == 'hillshade':
-            return file_vis_dir / f"{basename}_hillshade_multi.webp"
+            return file_vis_dir / f"{basename}_hillshade_multi.{ext}"
         else:
-            return file_vis_dir / f"{basename}_{name}.webp"
+            return file_vis_dir / f"{basename}_{name}.{ext}"
 
-    def generate_all_visualizations(self, dtm_file, basename, resolution=None):
+    def generate_all_visualizations(self, dtm_file, basename, resolution=None, vis_dir=None):
         """Generate all archaeological visualizations for one DTM file.
 
         Optimisation: SharedDEM is only computed if at least one visualization
@@ -219,8 +218,8 @@ class LidarArchaeoPipeline:
             resolution = self.resolution
         logger.info("  Génération visualisations:")
 
-        # Create per-file subdirectory
-        file_vis_dir = self.vis_dir / basename
+        # Use provided vis_dir (for multi-resolution subdirectories) or default
+        file_vis_dir = vis_dir if vis_dir else (self.vis_dir / basename)
         file_vis_dir.mkdir(exist_ok=True)
         total = len(self.viz_steps)
 
@@ -230,7 +229,7 @@ class LidarArchaeoPipeline:
             if self.force:
                 needs_generation[name] = True
             else:
-                expected_webp = self._expected_webp_path(name, basename, file_vis_dir)
+                expected_webp = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
                 needs_generation[name] = not expected_webp.exists()
 
         to_generate = [n for n, needed in needs_generation.items() if needed]
@@ -242,7 +241,7 @@ class LidarArchaeoPipeline:
             # Still need to return results dict for PDF check
             vis_results = {}
             for name, func in self.viz_steps:
-                vis_results[name] = self._expected_webp_path(name, basename, file_vis_dir)
+                vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
             return vis_results
 
         # Phase 2: compute SharedDEM only if needed
@@ -258,7 +257,7 @@ class LidarArchaeoPipeline:
         for idx, (name, func) in enumerate(self.viz_steps, 1):
             if not needs_generation[name]:
                 logger.info(f"    [{idx}/{total}] {name}: déjà existant, ignoré")
-                vis_results[name] = self._expected_webp_path(name, basename, file_vis_dir)
+                vis_results[name] = self._expected_output_path(name, basename, file_vis_dir, self.output_format)
                 continue
 
             # When --force, delete existing TIF to ensure clean regeneration
@@ -296,8 +295,9 @@ class LidarArchaeoPipeline:
             # Free GPU memory between visualizations to prevent OOM
             gpu_cleanup()
 
-        # Convert to WebP (only newly generated TIFs, not skipped ones)
-        logger.info("  Conversion images WebP:")
+        # Convert to output format (only newly generated TIFs, not skipped ones)
+        fmt_label = self.output_format.upper()
+        logger.info(f"  Conversion images {fmt_label}:")
         source_info = {
             'method': self.ground_method,
             'date': datetime.now().strftime('%Y-%m-%d'),
@@ -305,9 +305,9 @@ class LidarArchaeoPipeline:
         }
         for name, tif_file in vis_results.items():
             if tif_file and isinstance(tif_file, Path) and tif_file.suffix == '.tif' and tif_file.exists():
-                webp_file = tif_to_png(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, source_info=source_info, quality=self.quality)
-                if webp_file:
-                    logger.info(f"    ✓ {webp_file.name}")
+                img_file = tif_to_png(tif_file, file_vis_dir, resolution, keep_tif=self.keep_tif, source_info=source_info, quality=self.quality, output_format=self.output_format)
+                if img_file:
+                    logger.info(f"    ✓ {img_file.name}")
 
         # Clean up remaining TIF files unless --keep-tif
         if not self.keep_tif:
@@ -411,17 +411,15 @@ class LidarArchaeoPipeline:
             if len(self.resolutions) > 1:
                 logger.info(f"  --- Résolution {res}m/px ---")
 
-            # For additional resolutions, use suffixed subdirectory and basename
+            # For additional resolutions, use suffixed subdirectory
             if res_suffix:
                 vis_dir = self.vis_dir / f"{basename}{res_suffix}"
-                pdf_basename = f"{basename}{res_suffix}"
             else:
                 vis_dir = self.vis_dir / basename
-                pdf_basename = basename
 
             vis_dir.mkdir(exist_ok=True)
 
-            self.generate_all_visualizations(dtm_path, basename, actual_res)
+            self.generate_all_visualizations(dtm_path, basename, actual_res, vis_dir=vis_dir)
 
         t_total = time.time() - t_start
         logger.info(f"✓ {basename} terminé en {t_total:.1f}s")
@@ -452,29 +450,42 @@ class LidarArchaeoPipeline:
             logger.info(f"Traitement parallèle avec {self.workers} workers...")
             logger.info(f"Fichiers: {len(files)}")
             with ProcessPoolExecutor(max_workers=self.workers) as executor:
+                # Pass resolutions as comma-separated string for multiprocessing serialization
+                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), self.resolution, self.force, self.ground_method, self.force_classify, self.keep_tif, self.quality, self.only_viz, self.skip_viz): laz_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): laz_file
                     for laz_file in files
                 }
                 done = 0
-                for future in as_completed(future_to_file):
-                    laz_file = future_to_file[future]
-                    done += 1
-                    try:
-                        success = future.result()
-                        results[laz_file.name] = success
-                        status = "✓" if success else "✗"
-                        logger.info(f"  [{done}/{len(files)}] {status} {laz_file.name}")
-                    except Exception as e:
-                        logger.error(f"  [{done}/{len(files)}] ✗ {laz_file.name}: {e}")
-                        logger.debug(f"  Traceback:", exc_info=True)
-                        results[laz_file.name] = False
+                try:
+                    for future in as_completed(future_to_file):
+                        laz_file = future_to_file[future]
+                        done += 1
+                        try:
+                            success = future.result()
+                            results[laz_file.name] = success
+                            status = "✓" if success else "✗"
+                            logger.info(f"  [{done}/{len(files)}] {status} {laz_file.name}")
+                        except Exception as e:
+                            logger.error(f"  [{done}/{len(files)}] ✗ {laz_file.name}: {e}")
+                            logger.debug(f"  Traceback:", exc_info=True)
+                            results[laz_file.name] = False
+                except KeyboardInterrupt:
+                    logger.info("Interruption — annulation des travaux en cours...")
+                    for f in future_to_file:
+                        f.cancel()
+                    executor.shutdown(wait=False, cancel_futures=True)
+                    logger.info("Travaux annulés.")
+                    return
         else:
             total = len(files)
             for idx, laz_file in enumerate(files, 1):
                 logger.info(f"--- Fichier {idx}/{total} ---")
                 try:
                     results[laz_file.name] = self.process_file(laz_file)
+                except KeyboardInterrupt:
+                    logger.info("Interruption — arrêt immédiat.")
+                    return
                 except Exception as e:
                     logger.error(f"✗ Erreur traitement {laz_file.name}: {e}")
                     logger.debug("Traceback:", exc_info=True)
@@ -500,7 +511,6 @@ class LidarArchaeoPipeline:
         logger.info(f"\nRésultats dans: {self.output_dir}")
         logger.info(f"  • DTM           : {self.dtm_dir}")
         logger.info(f"  • Visualisations: {self.vis_dir}")
-        logger.info(f"  • Rapports PDF  : {self.pdf_dir}")
 
         # Clean up temporary files
         logger.info("Nettoyage des fichiers temporaires...")
@@ -516,7 +526,7 @@ 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=85, only_viz=None, skip_viz=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'):
     """Standalone function for multiprocessing — creates its own pipeline instance.
 
     Each worker gets its own temp directory to avoid file conflicts.
@@ -537,7 +547,7 @@ def _process_file_standalone(laz_file_str, input_dir, output_dir, resolution, fo
         worker_logger.addHandler(handler)
     worker_logger.addFilter(_file_filter)
 
-    pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif, quality=quality, only_viz=only_viz, skip_viz=skip_viz)
+    pipeline = LidarArchaeoPipeline(input_dir, output_dir, resolution=resolution, workers=1, force=force, ground_method=ground_method, force_classify=force_classify, keep_tif=keep_tif, quality=quality, only_viz=only_viz, skip_viz=skip_viz, output_format=output_format)
     basename = _file_basename(laz_file_str)
     pipeline.temp_dir = pipeline.output_dir / "temp" / basename
     pipeline.temp_dir.mkdir(exist_ok=True)
diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py
index 3279a7d..ab6ebef 100644
--- a/lidar_pipeline/rendering.py
+++ b/lidar_pipeline/rendering.py
@@ -1,8 +1,8 @@
-"""Rendering module: colormap registry, GeoTIFF-to-WebP conversion, and PDF report generation.
+"""Rendering module: colormap registry, GeoTIFF-to-image conversion, and PDF report generation.
 
 Contains:
 - COLORMAPS: registry mapping filename keywords to (cmap, title, legend, description)
-- tif_to_png(): convert a GeoTIFF to a WebP visualization with legend, scale bar, north arrow
+- tif_to_png(): convert a GeoTIFF to a WebP/AVIF visualization with legend, scale bar, north arrow
 - generate_pdf_report(): generate an A3 PDF report with all visualizations
 """
 
@@ -74,18 +74,10 @@ COLORMAPS = {
         'description': 'Détecte les ruptures de pente — utile pour bords de terrasses et levées',
         'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
     },
-    'svf': {
-        'cmap': 'viridis',
-        'title': 'Sky-View Factor (Ray-tracing 16 directions)',
-        'legend': 'Fraction de ciel visible depuis chaque point\nSombre = Encaissé (fossés, vallées, rues)\nClair = Dégagé (sommets, plateformes, plateaux)',
-        'description': 'Ray-tracing sur 16 azimuts — élimine l\'éclairage, détecte structures linéaires et enclos',
-        'vmin_mode': 'percentile', 'vmin_pct': 5,
-        'vmax_mode': 'percentile', 'vmax_pct': 95,
-    },
     'mslrm': {
         'cmap': 'RdBu_r',
-        'title': 'MSRM - Multi-Scale Relief Model (5 échelles)',
-        'legend': 'Relief combiné à 5 échelles (5m à 100m)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve, fossé)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = 5 échelles combinées\nMSRM détecte du micro au macro',
+        'title': 'MSRM - Multi-Scale Relief Model (échelles adaptatives)',
+        'legend': 'Relief combiné multi-échelles (adapté à la résolution)\nRouge = Surélévation (mur, tumulus, levée)\nBleu = Dépression (fossé, douve)\n\nDifférence avec LRM:\nLRM = 1 échelle (15m)\nMSRM = échelles combinées pondérées\nMSRM détecte du micro au macro',
         'description': 'Combine LRM à 5 échelles — détecte structures de 5m à 100m simultanément',
         'vmin_mode': 'symmetric', 'sym_pct': (2, 98),
     },
@@ -114,8 +106,8 @@ COLORMAPS = {
     },
     'tpi': {
         'cmap': 'BrBG',
-        'title': 'TPI - Topographic Position Index (2 échelles)',
-        'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\nCombine échelle fine 5m + large 100m',
+        'title': 'TPI - Topographic Position Index (4 échelles)',
+        'legend': 'Position dans le paysage\nBrun/Sombre = Plus bas que le voisinage (fossé, vallée)\nVert/Clair = Plus haut que le voisinage (crête, plateau)\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),
     },
@@ -128,23 +120,23 @@ COLORMAPS = {
     },
     'roughness': {
         'cmap': 'magma',
-        'title': 'Rugosité de Surface (Écart-type local 5m)',
-        'legend': 'Irrégularité du terrain dans un voisinage de 5m\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nMax: {vmax:.2f}m',
+        'title': 'Rugosité Multi-Échelle (3m + 15m)',
+        'legend': 'Irrégularité du terrain combinée fine + large\nSombre = Surface lisse (route, mur, sol plat)\nClair = Surface rugueuse (végétation, ruines, pierres)\nCombine rugosité fine 3m (70%) + large 15m (30%)',
         'description': 'Mesure la variabilité locale — surfaces anthropiques lisses vs naturelles rugueuses',
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 97,
     },
     'anomalies': {
         'cmap': 'coolwarm',
-        'title': 'Anomalies Statistiques (Z-score x Moran\'s I)',
-        'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine z-score (intensité) et\nMoran\'s I (regroupement spatial)',
+        'title': 'Anomalies Statistiques (MSRM multi-échelle + Moran\'s I)',
+        'legend': 'Anomalies topographiques significatives\nRouge vif = Surélévation anormale (mur, tumulus)\nBleu vif = Dépression anormale (fossé, doline)\nBlanc/gris = Normal\n\nCombine MSRM normalisé (intensité) et\nMoran\'s I (regroupement spatial)',
         'description': 'Détecte uniquement les anomalies statistiquement significatives — filtre le bruit de fond',
         'vmin_mode': 'symmetric', 'sym_pct': (5, 95),
     },
     'wavelet': {
         'cmap': 'cividis',
         'title': 'Ondelette Mexican Hat (CWT multi-échelle)',
-        'legend': 'Réponse de la transformée en ondelette à 5 échelles\nÉchelles: 2m, 5m, 10m, 20m, 50m\n\nClair = Structure détectée à cette échelle\nSombre = Pas de structure\n\nOptimisé pour formes circulaires:\ntumulus, enclos, fossés annulaires',
+        '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),
     },
@@ -156,14 +148,6 @@ COLORMAPS = {
         'vmin_mode': 'fixed', 'vmin_val': 0,
         'vmax_mode': 'percentile', 'vmax_pct': 98,
     },
-    'local_dominance': {
-        'cmap': 'RdYlBu_r',
-        'title': 'Dominance Locale (position relative dans le voisinage)',
-        'legend': 'Proportion du voisinage sous le point central\nRouge = Point dominant (sommet, crête)\nBleu = Point encaissé (fossé, vallée)\nRayon: 15m',
-        'description': 'Mesure la saillie locale — complémentaire de l\'openness',
-        'vmin_mode': 'percentile', 'vmin_pct': 2,
-        'vmax_mode': 'percentile', 'vmax_pct': 98,
-    },
 }
 
 # RGB entries (ortho/topo) are handled specially
@@ -271,24 +255,26 @@ def _nice_scale(extent_m):
     return chosen, f"{chosen} m"
 
 
-def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, quality=85):
-    """Convert GeoTIFF to visualization WebP with GPS coordinates, legend, and scale bar.
+def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None, quality=98, output_format='avif'):
+    """Convert GeoTIFF to visualization image (WebP or AVIF) with GPS coordinates, legend, and scale bar.
 
     Args:
         tif_file: Path to input GeoTIFF.
-        vis_dir: Output directory for the WebP file.
+        vis_dir: Output directory for the image file.
         resolution: Grid resolution in m/px.
         keep_tif: If True, keep the source TIFF after conversion.
         source_info: Dict with method/date/basename for metadata.
-        quality: WebP quality (1-100). Use 100 for lossless. Default 85.
+        quality: Image quality (1-100). Use 100 for lossless. Default 95.
+        output_format: Output format ('webp' or 'avif'). Default 'webp'.
 
     Returns:
-        Path to output WebP file, or None on failure.
+        Path to output image file, or None on failure.
     """
     if not tif_file or not tif_file.exists():
         return None
 
-    webp_file = vis_dir / f"{tif_file.stem}.webp"
+    ext = 'avif' if output_format == 'avif' else 'webp'
+    output_file = vis_dir / f"{tif_file.stem}.{ext}"
 
     try:
         with rasterio.open(tif_file) as src:
@@ -582,7 +568,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
 
             fig.patch.set_facecolor('white')
 
-            # Save as PNG then convert to WebP — fixed layout, no bbox_inches='tight'
+            # Save as PNG then convert to final format — fixed layout, no bbox_inches='tight'
             save_dpi = 200 if width > 3000 else 150
             png_temp = vis_dir / f"{tif_file.stem}_temp.png"
             try:
@@ -591,20 +577,21 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
                 plt.close()
 
             img = PILImage.open(str(png_temp))
+            pil_format = 'AVIF' if output_format == 'avif' else 'WEBP'
             if quality >= 100:
-                img.save(str(webp_file), format='WEBP', lossless=True)
+                img.save(str(output_file), format=pil_format, lossless=True)
             else:
-                img.save(str(webp_file), format='WEBP', quality=quality)
+                img.save(str(output_file), format=pil_format, quality=quality)
             png_temp.unlink(missing_ok=True)
 
             # Delete source TIFF (unless --keep-tif)
             if not keep_tif:
                 tif_file.unlink(missing_ok=True)
 
-            return webp_file
+            return output_file
 
     except Exception as e:
-        logger.error(f"    Erreur conversion WebP: {e}", exc_info=True)
+        logger.error(f"    Erreur conversion {ext.upper()}: {e}", exc_info=True)
         return None
 
 
diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py
index b6088a6..63e1040 100644
--- a/lidar_pipeline/visualizations.py
+++ b/lidar_pipeline/visualizations.py
@@ -250,7 +250,7 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
 def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
     """Generate multi-directional hillshade with contrast enhancement — GPU if available.
 
-    Combines 4-direction hillshade (NW, NE, SW, SE) with slope shading.
+    Combines 8-direction hillshade with slope shading for balanced illumination.
     Applies percentile normalization and gamma correction to restore
     contrast lost by averaging multiple azimuths.
     """
@@ -279,8 +279,9 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
             sin_slope = xp.sin(slope)
             cos_slope = xp.cos(slope)
 
-        azimuts = [315, 45, 225, 135]
-        altitude = 30
+        # 8 azimuths for balanced illumination (eliminates directional bias)
+        azimuts = [0, 45, 90, 135, 180, 225, 270, 315]
+        altitude = 35  # Higher altitude for better micro-relief detection
         hillshades = []
 
         alt_rad = xp.radians(xp.array(altitude))
@@ -297,7 +298,6 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
         combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
 
         # Contrast enhancement: percentile stretch + gamma
-        # Averaging 4 azimuths flattens contrast — this restores it
         combined_np = to_cpu(combined)
         nan_mask = shared.nan_mask if shared else np.isnan(to_cpu(dem_np) if HAS_GPU else dem_np)
         valid = combined_np[~nan_mask]
@@ -415,7 +415,11 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None):
 # ============================================================
 
 def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None):
-    """Local Relief Model - deviation from local mean (GPU if available)."""
+    """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 HAS_GPU else ""
     logger.info(f"  → Local Relief Model{gpu_tag}...")
     t0 = time.time()
@@ -429,7 +433,10 @@ def generate_lrm(dem_file, basename, vis_dir, resolution, shared=None):
         else:
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
-            local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15/resolution)
+            # 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)
@@ -473,7 +480,7 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
         angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
         dx_dir = np.cos(angles)
         dy_dir = np.sin(angles)
-        max_dist = int(50 / res)
+        max_dist = int(100 / res)
 
         padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
         svf = xp.zeros_like(dem)
@@ -520,6 +527,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
     - 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.
     """
     name = "positive_openness" if positive else "negative_openness"
     gpu_tag = " [GPU]" if HAS_GPU else ""
@@ -548,7 +556,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
         angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
         dx_dir = np.cos(angles)
         dy_dir = np.sin(angles)
-        max_dist = int(50 / res)
+        max_dist = int(100 / res)
 
         padded = xp.pad(dem, max_dist, mode='constant', constant_values=xp.nan)
         openness_sum = xp.zeros_like(dem)
@@ -646,7 +654,11 @@ def generate_local_dominance(dem_file, basename, vis_dir, resolution, shared=Non
 
 
 def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
-    """Multi-Scale Relief Model (MSRM) - LRM at 5 scales combined (GPU if available)."""
+    """Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).
+
+    Scales adapt to resolution. Std normalization per scale.
+    Weighted combination favoring archaeologically relevant scales (5-25m).
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Multi-Scale Relief Model (MSRM){gpu_tag}...")
     t0 = time.time()
@@ -662,7 +674,18 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
 
-        sigmas = [5, 10, 25, 50, 100]
+        # Adaptive scales: finer at higher resolution
+        min_scale = max(2.0, resolution * 4)
+        candidate_scales = [2, 5, 10, 20, 50, 100, 200]
+        sigmas = [s for s in candidate_scales if s >= min_scale]
+
+        # Archaeological weights: favor 5-25m range (ditches, enclosures, tumulus)
+        scale_weights = {
+            2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6, 200: 0.4,
+        }
+        weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
+
+        logger.info(f"    MSRM échelles: {sigmas}m")
         lrm_stack = []
 
         for sigma in sigmas:
@@ -673,16 +696,19 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
                 local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
             lrm = dem_np - local_mean
             lrm[nan_mask] = np.nan
+            # Std normalization: x / std — preserves sign and contrast better than z-score
             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))
+            lrm = lrm / lrm_std
+            lrm_stack.append(lrm.astype(np.float32))
 
+        # Weighted combination
         lrm_array = np.array(lrm_stack)
+        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')
-                mslrm = np.sqrt(np.nanmean(lrm_array ** 2, axis=0))
+                mslrm = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
         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}")
@@ -696,7 +722,8 @@ 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 fine (5m) and broad (100m) scales.
+    Computed at 4 scales with std normalization and weighted combination.
+    Weights favor fine and medium scales (archaeologically relevant).
     """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → TPI multi-échelle{gpu_tag}...")
@@ -713,29 +740,34 @@ def generate_tpi(dem_file, basename, vis_dir, resolution, shared=None):
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
 
-        fine_size = int(5 / resolution)
-        if fine_size % 2 == 0:
-            fine_size += 1
-        if shared:
-            fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
-        else:
-            fine_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=fine_size)
-        tpi_fine = dem_np - fine_mean
-        tpi_fine[nan_mask] = np.nan
+        # 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)
 
-        broad_size = int(100 / resolution)
-        if broad_size % 2 == 0:
-            broad_size += 1
-        if shared:
-            broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size)
-        else:
-            broad_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=broad_size)
-        tpi_broad = dem_np - broad_mean
-        tpi_broad[nan_mask] = np.nan
+        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))
 
-        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)
+        # 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)
@@ -756,7 +788,7 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
     """SAILORE - Self-Adaptive Improved Local Relief Model (GPU if available).
 
     Kernel size adapts to local slope: flat areas get larger kernels,
-    steep areas get smaller kernels.
+    steep areas get smaller kernels. Scales adapt to resolution.
     """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → SAILORE (LRM adaptatif){gpu_tag}...")
@@ -778,8 +810,13 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
             slope_deg = np.degrees(slope)
             slope_deg[nan_mask] = np.nan
 
-        sigma_min = 2.0 / resolution
-        sigma_max = 25.0 / resolution
+        # Adaptive scales: finer at higher resolution
+        sigma_min_m = max(1.0, 2.0 * 0.5 / resolution)  # 2m at 0.5, ~5m at 0.2
+        sigma_mid_m = max(5.0, 13.5 * 0.5 / resolution)  # 13.5m at 0.5, ~33m at 0.2
+        sigma_max_m = max(5.0, 25.0 * 0.5 / resolution)  # 25m at 0.5, ~62m at 0.2
+        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)
 
         if shared:
@@ -822,7 +859,11 @@ def generate_sailore(dem_file, basename, vis_dir, resolution, shared=None):
 # ============================================================
 
 def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
-    """Surface roughness - standard deviation of elevation in a window (GPU-accelerated)."""
+    """Surface roughness - multi-scale standard deviation (GPU-accelerated).
+
+    Combines fine (3m) and broad (15m) roughness for better detection
+    of archaeological features at multiple scales.
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Rugosité de surface{gpu_tag}...")
     t0 = time.time()
@@ -838,20 +879,43 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
 
-        window_size = int(5 / resolution)
-        if window_size % 2 == 0:
-            window_size += 1
+        # Fine roughness (3m window)
+        fine_size = max(3, int(3 / resolution))
+        if fine_size % 2 == 0:
+            fine_size += 1
 
-        # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (NaN-aware)
         if shared:
-            local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window_size)
-            # For local_mean_sq, we need to filter filled², not filled
-            local_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=window_size)
-            local_mean_sq[shared.nan_mask] = np.nan
+            fine_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=fine_size)
+            fine_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
+            fine_mean_sq[shared.nan_mask] = np.nan
         else:
-            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))
+            fine_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=fine_size)
+            fine_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=fine_size)
+        roughness_fine = np.sqrt(np.maximum(fine_mean_sq - fine_mean * fine_mean, 0))
+        roughness_fine[nan_mask] = np.nan
+
+        # Broad roughness (15m window)
+        broad_size = max(3, int(15 / resolution))
+        if broad_size % 2 == 0:
+            broad_size += 1
+
+        if shared:
+            broad_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=broad_size)
+            broad_mean_sq = _filter_nanaware(shared.filled.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
+            broad_mean_sq[shared.nan_mask] = np.nan
+        else:
+            broad_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=broad_size)
+            broad_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=broad_size)
+        roughness_broad = np.sqrt(np.maximum(broad_mean_sq - broad_mean * broad_mean, 0))
+        roughness_broad[nan_mask] = np.nan
+
+        # Std normalization per scale then weighted combination
+        fine_valid = roughness_fine[~nan_mask]
+        broad_valid = roughness_broad[~nan_mask]
+        fine_std = max(np.nanstd(fine_valid), 0.01) if len(fine_valid) > 0 else 0.01
+        broad_std = max(np.nanstd(broad_valid), 0.01) if len(broad_valid) > 0 else 0.01
+
+        roughness = 0.7 * roughness_fine / fine_std + 0.3 * roughness_broad / broad_std
         roughness[nan_mask] = np.nan
 
         roughness = to_cpu(roughness)
@@ -868,7 +932,11 @@ def generate_roughness(dem_file, basename, vis_dir, resolution, shared=None):
 # ============================================================
 
 def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
-    """Statistical anomaly detection - z-score of local relief + Local Moran's I — GPU if available."""
+    """Statistical anomaly detection - std-normalized multi-scale relief + Local Moran's I — GPU if available.
+
+    Uses MSRM (multi-scale LRM) instead of single-scale LRM for better detection
+    of anomalies at all scales.
+    """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Détection anomalies statistiques{gpu_tag}...")
     t0 = time.time()
@@ -880,30 +948,60 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
             crs = shared.crs
             dem_np = shared.dem_np
             nan_mask = shared.nan_mask
-            lrm = shared.lrm_15.copy()
         else:
             dem_np, transform, crs = _read_dem(dem_file)
             nan_mask = np.isnan(dem_np)
-            lrm_mean_val = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15 / resolution)
-            lrm = dem_np - lrm_mean_val
+
+        # Multi-scale LRM: compute MSRM-like combined relief
+        min_scale = max(2.0, resolution * 4)
+        candidate_scales = [2, 5, 10, 20, 50, 100]
+        sigmas = [s for s in candidate_scales if s >= min_scale]
+        lrm_stack = []
+
+        for sigma in sigmas:
+            sigma_px = sigma / resolution
+            if shared:
+                local_mean = _filter_nanaware_from_filled(shared, xp_gaussian_filter, sigma=sigma_px)
+            else:
+                local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
+            lrm = dem_np - local_mean
             lrm[nan_mask] = np.nan
+            # Std normalization — preserves contrast better than z-score
+            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
+            else:
+                lrm_norm = lrm
+            lrm_stack.append(lrm_norm.astype(np.float32))
 
-        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
+        # Weighted RMS combination (favor 5-25m scales)
+        scale_weights = {2: 0.8, 5: 2.0, 10: 1.8, 20: 1.5, 50: 1.0, 100: 0.6}
+        weights = np.array([scale_weights.get(s, 1.0) for s in sigmas])
+        lrm_array = np.array(lrm_stack)
+        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')
+                msrm = np.sqrt(np.nansum((lrm_array ** 2) * weights_3d, axis=0) / np.sum(weights))
+        msrm[nan_mask] = np.nan
 
-        window = int(10 / resolution)
+        # Std normalization of MSRM — preserves contrast better than z-score
+        valid_msrm = msrm[~nan_mask]
+        msrm_std = max(np.nanstd(valid_msrm), 0.01) if len(valid_msrm) > 0 else 0.01
+        z_score = msrm / msrm_std
+
+        # Local Moran's I for spatial clustering
+        window = max(3, int(10 / resolution))
         if window % 2 == 0:
             window += 1
 
         if shared:
-            local_mean = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window)
+            local_mean_z = _filter_nanaware_from_filled(shared, xp_uniform_filter, size=window)
         else:
-            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
+            local_mean_z = _filter_nanaware(z_score, xp_uniform_filter, size=window)
+        z_mean_global = np.nanmean(z_score[~nan_mask]) if np.any(~nan_mask) else 0
+        z_std_global = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
+        morans_i = z_score * (local_mean_z - z_mean_global) / z_std_global
         anomaly_score = np.abs(z_score) * np.sign(morans_i)
         anomaly_score[nan_mask] = np.nan
 
@@ -922,7 +1020,13 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution, shared=None):
 def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
     """Mexican Hat wavelet multi-scale analysis (GPU if available).
 
-    CWT 2D at multiple scales to detect circular features.
+    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).
     """
     gpu_tag = " [GPU]" if HAS_GPU else ""
     logger.info(f"  → Ondelette Mexican Hat multi-échelle{gpu_tag}...")
@@ -941,7 +1045,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))
 
-        scales = [2, 5, 10, 20, 50]
+        # Adapt scales to resolution: finer scales available at higher resolution
+        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]
+
+        # 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
+        }
+        weights = np.array([scale_weights.get(s, 1.0) for s in scales])
+
+        logger.info(f"    Échelles CWT: {scales}m (résolution {resolution}m/px)")
         wavelet_stack = []
 
         for scale_m in scales:
@@ -954,15 +1076,21 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution, shared=None):
                 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
+        stack = np.array(wavelet_stack)
+        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.nanmean(np.array(wavelet_stack) ** 2, axis=0))
+                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)
diff --git a/run.sh b/run.sh
index ec1d551..9dec617 100755
--- a/run.sh
+++ b/run.sh
@@ -3,18 +3,21 @@
 # Utilisation: ./run.sh [options]
 #
 # Options:
-#   -r RESOLUTION   Résolution en m/px (défaut: 0.5)
+#   -r RESOLUTION   Résolution en m/px, ou multiples séparées par virgules (défaut: 0.5, ex: 0.5,0.2)
 #   -w WORKERS      Nombre de workers parallèles (défaut: 1)
 #   -g              Activer l'accélération GPU
-#   -v              Mode verbeux (timestamps + niveaux)
-#   --debug         Mode debug (détails internes fichier:ligne)
-#   -f / --force    Régénérer tous les fichiers même si existants
-#   --keep-tif       Conserver les fichiers TIFF pour régénérer les WebP
+#   -v              Mode verbeux
+#   --debug         Mode debug
+#   -f / --force    Régénérer tous les fichiers
+#   --keep-tif      Conserver les fichiers TIFF
 #   --force-classification
-#                   Reclassifier le sol même si le fichier .las existe déjà
 #   --ground-classification {auto,smrf,csf}
-#                   Méthode de classification du sol (défaut: auto)
-#   --file NOM...   Traiter un ou plusieurs fichiers LAZ spécifiques
+#   --quality N     Qualité image 1-100 (défaut: 98)
+#   --lossless      Compression lossless
+#   --format FMT    Format de sortie : avif (défaut) ou webp
+#   --only VIZ...   Générer uniquement ces visualisations
+#   --skip VIZ...   Exclure ces visualisations
+#   --file NOM...   Traiter un ou plusieurs fichiers LAZ
 #   --test          Exécuter les tests unitaires
 #   -h              Afficher l'aide complète
 
@@ -32,7 +35,7 @@ if [ $# -eq 0 ]; then
     echo "Usage: $0 [options]"
     echo ""
     echo "Options:"
-    echo "  -r RESOLUTION   Résolution en m/px (défaut: 0.5)"
+    echo "  -r RESOLUTION   Résolution en m/px, ou multiples (défaut: 0.5, ex: 0.5,0.2)"
     echo "  -w WORKERS      Nombre de workers CPU parallèles (défaut: 1)"
     echo "  -g              Activer l'accélération GPU NVIDIA"
     echo "  -v              Mode verbeux (timestamps + niveaux)"
@@ -52,6 +55,7 @@ if [ $# -eq 0 ]; then
     echo "  $0 -g -w 4                               # GPU + 4 workers"
     echo "  $0 -g -v                                 # GPU + verbeux"
     echo "  $0 -g -r 0.2                             # Haute résolution"
+    echo "  $0 -g -r 0.5,0.2                          # Multi-résolution (0.5m + 0.2m)"
     echo "  $0 -g --force                            # Régénérer WebP (DTM conservé si --keep-tif)"
     echo "  $0 -g --force-classification             # Reclassifier le sol seulement"
     echo "  $0 -g --ground-classification csf        # Forcer CSF (terrain complexe)"
@@ -69,6 +73,7 @@ GROUND_METHOD=""
 FORCE_CLASSIFY_FLAG=""
 KEEP_TIF_FLAG=""
 QUALITY=""
+FORMAT_FLAG=""
 ONLY_FLAG=""
 SKIP_FLAG=""
 
@@ -88,6 +93,7 @@ while [ $# -gt 0 ]; do
         --ground-classification=*) GROUND_METHOD="${1#--ground-classification=}"; shift ;;
         --quality) QUALITY="--quality $2"; shift 2 ;;
         --lossless) QUALITY="--lossless"; shift ;;
+        --format) FORMAT_FLAG="--format $2"; shift 2 ;;
         --only) shift; ONLY_FLAG="--only"; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do ONLY_FLAG="$ONLY_FLAG $1"; shift; done ;;
         --skip) shift; SKIP_FLAG="--skip"; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do SKIP_FLAG="$SKIP_FLAG $1"; shift; done ;;
         --file) shift; while [ $# -gt 0 ] && [[ ! "$1" =~ ^- ]]; do FILE_ARGS="$FILE_ARGS $1"; shift; done ;;
@@ -109,8 +115,9 @@ while [ $# -gt 0 ]; do
             echo "  --keep-tif      Conserver les TIFF pour régénérer les WebP"
             echo "  --ground-classification {auto,smrf,csf}"
             echo "                  Méthode de classification du sol (défaut: auto)"
-            echo "  --quality N     Qualité WebP 1-100 (défaut: 85, 100=lossless)"
-            echo "  --lossless      Compression WebP lossless (équivalent à --quality 100)"
+            echo "  --quality N     Qualité image 1-100 (défaut: 98, 100=lossless)"
+            echo "  --lossless      Compression lossless (équivalent à --quality 100)"
+            echo "  --format FMT    Format de sortie : avif (défaut) ou webp"
             echo "  --only VIZ...   Générer uniquement ces visualisations"
             echo "  --skip VIZ...   Exclure ces visualisations"
             echo "  --file NOM...   Traiter un ou plusieurs fichiers LAZ"
@@ -118,14 +125,15 @@ while [ $# -gt 0 ]; do
             echo "  -h              Afficher cette aide"
             echo ""
             echo "Visualisations disponibles:"
-            echo "  hillshade slope aspect curvature svf lrm pos_open neg_open"
-            echo "  mslrm tpi sailore roughness anomalies wavelet flow local_dominance ortho topo"
+            echo "  hillshade slope aspect curvature lrm pos_open neg_open"
+            echo "  mslrm tpi sailore roughness anomalies wavelet flow ortho topo"
             echo ""
             echo "Exemples:"
             echo "  $0 -g                                    # GPU, auto"
             echo "  $0 -g -w 4                               # GPU + 4 workers"
             echo "  $0 -g -v                                 # GPU + verbeux"
             echo "  $0 -g -r 0.2                             # Haute résolution"
+    echo "  $0 -g -r 0.5,0.2                          # Multi-résolution (0.5m + 0.2m)"
             echo "  $0 -g --force                            # Régénérer WebP"
             echo "  $0 -g --only hillshade svf lrm           # Seulement 3 visualisations"
             echo "  $0 -g --skip ortho topo                  # Sans les overlays IGN"
@@ -178,7 +186,8 @@ echo "  Verbeux    : $([ -n "$VERBOSE_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Force      : $([ -n "$FORCE_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Force classif.: $([ -n "$FORCE_CLASSIFY_FLAG" ] && echo 'OUI' || echo 'non')"
 echo "  Keep TIFF    : $([ -n "$KEEP_TIF_FLAG" ] && echo 'OUI' || echo 'non')"
-echo "  Qualité WebP : $([ -n "$QUALITY" ] && echo "$QUALITY" || echo '85')"
+echo "  Qualité image : $([ -n "$QUALITY" ] && echo "$QUALITY" || echo '98')"
+echo "  Format      : $([ -n "$FORMAT_FLAG" ] && echo "${FORMAT_FLAG#--format }" || echo 'avif')"
 echo "  Classification sol : $([ -n "$GROUND_METHOD" ] && echo "$GROUND_METHOD" || echo 'auto')"
 if [ -n "$ONLY_FLAG" ]; then
     echo "  Visualisations: uniquement${ONLY_FLAG#--only}"
@@ -190,7 +199,7 @@ if [ -n "$FILE_ARGS" ]; then
 fi
 echo "============================================"
 
-CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG $QUALITY"
+CMD_ARGS="-o /data/output -r $RESOLUTION -w $WORKERS $VERBOSE_FLAG $FORCE_FLAG $FORCE_CLASSIFY_FLAG $KEEP_TIF_FLAG $QUALITY $FORMAT_FLAG"
 if [ -n "$GROUND_METHOD" ]; then
     CMD_ARGS="$CMD_ARGS --ground-classification $GROUND_METHOD"
 fi