diff --git a/AGENTS.md b/AGENTS.md
index 2be9c2a..269671b 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -25,14 +25,16 @@
 - **Sub-tuilage intégral** : `_CARTO_SUBTILED_VIZ` (vide dans `index.py`) découpe TOUTES les couches en quadrants 500 m à 0,2 m ; ortho/topo sont encodées en AVIF q75 (`_SUBTILE_DETAIL_VIZ`) contre q55 pour les rampes de couleur. Une couche qui échoue à la découpe retombe en dalle entière (`_fallback_full_dalle`) sans pénaliser les autres.
 - **Tuiles XYZ réutilisables hors du projet** : `/tiles/{layer}/{z}/{x}/{y}.png` suit le schéma OpenStreetMap (256 px, EPSG:3857, PNG RGBA transparent hors emprise, CORS `*`) — JOSM, iD, QGIS, uMap et MapLibre le consomment tel quel, avec découverte via TileJSON / WMTS / `josm.imagery.xml`. Le 512 px (`@2x.webp`) est réservé à l'interface interne (moitié moins de requêtes en HTTP/1.1). Toute évolution du gabarit d'URL casse des configurations clientes : la changer demande une décision explicite.
 - **Bilingual naming**: all code identifiers are English; every user-facing string, log message, argparse help, and comment is French.
-- **Adding a visualization requires 4 edits**: (1) `generate_X()` in `visualizations.py`, (2) entry in `VIZ_STEPS` in `pipeline.py`, (3) entry in `COLORMAPS` in `rendering.py`, (4) entry in `VIZ_LEGENDS` in `index.py` (title/legend/description + sampled cmap gradient — single text source merged into `COLORMAPS` at import, also served in `/api/map/meta` and le TileJSON). Missing any one breaks the pipeline.
+- **Adding a visualization requires 4 edits**: (1) `generate_X()` in `visualizations.py`, (2) entry in `VIZ_STEPS` in `pipeline.py`, (3) entry in `COLORMAPS` in `rendering.py` (or `RGB_LEGENDS` for an RGB output), (4) entry in `VIZ_LEGENDS` in `index.py` (title/legend/description + sampled cmap gradient — single text source merged into `COLORMAPS` at import, also served in `/api/map/meta` and le TileJSON). Missing any one breaks the pipeline.
 - **`generate_*` signature is strict**: `(dem_file, basename, vis_dir, resolution, shared=None)` returning `Path` on success, `None` on failure. IGN overlays (`ortho`, `topo`) omit `shared`.
 - **Return `None` on failure, never raise**: `dtm.py`, `visualizations.py`, and `ign.py` all return `None` to let the pipeline continue. Raising aborts the entire file.
 - **Logger is always `logging.getLogger("lidar")`**, never `__name__`. All modules route through this single logger so worker processes can configure it.
+- **Openness à échelle fixe** : `generate_openness` normalise par des références figées `OPENNESS_POS_REF` / `OPENNESS_NEG_REF` (degrés, médianes mesurées sur 15 dalles réparties sur le territoire) et plus par z-score de dalle — même ouverture = même couleur, mosaïque jointive.
 - **Filename special-cases** in `_expected_output_path()`: `pos_open` → `positive_openness`, `neg_open` → `negative_openness`, `hillshade` → `hillshade_multi`.
 - **Default output is AVIF**, not WebP. Use `--format webp` for WebP. Quality default is 60 (visually lossless on smooth color ramps, ~÷3 vs q98).
 - **Calage vertical des faisceaux de vol** : chaque tuile mélange plusieurs passes (1-2 `PointSourceId` par passe) parfois biaisées verticalement de quelques cm (±2,5 cm mesurés sur 1000_6882). `create_dtm_fast` mesure l'offset robuste de chaque faisceau (points sol, maille 1 m, surface médiane itérée 3×) et retranche les offsets ≥ 0,5 cm (`STRIP_ALIGN_THRESHOLD` dans `dtm.py`) avant rastérisation. Offsets calculés **par tuile** (ils dérivent le long d'une ligne de vol : jamais de table globale), mémoïsés par LAS sol, consignés dans `DTM/*_dtm*_stripalign.json` (sidecar de cache : absent, ou version/seuil/paramètres différents ⇒ régénération du DTM). Désactivable : `--no-strip-align`.
 - **Gigue intra-faisceau (2ᵉ passe du calage)** : les lignes de balayage successives d'une MÊME passe peuvent être décalées verticalement de façon aléatoire (vibration capteur / bruit haute fréquence de trajectoire) — un offset constant par faisceau n'y suffit pas. `_strip_jitter_offsets` découpe chaque faisceau en fenêtres de temps GPS (`STRIP_JITTER_BIN` = 0,1 s, origine de temps propre à chaque faisceau), mesure l'offset robuste de chaque fenêtre contre la surface médiane des AUTRES faisceaux (maille 1 m partagée, ≥ `STRIP_JITTER_MIN_CELLS` = 40 cellules), lisse la série (médiane glissante `STRIP_JITTER_SMOOTH` = 5 fenêtres), la borne à ± `STRIP_JITTER_MAX` (10 cm) puis l'interpole au temps GPS de chaque point (`_apply_strip_jitter`) ; en recouvrement à deux faisceaux, chacun reçoit une série (chacun absorbe sa part). Requiert la dimension `gps_time` (silencieusement ignorée sinon). Sidecar version 2 (séries dans `jitter`), couverte par `--no-strip-align`.
