script-volume-auto/python/render_images.py

#!/usr/bin/env python3
"""
Generation des images du dossier de resultats:
  - "heightmap" : construit une carte de hauteur (grille 2.5D) coloree en PNG directement a partir
                   d'un nuage ASCII XYZ (binning numpy). Le paquet apt CloudCompare (Debian) plante
                   sur l'export GeoTIFF (assertion GDAL manquante) - on reconstruit donc la grille
                   nous-memes plutot que de dependre de -RASTERIZE -OUTPUT_RASTER_Z.
  - "summary"   : genere les graphiques de synthese (volume/points/robustesse vs niveau de decimation)
                   a partir du report.json produit par le pipeline.

Usage:
    python3 render_images.py heightmap --input cloud.xyz --gridstep 0.05 --output heightmap.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 cmd_heightmap(args: argparse.Namespace) -> int:
    try:
        pts = np.loadtxt(args.input, dtype=np.float64, usecols=(0, 1, 2))
    except Exception as exc:
        print(json.dumps({"error": f"lecture xyz impossible: {exc}"}))
        return 1

    if pts.ndim == 1:
        pts = pts.reshape(1, -1)

    x, y, z = pts[:, 0], pts[:, 1], pts[:, 2]
    xmin, xmax = float(np.min(x)), float(np.max(x))
    ymin, ymax = float(np.min(y)), float(np.max(y))
    nx = max(2, int(np.ceil((xmax - xmin) / args.gridstep)) + 1) if xmax > xmin else 2
    ny = max(2, int(np.ceil((ymax - ymin) / args.gridstep)) + 1) if ymax > ymin else 2
    # Grille bornee pour eviter une image demesuree si le pas est mal renseigne
    nx, ny = min(nx, 2000), min(ny, 2000)

    arr = None
    if pts.shape[0] >= 1 and xmax > xmin and ymax > ymin:
        stat, _, _, _ = binned_statistic_2d(x, y, z, statistic="mean", bins=[nx, ny])
        arr = stat.T

    fig, ax = plt.subplots(figsize=(6, 5.2), dpi=130)
    if arr is None or not np.isfinite(arr).any():
        ax.text(0.5, 0.5, "Aucune donnee exploitable", ha="center", va="center")
    else:
        im = ax.imshow(arr, cmap="terrain", origin="lower", extent=[xmin, xmax, ymin, ymax], aspect="auto")
        cbar = fig.colorbar(im, ax=ax, shrink=0.85)
        cbar.set_label("Altitude (m, moyenne par cellule)")
    ax.set_title(args.title or Path(args.input).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]
    points = [lvl["pointCount"] for lvl in levels]
    matching = [lvl["matchingCellsPct"] for lvl in levels]
    v0 = volumes[0] if volumes else None

    # --- 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)
    ax.semilogy(xs, points, marker="o", color="#16a34a")
    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")
    fig.tight_layout()
    fig.savefig(outdir / "points_vs_level.png")
    plt.close(fig)

    # --- Ecart relatif de volume + robustesse (matching cells %) vs niveau ---
    fig, ax1 = plt.subplots(figsize=(6, 4), dpi=130)
    if v0:
        delta_pct = [100.0 * (v - v0) / v0 for v in volumes]
        ax1.plot(xs, delta_pct, marker="o", color="#dc2626", label="Ecart volume vs niveau 0 (%)")
    ax1.set_xlabel("Niveau de decimation")
    ax1.set_ylabel("Ecart volume (%)", color="#dc2626")
    ax1.set_xticks(xs)
    ax1.grid(alpha=0.3)

    ax2 = ax1.twinx()
    ax2.plot(xs, matching, marker="s", linestyle="--", color="#7c3aed", label="Matching cells (%)")
    ax2.set_ylabel("Matching cells (%)", color="#7c3aed")
    ax2.set_ylim(0, 105)

    fig.suptitle("Robustesse : ecart de volume et couverture de grille")
    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_heightmap = sub.add_parser("heightmap")
    p_heightmap.add_argument("--input", required=True)
    p_heightmap.add_argument("--gridstep", required=True, type=float)
    p_heightmap.add_argument("--output", required=True)
    p_heightmap.add_argument("--title", default=None)
    p_heightmap.set_defaults(func=cmd_heightmap)

    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())