script-volume-auto/python/render_images.py

#!/usr/bin/env python3
"""
Generation des images du dossier de resultats:
  - "diffmap" : construit UNE carte de comparaison (grille 2.5D, degrade bleu->rouge façon
                 CloudCompare) entre le nuage du dessus (ceil) et soit un plan Z constant, soit le
                 nuage du dessous (ground). Le paquet apt CloudCompare (Debian) plante sur l'export
                 GeoTIFF (assertion GDAL manquante) - on reconstruit donc la grille nous-memes a
                 partir d'exports ASCII, plutot que de dependre du rasterizer CloudCompare.
  - "summary" : genere les graphiques de synthese (volume/points vs niveau de decimation) a partir
                 du report.json produit par le pipeline.

Usage:
    python3 render_images.py diffmap --ceil top.xyz --ground bottom.xyz --gridstep 0.05 --output diff.png --title "Niveau 0"
    python3 render_images.py diffmap --ceil cloud.xyz --const-height 100.0 --gridstep 0.05 --output diff.png --title "Niveau 0"
    python3 render_images.py summary --report report.json --outdir images/
"""
import argparse
import json
import sys
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import binned_statistic_2d


def _load_xyz(path: str) -> np.ndarray:
    pts = np.loadtxt(path, dtype=np.float64, usecols=(0, 1, 2))
    if pts.ndim == 1:
        pts = pts.reshape(1, -1)
    return pts


def _grid_edges(xmin: float, xmax: float, ymin: float, ymax: float, step: float):
    nx = max(2, int(np.ceil((xmax - xmin) / step)) + 1) if xmax > xmin else 2
    ny = max(2, int(np.ceil((ymax - ymin) / step)) + 1) if ymax > ymin else 2
    nx, ny = min(nx, 2000), min(ny, 2000)  # borne pour eviter une image demesuree
    return np.linspace(xmin, xmax, nx), np.linspace(ymin, ymax, ny)


def cmd_diffmap(args: argparse.Namespace) -> int:
    try:
        ceil_pts = _load_xyz(args.ceil)
        ground_pts = _load_xyz(args.ground) if args.ground else None
    except Exception as exc:
        print(json.dumps({"error": f"lecture xyz impossible: {exc}"}))
        return 1

    if args.ground is None and args.const_height is None:
        print(json.dumps({"error": "il faut soit --ground, soit --const-height"}))
        return 1

    xmin = float(np.min(ceil_pts[:, 0]))
    xmax = float(np.max(ceil_pts[:, 0]))
    ymin = float(np.min(ceil_pts[:, 1]))
    ymax = float(np.max(ceil_pts[:, 1]))
    if ground_pts is not None:
        xmin = min(xmin, float(np.min(ground_pts[:, 0])))
        xmax = max(xmax, float(np.max(ground_pts[:, 0])))
        ymin = min(ymin, float(np.min(ground_pts[:, 1])))
        ymax = max(ymax, float(np.max(ground_pts[:, 1])))

    fig, ax = plt.subplots(figsize=(6.4, 5.2), dpi=130)
    if xmax <= xmin or ymax <= ymin:
        ax.text(0.5, 0.5, "Aucune donnee exploitable", ha="center", va="center")
    else:
        xedges, yedges = _grid_edges(xmin, xmax, ymin, ymax, args.gridstep)
        ceil_stat, _, _, _ = binned_statistic_2d(ceil_pts[:, 0], ceil_pts[:, 1], ceil_pts[:, 2], statistic="mean", bins=[xedges, yedges])
        if ground_pts is not None:
            ground_stat, _, _, _ = binned_statistic_2d(ground_pts[:, 0], ground_pts[:, 1], ground_pts[:, 2], statistic="mean", bins=[xedges, yedges])
            diff = (ceil_stat - ground_stat).T
        else:
            diff = (ceil_stat - args.const_height).T

        if not np.isfinite(diff).any():
            ax.text(0.5, 0.5, "Aucune donnee exploitable", ha="center", va="center")
        else:
            cmap = plt.get_cmap("jet").copy()
            cmap.set_bad("#e5e7eb")
            im = ax.imshow(diff, cmap=cmap, origin="lower", extent=[xmin, xmax, ymin, ymax], aspect="auto")
            cbar = fig.colorbar(im, ax=ax, shrink=0.85)
            cbar.set_label("Ecart de hauteur (m)")

    ax.set_title(args.title or Path(args.ceil).stem)
    ax.set_xlabel("X (m)")
    ax.set_ylabel("Y (m)")
    fig.tight_layout()
    fig.savefig(args.output)
    plt.close(fig)
    print(json.dumps({"ok": True, "output": args.output}))
    return 0


