diff --git a/lidar_pipeline/export_pdf.py b/lidar_pipeline/export_pdf.py index 5da5edc..21ddf79 100644 --- a/lidar_pipeline/export_pdf.py +++ b/lidar_pipeline/export_pdf.py @@ -192,3 +192,116 @@ def compose_l93(output_dir, bbox, px_size, layer=LAYER): if mask.getextrema() != (255, 255): canvas.paste(_hatch(canvas.size), (0, 0), ImageOps.invert(mask)) return canvas, mask, sorted(cells) + + +ROSE_L = 64.0 +ROSE_CHROMA = 60.0 # = visualizations.RELIEF_CHROMA +# Classes de densité de points sol (pts/m²) : rouge = donnée faible. +DENSITY_CLASSES = [ + (0.0, "moins de 1", "#d7301f"), + (1.0, "1 à 3", "#fc8d59"), + (3.0, "3 à 6", "#fdcc8a"), + (6.0, "6 à 10", "#a1d99b"), + (10.0, "10 et plus", "#31a354"), +] + + +def lab_to_rgb(L, a, b): + """CIELAB (D65) → sRGB 8 bits, même formule que la carte (labToRgb).""" + fy = (L + 16) / 116 + fx, fz = fy + a / 500, fy - b / 200 + + def finv(t): + return t ** 3 if t > 6 / 29 else 3 * (6 / 29) ** 2 * (t - 4 / 29) + + X, Y, Z = 0.95047 * finv(fx), finv(fy), 1.08883 * finv(fz) + lin = (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) + out = [] + for c in lin: + v = 12.92 * c if c <= 0.0031308 else 1.055 * max(c, 0.0) ** (1 / 2.4) - 0.055 + out.append(int(round(min(1.0, max(0.0, v)) * 255))) + return tuple(out) + + +def rose_color(compass_deg): + """Couleur d'une orientation de pente (0 = N, sens horaire), comme la rose de la carte.""" + chroma = ROSE_CHROMA * min(1.0, ROSE_L * (100 - ROSE_L) / 2500) + h = math.radians((compass_deg + 90) % 360) + return lab_to_rgb(ROSE_L, chroma * math.cos(h), chroma * math.sin(h)) + + +def density_color(v): + """Couleur de classe d'une densité sol (pts/m²).""" + color = DENSITY_CLASSES[0][2] + for low, _label, c in DENSITY_CLASSES: + if v >= low: + color = c + return color + + +def _pdf_text(s): + """Texte encodable par les polices standard (WinAnsi/cp1252).""" + return "".join(ch if ch.encode("cp1252", "ignore") else "?" for ch in str(s)) + + +def _cells_of_bbox(bbox): + """Dalles LHD 1 km intersectant une emprise L93.""" + c0, c1 = int(math.floor(bbox[0] / 1000)), int(math.ceil(bbox[2] / 1000)) + r0, r1 = int(math.floor(bbox[1] / 1000)) + 1, int(math.ceil(bbox[3] / 1000)) + return [(c, r) for c in range(c0, c1) for r in range(r1, r0 - 1, -1)] + + +def zone_quality(bbox, table, cells_with_relief): + """Agrège la qualité des dalles sur une emprise (pondérée par la surface).""" + from .index import parse_basename_coords + from .quality import DENSITY_CELL_M + + by_cell = {} + for base, data in table.items(): + coords = parse_basename_coords(base) + if coords is not None: + by_cell[tuple(coords)] = data + cells = _cells_of_bbox(bbox) + relief = set(map(tuple, cells_with_relief)) + total_w = dens_w = empty_w = 0.0 + dmin = None + starts, ends, sources = [], [], set() + grid_cells, missing_q = [], [] + for col, row in cells: + q = by_cell.get((col, row)) + if q is None: + missing_q.append((col, row)) + continue + x0, y1 = col * 1000.0, row * 1000.0 + grid = q.get("density_grid") or [] + for j, line in enumerate(grid): + for i, v in enumerate(line): + gx0, gx1 = x0 + i * DENSITY_CELL_M, x0 + (i + 1) * DENSITY_CELL_M + gy1, gy0 = y1 - j * DENSITY_CELL_M, y1 - (j + 1) * DENSITY_CELL_M + ox = min(gx1, bbox[2]) - max(gx0, bbox[0]) + oy = min(gy1, bbox[3]) - max(gy0, bbox[1]) + if ox <= 0 or oy <= 0: + continue + area = ox * oy + total_w += area + dens_w += v * area + empty_w += float(q.get("empty_fraction") or 0.0) * area + dmin = v if dmin is None else min(dmin, v) + grid_cells.append((gx0, gy0, gx1, gy1, v)) + if q.get("acq_start"): + starts.append(q["acq_start"]); ends.append(q["acq_end"] or q["acq_start"]) + sources.add(q.get("acq_source")) + return { + "density_mean": dens_w / total_w if total_w else None, + "density_min": dmin, + "empty_fraction": empty_w / total_w if total_w else None, + "acq_start": min(starts) if starts else None, + "acq_end": max(ends) if ends else None, + "acq_sources": sources, + "cells": cells, + "missing_relief": [c for c in cells if c not in relief], + "missing_quality": missing_q, + "grid_cells": grid_cells, + } diff --git a/lidar_pipeline/tests/test_export_pdf.py b/lidar_pipeline/tests/test_export_pdf.py index 0cb21f9..346411e 100644 --- a/lidar_pipeline/tests/test_export_pdf.py +++ b/lidar_pipeline/tests/test_export_pdf.py @@ -119,3 +119,59 @@ def test_compose_l93_no_data_returns_empty_cells(tmp_path): from lidar_pipeline.export_pdf import compose_l93 img, mask, cells = compose_l93(tmp_path, (0.0, 0.0, 100.0, 100.0), 1.0) assert cells == [] and mask.getextrema() == (0, 0) + + +def test_lab_to_rgb_matches_numpy_reference(): + import numpy as np + from lidar_pipeline.export_pdf import lab_to_rgb + from lidar_pipeline.visualizations import _lab_to_srgb + for L, a, b in [(64, 30, -20), (20, 0, 0), (90, -10, 40), (50, 60, 0)]: + ref = tuple(int(round(v * 255)) for v in _lab_to_srgb( + np.array(L, float), np.array(a, float), np.array(b, float))) + assert all(abs(p - q) <= 1 for p, q in zip(lab_to_rgb(L, a, b), ref)) + + +def test_rose_color_distinct_orientations(): + from lidar_pipeline.export_pdf import rose_color + colors = {rose_color(c) for c in (0, 90, 180, 270)} + assert len(colors) == 4 + + +def test_density_color_classes(): + from lidar_pipeline.export_pdf import DENSITY_CLASSES, density_color + assert density_color(0.0) == DENSITY_CLASSES[0][2] + assert density_color(50.0) == DENSITY_CLASSES[-1][2] + + +def test_pdf_text_replaces_unencodable(): + from lidar_pipeline.export_pdf import _pdf_text + assert _pdf_text("Relief orienté — 1:2 000 ≥ 🚀") == "Relief orienté — 1:2 000 ? ?" + + +def _q(density, start="2023-03-15", end="2023-03-17", empty=0.1): + return {"version": 1, "ground_density": density, + "density_grid": [[density] * 20 for _ in range(20)], + "empty_fraction": empty, "acq_start": start, "acq_end": end, + "acq_source": "gps"} + + +def test_zone_quality_aggregates_two_cells(): + from lidar_pipeline.export_pdf import zone_quality + table = {"LHD_FXX_1054_6882_PTS_LAMB93_IGN69": _q(4.0, "2023-03-15", "2023-03-15"), + "LHD_FXX_1055_6882_PTS_LAMB93_IGN69": _q(8.0, "2023-04-02", "2023-04-03", 0.3)} + bbox = (1054500.0, 6881500.0, 1055500.0, 6881600.0) # moitié de chaque dalle + z = zone_quality(bbox, table, [(1054, 6882), (1055, 6882)]) + assert abs(z["density_mean"] - 6.0) < 1e-6 and z["density_min"] == 4.0 + assert abs(z["empty_fraction"] - 0.2) < 1e-6 + assert (z["acq_start"], z["acq_end"]) == ("2023-03-15", "2023-04-03") + assert z["missing_relief"] == [] and z["missing_quality"] == [] + assert len(z["grid_cells"]) == 2 * 10 * 2 # 10 mailles en x × 2 en y par dalle + + +def test_zone_quality_missing_cells(): + from lidar_pipeline.export_pdf import zone_quality + bbox = (1054500.0, 6881500.0, 1055500.0, 6881600.0) + z = zone_quality(bbox, {}, [(1054, 6882)]) + assert z["density_mean"] is None and z["acq_start"] is None + assert z["missing_quality"] == [(1054, 6882), (1055, 6882)] + assert z["missing_relief"] == [(1055, 6882)]