+- **Relief orienté (`relief_oriente`, couche par défaut de la carte)** : image RGB unique (GeoTIFF uint8 3 bandes, rendue telle quelle comme ortho/topo via `RGB_KEYWORDS` dans `rendering.py`) — clarté CIELAB = openness positive **locale** (MNT − gaussienne `RELIEF_DETREND_M` = 10 m, rayons `RELIEF_RADII_M` = 5/10/20 m, 16 directions) 65 % + ombrage 35 % ; teinte = aspect, chroma fixe (`RELIEF_CHROMA`). Échelle log fixe `RELIEF_OPEN_RANGE` (pas de statistique par dalle) et support total 40 m < bande de raccord 100 m : dalles jointives. Rapide : détendance + rayons sur grille décimée ~0,8 m (`RELIEF_GRID_M`), noyau dédié qui n'accumule que la moyenne des angles (`_mean_horizon_*` : CuPy `RawKernel` sur GPU, numba parallèle sur CPU, numpy en repli), colorisation fusionnée (numba) ou vectorisée sans trigonométrie (CuPy) via une table L* × teinte (`_relief_lut`). ~5 s de calcul hors préparation sur CPU 12 cœurs. Tout changement de constante change le rendu : régénérer les dalles (`--only relief_oriente --force`).
 - **Openness sous-échantillonnée** : `generate_openness` calcule le lancé de rayons (l'étape la plus coûteuse : 532 s/tuile à 0,2 m sur CPU) sur une grille décimée par blocs (`OPENNESS_DOWNSAMPLE = 2` : max par bloc en positive, min en négative — préserve les reliefs qui bornent l'horizon) puis rééchantillonne en bilinéaire. Coût ÷ facteur³ : 532 s → 40 s (×13). Signal archéologique préservé (corr. 0,93 après lissage) ; la texture de bruit sub-métrique disparaît. `--openness-downsample 1` = pleine résolution. SVF et openness anisotrope ne sont PAS concernées.
 - **Raccord des bords entre tuiles** : les rendus à grand noyau (openness/SVF : rayons 100 m ; LRM : 15 m) tronquent leur fenêtre au bord de dalle — bandes d'artefacts à chaque changement de tuile. `--edge-buffer N` (défaut 0 = off ; case « Raccord des bords » de la carte, `EDGE_BUFFER_METERS` = 100 m dans `mapserve.py`) fait rastériser le DTM sur la **dalle nominale 1 km alignée sur la grille** plus une bande de N m remplie avec les points sol des 8 LAZ voisines (`_neighbor_ground_points` dans `dtm.py` : PDAL en flux, découpe + filtre de classes IGN ; voisine absente = téléchargement automatique depuis le catalogue IGN avant le run, **isolée dans `input/edge_neighbors/`** pour ne pas gonfler le corpus des passes globales (dédupliqué sur tout le lot, `_fetch_edge_neighbors` dans `pipeline.py`) ; introuvable ou échec = bande vide). Les visualisations calculent sur l'emprise étendue puis `rendering.py` (`_core_tile_window`, via `tif_to_crop`/`tif_to_png`) recadre les sorties sur la dalle 1 km exacte lue dans le nom LHD — les AVIF restent des carrés 1 km alignés dans la mosaïque. Tampon consigné dans le tag GeoTIFF `LIDAR_EDGE_BUFFER` du DTM : changer `--edge-buffer` invalide le cache DTM automatiquement (tag absent = 0). Bandes voisines non calées par faisceaux (contexte seul, recadrée hors image finale). Coût : ~7 s de lecture par voisine + ~44 % de pixels en plus à 100 m/0,2 m. Nom hors pattern LHD : option ignorée (bornes d'en-tête, pas de recadrage).
 - **Tests use lazy imports inside each test function**, never at module top, to avoid importing CuPy/GDAL at import time.
diff --git a/README.md b/README.md
index 3f63d0f..bbe7804 100644
--- a/README.md
+++ b/README.md
@@ -4,6 +4,9 @@ Workflow automatisé pour générer des visualisations exploitables à partir de
 
 ## Visualisations (18 par fichier)
 
+### Relief orienté (couche par défaut)
+Une seule image fusionne le micro-relief et l'orientation des pentes : la **clarté** porte le relief local (openness positive sur MNT détendancé, rayons 5–20 m, plus un léger ombrage), la **teinte** porte l'orientation (aspect, cercle CIELAB à clarté constante : aucune couleur ne crée de faux relief). Échelle fixe : les dalles voisines se raccordent sans couture.
+
 ### Visualisations principales
 | # | Visualisation | Utilité archéologique |
 |---|--------------|----------------------|
diff --git a/lidar_pipeline/gpu.py b/lidar_pipeline/gpu.py
index 85e3a7d..f3dc6b4 100644
--- a/lidar_pipeline/gpu.py
+++ b/lidar_pipeline/gpu.py
@@ -399,6 +399,18 @@ def xp_uniform_filter(arr, size):
     return ndimage.uniform_filter(arr, size)
 
 
+def xp_zoom(arr, factor, order=1):
+    """Agrandissement aligné sur les centres de pixels (grid_mode) : chaque
+    pixel source couvre exactement factor×factor pixels, sans décalage."""
+    if _cp is not None and isinstance(arr, _cp.ndarray):
+        try:
+            return _cp_ndimage.zoom(arr, factor, order=order, mode='nearest', grid_mode=True)
+        except Exception as e:
+            logger.warning(f"Zoom GPU échoué ({e}) — repli CPU")
+            arr = to_cpu(arr)
+    return ndimage.zoom(arr, factor, order=order, mode='nearest', grid_mode=True)
+
+
 def xp_minimum_filter(arr, footprint=None, size=None):
     if _cp is not None and isinstance(arr, _cp.ndarray):
         try:
diff --git a/lidar_pipeline/index.py b/lidar_pipeline/index.py
index 3ba8e12..3e49404 100644
--- a/lidar_pipeline/index.py
+++ b/lidar_pipeline/index.py
@@ -45,6 +45,7 @@ VIZ_LABELS = {
     'flow_acc':          'Accumulation d\'écoulement',
     'solar':             'Éclairage solaire',
     'anomaly':           'Carte d\'anomalies',
+    'relief_oriente':    'Relief orienté',
     'ortho':             'Orthophoto IGN',
     'topo':              'Carte topographique IGN',
 }
@@ -106,7 +107,7 @@ VIZ_LEGENDS = {
     },
     'positive_openness': {
         'title': 'Openness Positive (ouverture vers le haut)',
-        'legend': 'Angle d\'ouverture vers le ciel (z-score local)\nClair = Vue dégagée vers le ciel (sommets, plateaux)\nSombre = Vue bloquée (vallées encaissées)\nÉchelle fixe ±3σ — couleurs homogènes entre tuiles',
+        'legend': 'Angle d\'ouverture vers le ciel (écart à une référence nationale figée)\nClair = Vue dégagée vers le ciel (sommets, plateaux)\nSombre = Vue bloquée (vallées encaissées)\nMême angle = même couleur sur toutes les tuiles',
         'description': 'Ray-tracing 8 directions, multi-rayon — détecte crêtes et sommets',
         'cmap': 'YlOrBr',
         'gradient': ('#ffffe5', '#fff7bc', '#fee390', '#fec34f', '#fe9829',
@@ -115,7 +116,7 @@ VIZ_LEGENDS = {
     },
     'negative_openness': {
         'title': 'Openness Negative (ouverture vers le bas)',
-        'legend': 'Angle d\'ouverture vers le bas (z-score local)\nClair = Surplomb (bords de fossé, grottes)\nSombre = Terrain plat (fonds de vallée)\nÉchelle fixe ±3σ — couleurs homogènes entre tuiles\nMeilleur détecteur de cavités et dolines',
+        'legend': 'Angle d\'ouverture vers le bas (écart à une référence nationale figée)\nClair = Surplomb (bords de fossé, grottes)\nSombre = Terrain plat (fonds de vallée)\nMême angle = même couleur sur toutes les tuiles\nMeilleur détecteur de cavités et dolines',
         'description': 'Ray-tracing 8 directions, multi-rayon — détecte fossés, dolines, souterrains',
         'cmap': 'PuBu',
         'gradient': ('#fff7fb', '#ece7f2', '#d0d1e6', '#a5bddb', '#73a9cf',
@@ -176,6 +177,14 @@ VIZ_LEGENDS = {
                      '#fc4d2a', '#e2191c', '#bb0026', '#800026'),
         'ticks': ('0', '1'),
     },
+    'relief_oriente': {
+        'title': 'Relief orienté (openness locale × orientation)',
+        'legend': 'Clarté = micro-relief (openness locale 5–20 m + ombrage)\nClair = bosse, crête | Sombre = creux, fossé\nTeinte = orientation de la pente\nÉchelle fixe — couleurs homogènes entre tuiles',
+        'description': 'Openness sur MNT détendancé (σ 10 m), rayons 5/10/20 m, 16 directions ; teinte CIELAB = aspect',
+        'cmap': None,
+        'gradient': None,
+        'ticks': None,
+    },
     'ortho': {
         'title': 'Photographie Aérienne IGN',
         'legend': 'Orthophotographie\nImage aérienne',
@@ -197,11 +206,12 @@ VIZ_LEGENDS = {
 # Couches activées par défaut à l'ouverture de la carte : seules celles-ci
 # sont allumées, toutes les autres sont désactivées (mais restent accessibles
 # via le panneau). Ordre = pile basse → haute : une couche opaque (opacité 1)
-# doit rester en bas pour ne pas masquer les autres.
-DEFAULT_LAYERS = ('aspect', 'positive_openness')
+# doit rester en bas pour ne pas masquer les autres. Le relief orienté
+# fusionne déjà openness et aspect en une seule image (plus d'empilement).
+DEFAULT_LAYERS = ('relief_oriente',)
 
 # Opacité par défaut des couches allumées (0–1) ; les absentes restent à 1.
-DEFAULT_OPACITY = {'aspect': 0.74}
+DEFAULT_OPACITY = {}
 
 # Mode de fusion (mix-blend-mode CSS) par couche : les couches d'ombrage en
 # niveaux de gris fondent en « multiply » — elles assombrissent le relief sans
@@ -218,7 +228,7 @@ DEFAULT_VIZ = 'positive_openness'
 # référencées). None = toutes les visualisations présentes sur disque.
 # Le sélecteur de génération/régénération est piloté par le même registre
 # La pente y reste proposée mais éteinte par défaut.
-PANEL_VIZ = ('slope', 'aspect', 'positive_openness')
+PANEL_VIZ = ('relief_oriente', 'slope', 'aspect', 'positive_openness')
 
 # Correspondance mot-clé de fichier de sortie → nom d'étape --only du pipeline
 # (les trois visualisations dont le nom de sortie diffère du nom d'étape,
@@ -306,7 +316,7 @@ _MID_THUMB_SIZE = 640
 _VIZ_FALLBACK_ORDER = [
     'hillshade_multi', 'svf', 'slope', 'mslrm', 'positive_openness',
     'negative_openness', 'aspect', 'sailore', 'roughness',
-    'wavelet', 'flow_acc', 'solar', 'anomaly', 'ortho', 'topo',
+    'wavelet', 'flow_acc', 'solar', 'anomaly', 'relief_oriente', 'ortho', 'topo',
 ]
 
 # Regex pour parser les coordonnées tuile dans le basename LHD.
diff --git a/lidar_pipeline/pipeline.py b/lidar_pipeline/pipeline.py
index cb44a55..263fb17 100644
--- a/lidar_pipeline/pipeline.py
+++ b/lidar_pipeline/pipeline.py
@@ -87,6 +87,7 @@ from .visualizations import (
     generate_svf,
     generate_flow_accumulation,
     generate_anomaly_mask,
+    generate_relief_oriente,
 )
 from .gpu import gpu_cleanup, num_gpus, available_gpu_ids, restrict_gpus, safe_gpu_call
 from .ign import generate_ign_overlay
@@ -110,6 +111,7 @@ VIZ_STEPS = [
     ('flow_acc',   generate_flow_accumulation),
     ('solar',      generate_solar),
     ('anomaly',    generate_anomaly_mask),
+    ('relief_oriente', generate_relief_oriente),
     ('ortho',      lambda d, b, v, r: generate_ign_overlay(
                     d, b, v, r,
                     layer='ORTHOIMAGERY.ORTHOPHOTOS',
diff --git a/lidar_pipeline/rendering.py b/lidar_pipeline/rendering.py
index 85b1b5a..daaaa81 100644
--- a/lidar_pipeline/rendering.py
+++ b/lidar_pipeline/rendering.py
@@ -166,7 +166,9 @@ COLORMAPS = {
 RGB_LEGENDS = {
     'ortho': {},
     'topo': {},
+    'relief_oriente': {},   # RGB calculé (visualizations.generate_relief_oriente)
 }
+RGB_KEYWORDS = tuple(RGB_LEGENDS)
 
 # Fusion des textes de légende (titre / lecture du rendu / méthode de calcul)
 # depuis la source unique VIZ_LEGENDS (index.py, sans dépendance lourde).
@@ -406,7 +408,7 @@ def tif_to_png(tif_file, vis_dir, resolution, keep_tif=False, source_info=None,
 
     try:
         with rasterio.open(tif_file) as src:
-            is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo'))
+            is_rgb = src.count >= 3 and any(k in str(tif_file) for k in RGB_KEYWORDS)
 
             if is_rgb:
                 data = src.read([1, 2, 3])
@@ -840,7 +842,7 @@ def tif_to_crop(tif_file, vis_dir, resolution, keep_tif=False, quality=60, outpu
 
     try:
         with rasterio.open(tif_file) as src:
-            is_rgb = src.count >= 3 and any(k in str(tif_file) for k in ('ortho', 'topo'))
+            is_rgb = src.count >= 3 and any(k in str(tif_file) for k in RGB_KEYWORDS)
 
             if is_rgb:
                 data = src.read([1, 2, 3])
diff --git a/lidar_pipeline/tests/test_mapserve.py b/lidar_pipeline/tests/test_mapserve.py
index 61e86de..59fe146 100644
--- a/lidar_pipeline/tests/test_mapserve.py
+++ b/lidar_pipeline/tests/test_mapserve.py
@@ -198,15 +198,15 @@ def test_wmts_capabilities(tmp_path, monkeypatch):
 
 def test_map_meta(tmp_path, monkeypatch):
     """/api/map/meta décrit les couches et le contrat de tuilage de l'interface."""
-    mapserve = _setup(tmp_path, monkeypatch)
+    mapserve = _setup(tmp_path, monkeypatch, layers=("relief_oriente", "aspect", "slope"))
     meta = mapserve.map_meta()
-    assert {l["key"] for l in meta["layers"]} == {"aspect", "slope"}
+    assert {l["key"] for l in meta["layers"]} == {"relief_oriente", "aspect", "slope"}
     assert all(l["label"] for l in meta["layers"])
     assert meta["tile_url"] == "tiles/{layer}/{z}/{x}/{y}@2x.webp"
     assert meta["tile_size"] == 512 and meta["zoom_offset"] == -1
     assert meta["max_native_zoom"] == 18
     assert meta["bounds"] and meta["stamp"] > 0
-    assert "aspect" in meta["default_layers"]
+    assert meta["default_layers"] == ["relief_oriente"]
 
 
 def test_map_tile_info(tmp_path, monkeypatch):
diff --git a/lidar_pipeline/tests/test_pipeline.py b/lidar_pipeline/tests/test_pipeline.py
index c1e02c8..7af6cf7 100644
--- a/lidar_pipeline/tests/test_pipeline.py
+++ b/lidar_pipeline/tests/test_pipeline.py
@@ -20,9 +20,9 @@ class TestVizSteps:
         assert len(names) == len(set(names)), "VIZ_STEPS has duplicate names"
 
     def test_expected_visualization_count(self):
-        """Should have 15 visualizations (13 terrain + ortho + topo)."""
+        """Should have 16 visualizations (14 terrain + ortho + topo)."""
         from lidar_pipeline.pipeline import VIZ_STEPS
-        assert len(VIZ_STEPS) == 15
+        assert len(VIZ_STEPS) == 16
 
     def test_ortho_and_topo_present(self):
         from lidar_pipeline.pipeline import VIZ_STEPS
diff --git a/lidar_pipeline/tests/test_rendering.py b/lidar_pipeline/tests/test_rendering.py
index be24fd9..0993fcc 100644
--- a/lidar_pipeline/tests/test_rendering.py
+++ b/lidar_pipeline/tests/test_rendering.py
@@ -39,8 +39,9 @@ class TestColormaps:
             'neg_open': 'negative_openness',
             'hillshade': 'hillshade_multi',
         }
-        # IGN overlays (ortho, topo) are RGB images — no colormap needed
-        skip = {'ortho', 'topo'}
+        # Images RGB (fonds IGN, relief orienté) — pas de colormap
+        from lidar_pipeline.rendering import RGB_KEYWORDS
+        skip = set(RGB_KEYWORDS)
         for name, _ in VIZ_STEPS:
             if name in skip:
                 continue
diff --git a/lidar_pipeline/tests/test_visualizations.py b/lidar_pipeline/tests/test_visualizations.py
index f7fecf8..0cc4be4 100644
--- a/lidar_pipeline/tests/test_visualizations.py
+++ b/lidar_pipeline/tests/test_visualizations.py
@@ -126,6 +126,146 @@ class TestOpenness:
         corr = np.corrcoef(a[m], b[m])[0, 1]
         assert corr > 0.97, f"corrélation openness décimée/native trop faible : {corr:.3f}"
 
+    @pytest.mark.parametrize("positive", [True, False])
+    def test_scale_independent_of_rest_of_tile(self, tmp_path, tmp_output_dir, positive):
+        """Même relief local = même valeur, quel que soit le reste de la dalle.
+
+        Deux MNT identiques sur leur moitié ouest ; l'un porte en plus une
+        colline à l'est, à plus de 100 m (rayon max) de la zone comparée. Une
+        normalisation par dalle (z-score) décalait toute l'échelle et rendait
+        les mosaïques non jointives ; avec les références figées, la moitié
+        ouest doit sortir identique.
+        """
+        import rasterio
+        from rasterio.transform import from_bounds
+        from lidar_pipeline.visualizations import generate_openness
+
+        size = 300
+        x = np.arange(size, dtype=float)
+        X, Y = np.meshgrid(x, x)
+        flat = 100.0 + 2.0 * np.exp(-((X - 60)**2 + (Y - 150)**2) / (2 * 15**2))
+        hill = flat + 40.0 * np.exp(-((X - 260)**2 + (Y - 150)**2) / (2 * 20**2))
+        results = []
+        for name, dem in (("flat", flat), ("hill", hill)):
+            f = tmp_path / f"{name}.tif"
+            with rasterio.open(f, 'w', driver='GTiff', height=size, width=size,
+                               count=1, dtype='float32', crs='EPSG:2154',
+                               transform=from_bounds(660000, 6700000, 660300, 6700300,
+                                                     size, size)) as dst:
+                dst.write(dem.astype('float32'), 1)
+            out = generate_openness(f, name, tmp_output_dir, 1.0, positive=positive)
+            with rasterio.open(out) as src:
+                results.append(src.read(1))
+        # Tolérance : résidu d'interpolation de la grille décimée (≈0,01) ;
+        # un z-score par dalle décalerait toute la zone de plusieurs dixièmes.
+        west = np.s_[:, :100]
+        np.testing.assert_allclose(results[0][west], results[1][west], atol=0.05)
+
+
+def _write_dem(path, dem, res=1.0):
+    import rasterio
+    from rasterio.transform import from_origin
+    with rasterio.open(path, 'w', driver='GTiff', height=dem.shape[0], width=dem.shape[1],
+                       count=1, dtype='float32', crs='EPSG:2154',
+                       transform=from_origin(660000, 6700300, res, res)) as dst:
+        dst.write(dem.astype('float32'), 1)
+    return path
+
+
+class TestReliefOriente:
+    def test_generates_rgb_uint8(self, synthetic_dem, tmp_output_dir):
+        import rasterio
+        from lidar_pipeline.visualizations import generate_relief_oriente
+        result = generate_relief_oriente(synthetic_dem, "test", tmp_output_dir, 5.0)
+        assert result is not None and result.name == "test_relief_oriente.tif"
+        with rasterio.open(result) as src, rasterio.open(synthetic_dem) as dem:
+            assert src.count == 3 and src.dtypes[0] == 'uint8'
+            assert (src.height, src.width) == (dem.height, dem.width)
+            rgb = src.read()
+        assert rgb.std() > 5, "image uniforme : ni relief ni orientation rendus"
+
+    def test_nodata_gets_fixed_color(self, tmp_path, tmp_output_dir):
+        import rasterio
+        from lidar_pipeline.visualizations import generate_relief_oriente, RELIEF_NODATA_RGB
+        x = np.arange(120, dtype=float)
+        dem = 100 + 0.05 * x[None, :] + 0.02 * x[:, None]
+        dem[10:20, 10:20] = np.nan
+        out = generate_relief_oriente(_write_dem(tmp_path / "d.tif", dem), "t", tmp_output_dir, 1.0)
+        with rasterio.open(out) as src:
+            rgb = src.read()
+        assert tuple(rgb[:, 15, 15]) == RELIEF_NODATA_RGB
+
+    def test_horizon_kernels_agree(self):
+        """numba et numpy (et le noyau CUDA, même code) donnent le même angle."""
+        pytest.importorskip("numba")
+        from lidar_pipeline.visualizations import (
+            _horizon_rays, _mean_horizon_numba, _mean_horizon_numpy)
+        rng = np.random.default_rng(0)
+        dem = rng.normal(0, 0.3, (60, 70)).astype(np.float32)
+        dem[20:30, 30:40] += 2.0
+        offs, dist, cps = _horizon_rays(0.8, 16, (5, 10, 20))
+        a = _mean_horizon_numba(dem, offs, dist, cps)
+        b = _mean_horizon_numpy(dem, offs, dist, cps)
+        np.testing.assert_allclose(a, b, atol=1e-5)
+
+    def test_matches_legacy_ray_trace_geometry(self):
+        """Même géométrie de rayons que _ray_trace_horizons (8 directions)."""
+        from lidar_pipeline.visualizations import (
+            _horizon_rays, _mean_horizon_numpy, _ray_trace_horizons)
+        rng = np.random.default_rng(1)
+        dem = rng.normal(0, 0.5, (50, 50)).astype(np.float32)
+        radii = (5, 10, 20)
+        pos, _ = _ray_trace_horizons(dem, 50, 50, 1.0, 8, 20, list(radii))
+        legacy = np.mean(pos, axis=(0, 1))
+        offs, dist, cps = _horizon_rays(1.0, 8, radii)
+        np.testing.assert_allclose(_mean_horizon_numpy(dem, offs, dist, cps), legacy, atol=1e-5)
+
+    def test_seamless_independent_of_rest_of_tile(self, tmp_path, tmp_output_dir):
+        """Même relief local = mêmes couleurs, quel que soit le reste de la dalle
+        (support : rayon 20 m + détendance 4σ = 40 m, loin sous la bande de 100 m)."""
+        import rasterio
+        from lidar_pipeline.visualizations import generate_relief_oriente
+        size = 300
+        X, Y = np.meshgrid(np.arange(size, dtype=float), np.arange(size, dtype=float))
+        flat = 100.0 + 1.5 * np.exp(-((X - 60)**2 + (Y - 150)**2) / (2 * 6**2))
+        hill = flat + 40.0 * np.exp(-((X - 260)**2 + (Y - 150)**2) / (2 * 10**2))
+        imgs = []
+        for name, dem in (("flat", flat), ("hill", hill)):
+            out = generate_relief_oriente(_write_dem(tmp_path / f"{name}.tif", dem),
+                                          name, tmp_output_dir, 1.0)
+            with rasterio.open(out) as src:
+                imgs.append(src.read().astype(int))
+        diff = np.abs(imgs[0][:, :, :150] - imgs[1][:, :, :150])
+        assert diff.max() <= 1, f"écart de couleur loin de la colline : {diff.max()}"
+
+    def test_numba_colorize_matches_vectorized(self):
+        """Noyau CPU fusionné = chemin vectorisé (CuPy/numpy) à l'arrondi près."""
+        pytest.importorskip("numba")
+        from lidar_pipeline.gpu import xp_zoom
+        from lidar_pipeline.visualizations import (
+            _relief_colorize_numba, _relief_colorize_xp, _relief_lut, _pad_to)
+        rng = np.random.default_rng(2)
+        open_c = rng.uniform(0.3, 15, (30, 25)).astype(np.float32)
+        dx = rng.normal(0, 0.3, (120, 100)).astype(np.float32)
+        dy = rng.normal(0, 0.3, (120, 100)).astype(np.float32)
+        lut = _relief_lut()
+        a = _relief_colorize_numba(open_c, 4, dx, dy, lut).astype(int)
+        b = _relief_colorize_xp(np, _pad_to(np, xp_zoom(open_c, 4), 120, 100), dx, dy, lut).astype(int)
+        diff = np.abs(a - b).max(axis=2)
+        assert (diff > 3).mean() < 0.01, f"{(diff > 3).mean():.3%} pixels divergent"
+
+    def test_lut_lightness_is_monotonic_and_hue_neutral(self):
+        """Clarté croissante avec L* ; à L* fixé, toutes les teintes ont la même
+        luminance perçue (pas de faux relief dû à la couleur)."""
+        from lidar_pipeline.visualizations import _relief_lut
+        lut = _relief_lut().astype(float) / 255
+        lin = np.where(lut <= 0.04045, lut / 12.92, ((lut + 0.055) / 1.055) ** 2.4)
+        Y = lin @ np.array([0.2126, 0.7152, 0.0722])          # (256, 360)
+        assert np.all(np.diff(Y.mean(axis=1)) >= -1e-6)
+        Lstar = 116 * np.cbrt(Y[100:180]) - 16                # L* ≈ 40–70
+        # Écart résiduel : écrêtage de gamme sRGB de quelques teintes à C* 60
+        assert (Lstar.max(axis=1) - Lstar.min(axis=1)).max() < 6
+
 
 class TestMSLRM:
     def test_generates_tif(self, synthetic_dem, tmp_output_dir):
diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py
index 276cf5b..218c807 100644
--- a/lidar_pipeline/visualizations.py
+++ b/lidar_pipeline/visualizations.py
@@ -638,6 +638,19 @@ def generate_svf(dem_file, basename, vis_dir, resolution, shared=None):
 # ; facteur 2 ≈ ×8 plus rapide, rendu quasi identique. 1 = pleine résolution.
 OPENNESS_DOWNSAMPLE = 2
 
+# Références de normalisation de l'openness (degrés : moyenne, écart-type).
+# Le z-score par tuile rendait l'échelle non jointive — même ouverture
+# physique, couleur différente d'une dalle à l'autre selon le relief
+# environnant. Médianes inter-tuiles mesurées à 0,2 m (décimation ×2) sur
+# des dalles réparties sur le territoire (plaine, bocage, forêt, montagne,
+# volcans, delta, urbain, littoral). Références FIGÉES : même ouverture =
+# même couleur sur toutes les tuiles.
+# Rayons du lancé de rayons (mètres), moyennés à poids égaux.
+OPENNESS_RADII_M = (25, 50, 100)
+
+OPENNESS_POS_REF = (5.571, 4.112)
+OPENNESS_NEG_REF = (5.974, 4.070)
+
 
 def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, shared=None):
     """Positive/Negative Openness - multi-radius ray-tracing with std normalization.
@@ -645,7 +658,8 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
     Traces rays in 8 directions at 3 radii (25, 50, 100m) on a block-decimated
     grid (cf. OPENNESS_DOWNSAMPLE), then bilinearly resamples the result back
     to the requested resolution. Results are combined with equal weight across
-    radii, then normalized by standard deviation for cross-tile comparability.
+    radii, then normalized by FIXED references (OPENNESS_POS_REF /
+    OPENNESS_NEG_REF, degrees) so that adjacent tiles share one colour scale.
     """
     name = "positive_openness" if positive else "negative_openness"
     gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
@@ -657,7 +671,7 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
         dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
             _prepare_dem_for_raycast(dem_file, shared, resolution)
 
-        radii_m = [25, 50, 100]
+        radii_m = list(OPENNESS_RADII_M)
         n_dirs = 8
         full_rows, full_cols = rows, cols
 
@@ -699,12 +713,10 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
             openness_result = zoomed.astype(np.float32)
         openness_result[nan_mask] = np.nan
 
-        # Z-score (écarts locaux en sigmas) : unités comparables entre tuiles,
-        # plage de rendu fixe → mosaïque de couleur homogène
-        valid = openness_result[~nan_mask]
-        if len(valid) > 0:
-            std_val = max(np.nanstd(valid), 0.01)
-            openness_result = (openness_result - np.nanmean(valid)) / std_val
+        # Écart aux références figées, en sigmas : même ouverture physique =
+        # même valeur sur toutes les tuiles → mosaïque de couleur jointive
+        ref_mean, ref_std = OPENNESS_POS_REF if positive else OPENNESS_NEG_REF
+        openness_result = (openness_result - ref_mean) / ref_std
 
         _save_tif(output, openness_result, transform, crs)
         logger.info(f"  ✓ {name} terminé ({time.time()-t0:.1f}s){' [GPU]' if _gpu_mod.HAS_GPU else ''}")
@@ -714,6 +726,384 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
         return None
 
 
+# ============================================================
+# Relief orienté : openness locale × aspect en une seule image RGB
+# ============================================================
+#
+# Clarté (CIELAB L*) = micro-relief : openness positive calculée sur le MNT
+# détendancé (MNT − gaussienne 10 m), rayons courts 5/10/20 m, plus un léger
+# ombrage directionnel. Teinte = orientation de la pente (aspect) sur le
+# cercle CIELAB, à clarté constante : aucune couleur ne crée de faux relief.
+# Pas de statistique par dalle (échelle log FIXE) et un support total
+# (rayon 20 m + lissage 4σ = 40 m) inférieur à la bande de raccord de 100 m :
+# les dalles adjacentes se raccordent sans couture.
+#
+# Coût maîtrisé : détendance et lancé de rayons tournent sur une grille
+# décimée à ~0,8 m (bloc max, cf. OPENNESS_DOWNSAMPLE) — seule l'openness
+# finale est agrandie ; noyau dédié qui n'accumule que la moyenne des angles
+# (CuPy RawKernel sur GPU, numba parallèle sur CPU, numpy en dernier
+# recours) ; couleurs par table pré-calculée (L* × teinte) au lieu d'une
+# conversion Lab → sRGB par pixel.
+RELIEF_DETREND_M = 10.0          # σ de la gaussienne retirée au MNT
+RELIEF_RADII_M = (5, 10, 20)     # rayons du lancé de rayons, poids égaux
+RELIEF_N_DIRS = 16               # 16 directions : pas de losanges à 8 branches
+RELIEF_GRID_M = 0.8              # pas de la grille de calcul décimée
+RELIEF_OPEN_RANGE = (0.5, 12.0)  # degrés, échelle log fixe (jointive)
+RELIEF_CHROMA = 60.0             # chroma CIELAB max (atteint vers L* = 50)
+RELIEF_SHADE_WEIGHT = 0.35       # part de l'ombrage directionnel dans L*
+# Azimut de l'ombrage, exprimé dans le repère de l'aspect ci-dessous
+# (arctan2(dy, dx), dy vers le sud) — réglage validé sur la galerie de rendus.
+RELIEF_SHADE_AZIMUTH = 315.0
+RELIEF_SHADE_ALTITUDE = 45.0
+RELIEF_NODATA_RGB = (38, 38, 41)
+
+_RELIEF_LUT = None
+
+
+def _lab_to_srgb(L, a, b):
+    """CIELAB (D65) → sRGB [0, 1], écrêté."""
+    fy = (L + 16) / 116
+    fx = fy + a / 500
+    fz = fy - b / 200
+    finv = lambda t: np.where(t > 6 / 29, t ** 3, 3 * (6 / 29) ** 2 * (t - 4 / 29))
+    X, Y, Z = 0.95047 * finv(fx), finv(fy), 1.08883 * finv(fz)
+    rgb = np.stack([3.2406 * X - 1.5372 * Y - 0.4986 * Z,
+                    -0.9689 * X + 1.8758 * Y + 0.0415 * Z,
+                    0.0557 * X - 0.2040 * Y + 1.0570 * Z], -1)
+    rgb = np.where(rgb <= 0.0031308, 12.92 * rgb,
+                   1.055 * np.clip(rgb, 0, None) ** (1 / 2.4) - 0.055)
+    return np.clip(rgb, 0, 1)
+
+
+def _relief_lut():
+    """Table (256 niveaux de L*, 360 teintes) → sRGB uint8, calculée une fois.
+
+    La chroma décroît vers le noir et le blanc (L*(100−L*)/2500) : les
+    couleurs restent dans la gamme sRGB au lieu d'être écrêtées.
+    """
+    global _RELIEF_LUT
+    if _RELIEF_LUT is None:
+        L = np.linspace(0, 100, 256)[:, None]
+        h = np.radians(np.arange(360))[None, :]
+        C = RELIEF_CHROMA * np.clip(L * (100 - L) / 2500.0, 0, 1)
+        _RELIEF_LUT = (_lab_to_srgb(np.broadcast_to(L, (256, 360)), C * np.cos(h), C * np.sin(h))
+                       * 255 + 0.5).astype(np.uint8)
+    return _RELIEF_LUT
+
+
+def _horizon_rays(res, n_dirs, radii_m):
+    """Décalages (col, ligne) de chaque pas de rayon, distances et checkpoints.
+
+    Même géométrie que _ray_trace_horizons_core : direction k à l'angle
+    2πk/n, pas arrondis au pixel, distance = pas × résolution.
+    """
+    max_step = max(1, int(max(radii_m) / res))
+    angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
+    steps = np.arange(1, max_step + 1)
+    offs = np.empty((n_dirs, max_step, 2), dtype=np.int32)
+    offs[:, :, 0] = np.rint(np.cos(angles)[:, None] * steps)
+    offs[:, :, 1] = np.rint(np.sin(angles)[:, None] * steps)
+    dist = (steps * res).astype(np.float32)
+    cps = np.array(sorted(max(1, min(int(r / res), max_step)) for r in radii_m), dtype=np.int32)
+    return offs, dist, cps
+
+
+_HORIZON_CUDA_SRC = r'''
+extern "C" __global__
+void mean_horizon(const float* dem, const int* offs, const float* dist,
+                  const int* cps, const int rows, const int cols,
+                  const int nd, const int ns, const int nr, float* out) {
+    long idx = (long)blockDim.x * blockIdx.x + threadIdx.x;
+    if (idx >= (long)rows * cols) return;
+    int i = idx / cols, j = idx % cols;
+    float z0 = dem[idx], acc = 0.f;
+    for (int d = 0; d < nd; ++d) {
+        float run = 0.f;
+        int k = 0;
+        for (int s = 0; s < ns; ++s) {
+            int ii = i + offs[(d * ns + s) * 2 + 1];
+            int jj = j + offs[(d * ns + s) * 2];
+            if (ii >= 0 && ii < rows && jj >= 0 && jj < cols)
+                run = fmaxf(run, (dem[(long)ii * cols + jj] - z0) / dist[s]);
+            while (k < nr && s + 1 >= cps[k]) { acc += atanf(run); ++k; }
+        }
+        while (k < nr) { acc += atanf(run); ++k; }
+    }
+    out[idx] = acc / (float)(nd * nr);
+}
+'''
+_horizon_cuda_kernel = None
+_horizon_numba_kernel = None
+
+
+def _mean_horizon_gpu(dem, offs, dist, cps):
+    """Un thread CUDA par pixel ; tout reste en registres (aucun tableau
+    intermédiaire par direction ou par rayon)."""
+    global _horizon_cuda_kernel
+    cp = _gpu_mod._cp
+    if _horizon_cuda_kernel is None:
+        _horizon_cuda_kernel = cp.RawKernel(_HORIZON_CUDA_SRC, 'mean_horizon')
+    rows, cols = dem.shape
+    d_dem = cp.ascontiguousarray(cp.asarray(dem, dtype=cp.float32))
+    out = cp.empty((rows, cols), dtype=cp.float32)
+    n = rows * cols
+    threads = 256
+    _horizon_cuda_kernel(((n + threads - 1) // threads,), (threads,),
+                         (d_dem, cp.asarray(offs), cp.asarray(dist), cp.asarray(cps),
+                          np.int32(rows), np.int32(cols), np.int32(offs.shape[0]),
+                          np.int32(offs.shape[1]), np.int32(len(cps)), out))
+    return out
+
+
+def _mean_horizon_numba(dem, offs, dist, cps):
+    """Même noyau en numba parallèle (une ligne de pixels par thread CPU)."""
+    global _horizon_numba_kernel
+    if _horizon_numba_kernel is None:
+        from numba import njit, prange
+
+        @njit(parallel=True, cache=True, fastmath=True)
+        def _kernel(dem, offs, dist, cps):
+            rows, cols = dem.shape
+            nd, ns, nr = offs.shape[0], offs.shape[1], cps.shape[0]
+            out = np.empty((rows, cols), dtype=np.float32)
+            for i in prange(rows):
+                for j in range(cols):
+                    z0 = dem[i, j]
+                    acc = 0.0
+                    for d in range(nd):
+                        run = 0.0
+                        k = 0
+                        for s in range(ns):
+                            ii = i + offs[d, s, 1]
+                            jj = j + offs[d, s, 0]
+                            if ii >= 0 and ii < rows and jj >= 0 and jj < cols:
+                                t = (dem[ii, jj] - z0) / dist[s]
+                                if t > run:
+                                    run = t
+                            while k < nr and s + 1 >= cps[k]:
+                                acc += math.atan(run)
+                                k += 1
+                        while k < nr:
+                            acc += math.atan(run)
+                            k += 1
+                    out[i, j] = acc / (nd * nr)
+            return out
+        _horizon_numba_kernel = _kernel
+    return _horizon_numba_kernel(np.ascontiguousarray(dem, dtype=np.float32), offs, dist, cps)
+
+
+def _mean_horizon_numpy(dem, offs, dist, cps):
+    """Repli vectorisé (sans numba ni GPU) : un décalage de tableau par pas."""
+    rows, cols = dem.shape
+    pad = int(np.abs(offs).max())
+    padded = np.pad(dem.astype(np.float32), pad, constant_values=np.nan)
+    acc = np.zeros((rows, cols), dtype=np.float64)
+    for d in range(offs.shape[0]):
+        run = np.zeros((rows, cols), dtype=np.float32)
+        k = 0
+        for s in range(offs.shape[1]):
+            px, py = offs[d, s]
+            view = padded[pad + py:pad + py + rows, pad + px:pad + px + cols]
+            run = np.fmax(run, (view - dem) / dist[s])
+            while k < len(cps) and s + 1 >= cps[k]:
+                acc += np.arctan(run)
+                k += 1
+        while k < len(cps):
+            acc += np.arctan(run)
+            k += 1
+    return (acc / (offs.shape[0] * len(cps))).astype(np.float32)
+
+
+def _mean_horizon_angle(dem, res, n_dirs, radii_m):
+    """Moyenne (directions × rayons) de l'angle d'horizon positif, en radians.
+
+    GPU (CuPy RawKernel) si disponible — le résultat reste alors sur le GPU —,
+    sinon numba, sinon numpy. Un échec GPU (compilation NVRTC, VRAM) bascule
+    sur le CPU pour ce calcul sans perdre la visualisation.
+    """
+    offs, dist, cps = _horizon_rays(res, n_dirs, radii_m)
+    if _gpu_mod.is_gpu_active():
+        try:
+            return _mean_horizon_gpu(dem, offs, dist, cps), "GPU"
+        except Exception as e:
+            logger.warning(f"    ⚠ Noyau GPU relief orienté indisponible ({e}) — repli CPU")
+    dem = to_cpu(dem)
+    try:
+        return _mean_horizon_numba(dem, offs, dist, cps), "numba"
+    except ImportError:
+        return _mean_horizon_numpy(dem, offs, dist, cps), "numpy"
+
+
+def _pad_to(m, arr, rows, cols):
+    """Complète par répétition du bord (dimensions non multiples du facteur)."""
+    if arr.shape == (rows, cols):
+        return arr
+    return m.pad(arr[:rows, :cols], ((0, max(0, rows - arr.shape[0])), (0, max(0, cols - arr.shape[1]))),
+                 mode='edge')
+
+
+def _relief_params():
+    """Constantes scalaires du rendu (ombrage sans trigonométrie par pixel).
+
+    cos(pente) = 1/√(1+g²) et sin(pente)·cos(az − aspect) = (cos az·dx +
+    sin az·dy)/√(1+g²), avec aspect = arctan2(dy, dx) : l'ombrage ne coûte
+    qu'une racine par pixel.
+    """
+    zen = math.radians(90.0 - RELIEF_SHADE_ALTITUDE)
+    az = math.radians(RELIEF_SHADE_AZIMUTH)
+    lo, hi = RELIEF_OPEN_RANGE
+    return (math.cos(zen), math.sin(zen) * math.cos(az), math.sin(zen) * math.sin(az),
+            lo, 1.0 / math.log(hi / lo), RELIEF_SHADE_WEIGHT)
+
+
+def _relief_colorize_xp(m, openness, dx, dy, lut):
+    """Colorisation vectorisée (CuPy sur GPU, numpy en repli)."""
+    cz, sa_x, sa_y, lo, inv_log, w = _relief_params()
+    inv = 1.0 / m.sqrt(1.0 + dx * dx + dy * dy)
+    shade = m.clip((cz + sa_x * dx + sa_y * dy) * inv / cz, 0, 1.25) / 1.25
+    t = m.clip(m.log(m.maximum(openness, 1e-3) / lo) * inv_log, 0, 1)
+    L = 12.0 + 84.0 * ((1 - w) * (1 - t) + w * shade)
+    del inv, shade, t
+    li = m.clip(m.rint(L * 2.55), 0, 255).astype(m.int32)
+    hue = m.mod(m.rint(m.degrees(m.arctan2(dy, dx))), 360).astype(m.int32)
+    return lut[li, hue]
+
+
+_relief_color_numba_kernel = None
+
+
+def _relief_colorize_numba(open_c, f, dx, dy, lut):
+    """Colorisation CPU en une passe (numba parallèle) : l'openness est lue
+    sur la grille décimée par interpolation bilinéaire alignée sur les centres
+    de pixels (équivalent de xp_zoom), sans tableau pleine résolution
+    intermédiaire."""
+    global _relief_color_numba_kernel
+    if _relief_color_numba_kernel is None:
+        from numba import njit, prange
+
+        @njit(parallel=True, cache=True, fastmath=True)
+        def _kernel(open_c, f, dx, dy, lut, cz, sa_x, sa_y, lo, inv_log, w):
+            rows, cols = dx.shape
+            rc, cc = open_c.shape
+            out = np.empty((rows, cols, 3), dtype=np.uint8)
+            for i in prange(rows):
+                u = (i + 0.5) / f - 0.5
+                u = min(max(u, 0.0), rc - 1.0)
+                i0 = int(u)
+                i1 = min(i0 + 1, rc - 1)
+                fu = u - i0
+                for j in range(cols):
+                    v = (j + 0.5) / f - 0.5
+                    v = min(max(v, 0.0), cc - 1.0)
+                    j0 = int(v)
+                    j1 = min(j0 + 1, cc - 1)
+                    fv = v - j0
+                    o = ((1 - fu) * ((1 - fv) * open_c[i0, j0] + fv * open_c[i0, j1])
+                         + fu * ((1 - fv) * open_c[i1, j0] + fv * open_c[i1, j1]))
+                    gx = dx[i, j]
+                    gy = dy[i, j]
+                    inv = 1.0 / math.sqrt(1.0 + gx * gx + gy * gy)
+                    sh = (cz + sa_x * gx + sa_y * gy) * inv / cz
+                    sh = min(max(sh, 0.0), 1.25) / 1.25
+                    t = math.log(max(o, 1e-3) / lo) * inv_log
+                    t = min(max(t, 0.0), 1.0)
+                    L = 12.0 + 84.0 * ((1 - w) * (1 - t) + w * sh)
+                    li = min(max(int(round(L * 2.55)), 0), 255)
+                    h = int(round(math.degrees(math.atan2(gy, gx)))) % 360
+                    out[i, j, 0] = lut[li, h, 0]
+                    out[i, j, 1] = lut[li, h, 1]
+                    out[i, j, 2] = lut[li, h, 2]
+            return out
+        _relief_color_numba_kernel = _kernel
+    return _relief_color_numba_kernel(np.ascontiguousarray(open_c), float(f),
+                                      np.ascontiguousarray(dx, dtype=np.float32),
+                                      np.ascontiguousarray(dy, dtype=np.float32),
+                                      lut, *_relief_params())
+
+
+def generate_relief_oriente(dem_file, basename, vis_dir, resolution, shared=None):
+    """Relief orienté : image RGB unique fusionnant openness locale et aspect.
+
+    L* (clarté) = openness positive locale (MNT détendancé, rayons courts)
+    65 % + ombrage directionnel 35 % ; teinte = orientation de la pente ;
+    chroma CIELAB fixe. Sortie : GeoTIFF RGB uint8 (rendu tel quel, comme
+    les fonds IGN).
+    """
+    gpu_tag = " [GPU]" if _gpu_mod.HAS_GPU else ""
+    logger.info(f"  → Relief orienté (openness locale × aspect){gpu_tag}...")
+    t0 = time.time()
+    output = vis_dir / f"{basename}_relief_oriente.tif"
+
+    try:
+        dem, dem_np, rows, cols, res, nan_mask, transform, crs = \
+            _prepare_dem_for_raycast(dem_file, shared, resolution)
+        res = float(res)
+        use_gpu = _gpu_mod.is_gpu_active()
+        filled = (shared.filled_gpu if shared is not None and use_gpu else None)
+        if filled is None:
+            filled = to_gpu(dem) if use_gpu else dem.astype(np.float32)
+
+        # 1. Grille décimée : bloc max (préserve les reliefs qui bornent
+        #    l'horizon) et bloc moyen (support de la tendance à retirer)
+        f = max(1, int(round(RELIEF_GRID_M / res)))
+        if rows < 2 * f or cols < 2 * f:
+            f = 1
+        r2, c2 = (rows // f) * f, (cols // f) * f
+        blocks = filled[:r2, :c2].reshape(r2 // f, f, c2 // f, f)
+        coarse_max = blocks.max(axis=(1, 3))
+        coarse_mean = blocks.mean(axis=(1, 3))
+        res_c = res * f
+
+        # 2. Détendance (MNT − gaussienne) sur la grille décimée
+        trend = xp_gaussian_filter(coarse_mean, RELIEF_DETREND_M / res_c)
+        detrended = coarse_max - trend
+        del blocks, coarse_max, coarse_mean, trend
+        t_prep = time.time() - t0
+
+        # 3. Openness locale : angle d'horizon moyen (radians) sur la grille décimée
+        t1 = time.time()
+        angle, engine = _mean_horizon_angle(detrended, res_c, RELIEF_N_DIRS, RELIEF_RADII_M)
+        t_ray = time.time() - t1
+        on_gpu = engine == "GPU"
+        del detrended
+
+        # 4. Couleurs : agrandissement de l'openness, ombrage, table L* × teinte
+        t2 = time.time()
+        if shared is not None:
+            dy, dx = shared.dy, shared.dx
+        else:
+            dy, dx = np.gradient(to_cpu(filled).astype(np.float32), res)
+        lut = _relief_lut()
+        if on_gpu:
+            openness = _gpu_mod.xp_zoom(xp.degrees(angle), f, order=1) if f > 1 else xp.degrees(angle)
+            openness = _pad_to(xp, xp.asarray(openness), rows, cols)
+            rgb = to_cpu(_relief_colorize_xp(xp, openness, to_gpu(dx), to_gpu(dy), xp.asarray(lut)))
+            del openness
+            gpu_cleanup()
+            engine_color = "GPU"
+        else:
+            angle = to_cpu(angle)
+            try:
+                rgb = _relief_colorize_numba(np.degrees(angle).astype(np.float32), f, dx, dy, lut)
+                engine_color = "numba"
+            except ImportError:
+                openness = _gpu_mod.xp_zoom(np.degrees(angle), f, order=1) if f > 1 else np.degrees(angle)
+                rgb = _relief_colorize_xp(np, _pad_to(np, openness, rows, cols), dx, dy, lut)
+                engine_color = "numpy"
+        del angle, dx, dy
+        t_color = time.time() - t2
+
+        rgb[nan_mask] = RELIEF_NODATA_RGB
+        _save_tif(output, np.moveaxis(rgb, -1, 0), transform, crs, dtype='uint8', count=3)
+        logger.info(f"  ✓ Relief orienté terminé ({time.time()-t0:.1f}s : "
+                    f"préparation {t_prep:.1f}s, rayons {t_ray:.1f}s [{engine}], "
+                    f"couleurs {t_color:.1f}s [{engine_color}])")
+        return output
+    except Exception as e:
+        logger.error(f"  ✗ Erreur relief orienté: {e}", exc_info=True)
+        return None
+
+
 def generate_mslrm(dem_file, basename, vis_dir, resolution, shared=None):
     """Multi-Scale Relief Model (MSRM) - LRM at adaptive scales combined (GPU if available).