def cmd_summary(args: argparse.Namespace) -> int:
    with open(args.report, "r", encoding="utf-8") as f:
        report = json.load(f)

    levels = [lvl for lvl in report.get("levels", []) if lvl.get("status") == "ok"]
    if not levels:
        print(json.dumps({"error": "aucun niveau exploitable dans le rapport"}))
        return 1

    outdir = Path(args.outdir)
    outdir.mkdir(parents=True, exist_ok=True)

    xs = [lvl["level"] for lvl in levels]
    volumes = [lvl["volume"] for lvl in levels]
    surfaces = [lvl["surface"] for lvl in levels]
    v0 = volumes[0] if volumes else None
    s0 = surfaces[0] if surfaces else None

    # En mode cloud_compare chaque niveau a 2 nuages (top/bottom) ; en mode const_height, 1 seul.
    point_series = {}
    for lvl in levels:
        for c in lvl.get("clouds", []):
            point_series.setdefault(c["role"], []).append(c["pointCount"])
    role_colors = {"unique": "#16a34a", "top": "#16a34a", "bottom": "#ea580c"}
    role_labels = {"unique": "Points", "top": "Points (haut)", "bottom": "Points (bas)"}

    # --- Volume vs niveau ---
    fig, ax = plt.subplots(figsize=(6, 4), dpi=130)
    ax.plot(xs, volumes, marker="o", color="#2563eb")
    ax.set_xlabel("Niveau de decimation")
    ax.set_ylabel("Volume (m3)")
    ax.set_title("Volume calcule vs niveau de decimation")
    ax.set_xticks(xs)
    ax.grid(alpha=0.3)
    fig.tight_layout()
    fig.savefig(outdir / "volume_vs_level.png")
    plt.close(fig)

    # --- Nombre de points (log) vs niveau ---
    fig, ax = plt.subplots(figsize=(6, 4), dpi=130)
    for role, series in point_series.items():
        ax.semilogy(xs[: len(series)], series, marker="o", color=role_colors.get(role, "#16a34a"), label=role_labels.get(role, role))
    ax.set_xlabel("Niveau de decimation")
    ax.set_ylabel("Nombre de points (log)")
    ax.set_title("Nombre de points vs niveau de decimation")
    ax.set_xticks(xs)
    ax.grid(alpha=0.3, which="both")
    if len(point_series) > 1:
        ax.legend()
    fig.tight_layout()
    fig.savefig(outdir / "points_vs_level.png")
    plt.close(fig)

    # --- Ecart relatif volume + surface vs niveau 0 ---
    fig, ax = plt.subplots(figsize=(6, 4), dpi=130)
    if v0:
        delta_v = [100.0 * (v - v0) / v0 for v in volumes]
        ax.plot(xs, delta_v, marker="o", color="#dc2626", label="Volume")
    if s0:
        delta_s = [100.0 * (s - s0) / s0 for s in surfaces]
        ax.plot(xs, delta_s, marker="s", linestyle="--", color="#0ea5e9", label="Surface")
    ax.set_xlabel("Niveau de decimation")
    ax.set_ylabel("Ecart vs niveau 0 (%)")
    ax.set_title("Robustesse : ecart de volume et de surface vs niveau 0")
    ax.set_xticks(xs)
    ax.grid(alpha=0.3)
    ax.axhline(0, color="#94a3b8", linewidth=1)
    ax.legend()
    fig.tight_layout()
    fig.savefig(outdir / "robustness_vs_level.png")
    plt.close(fig)

    print(json.dumps({"ok": True, "outdir": str(outdir)}))
    return 0


def main() -> int:
    parser = argparse.ArgumentParser()
    sub = parser.add_subparsers(dest="command", required=True)

    p_diffmap = sub.add_parser("diffmap")
    p_diffmap.add_argument("--ceil", required=True)
    p_diffmap.add_argument("--ground", default=None)
    p_diffmap.add_argument("--const-height", dest="const_height", default=None, type=float)
    p_diffmap.add_argument("--gridstep", required=True, type=float)
    p_diffmap.add_argument("--output", required=True)
    p_diffmap.add_argument("--title", default=None)
    p_diffmap.set_defaults(func=cmd_diffmap)

    p_summary = sub.add_parser("summary")
    p_summary.add_argument("--report", required=True)
    p_summary.add_argument("--outdir", required=True)
    p_summary.set_defaults(func=cmd_summary)

    args = parser.parse_args()
    return args.func(args)


if __name__ == "__main__":
    sys.exit(main())