#!/usr/bin/env python3
"""
Self-contained reproduction of the ODM DEM/mesh slowdown caused by a handful of
gross point-cloud outliers inflating the DEM raster extent.

Needs only what is already inside the opendronemap/odm image (numpy, pdal, and
opendm itself), so a maintainer can run:

    docker run --rm -v "$PWD":/w --entrypoint python3 \
        opendronemap/odm:latest /w/reproduce_dem_outlier_blowup.py

It synthesises a point cloud that mimics a weak-baseline reconstruction: a dense
~80 x 80 m planar scene plus a few dozen stray "flyer" points a couple of km away
(exactly what SfM/MVS produces when scale is poorly constrained). It then runs
the REAL DEM pipeline entry point, opendm.dem.commands.create_dem() -- the same
function odm_dem and odm_meshing call -- twice, with and without the flyers, and
reports the raster size and wall-clock time for each.

NOTE: it deliberately drives the full create_dem chain (renderdem + gdalbuildvrt
+ gdal_translate + gdal_fillnodata + fastrasterfilter), not renderdem alone.
renderdem skips empty tiles quickly for LAS/LAZ input, but the merged
gdal_translate and fastrasterfilter still process the whole inflated raster, so
that is where the time goes.

Expected: WITH flyers -> multi-gigapixel raster, tens of seconds to minutes.
          WITHOUT flyers -> a few-megapixel raster, ~1-2 s.
Same real geometry, driven purely by <0.1% of the points.
"""
import os, sys, json, time, shutil
sys.path.insert(0, "/code")
import numpy as np
import pdal
from opendm.dem import commands

RESOLUTION = 0.05          # DEM resolution in meters
SCENE_M    = 80.0          # real scene is ~80 x 80 m
N_MAIN     = 200_000       # dense points on the real scene
N_FLYERS   = 30            # a tiny number of gross outliers
FLYER_M    = 1500.0        # flyers scattered up to ~1.5 km away


def make_cloud(with_flyers):
    rng = np.random.default_rng(42)
    x = rng.uniform(0, SCENE_M, N_MAIN)
    y = rng.uniform(0, SCENE_M, N_MAIN)
    z = rng.normal(0, 1.5, N_MAIN)
    if with_flyers:
        fx = rng.uniform(-FLYER_M, FLYER_M, N_FLYERS)
        fy = rng.uniform(-FLYER_M, FLYER_M, N_FLYERS)
        fz = rng.uniform(-FLYER_M, FLYER_M, N_FLYERS)
        x = np.concatenate([x, fx]); y = np.concatenate([y, fy]); z = np.concatenate([z, fz])
    arr = np.empty(len(x), dtype=[('X', 'f8'), ('Y', 'f8'), ('Z', 'f8')])
    arr['X'], arr['Y'], arr['Z'] = x, y, z
    return arr


def write_las(arr, path):
    if os.path.exists(path):
        os.remove(path)
    pdal.Pipeline(json.dumps([{"type": "writers.las", "filename": path}]), arrays=[arr]).execute()


def main():
    workdir = "/w" if os.path.isdir("/w") else "."
    print("ODM version:", open("/code/VERSION").read().strip() if os.path.exists("/code/VERSION") else "?")
    print("resolution=%.3fm  scene=%.0fm  main_pts=%d  flyers=%d (up to %.0fm away)\n"
          % (RESOLUTION, SCENE_M, N_MAIN, N_FLYERS, FLYER_M))

    results = {}
    for label, with_flyers in [("WITH flyers", True), ("WITHOUT flyers", False)]:
        pc = os.path.join(workdir, "synthetic_%s.las" % ("flyers" if with_flyers else "clean"))
        arr = make_cloud(with_flyers)
        write_las(arr, pc)
        dx = float(arr['X'].max() - arr['X'].min()); dy = float(arr['Y'].max() - arr['Y'].min())
        mp = dx * dy / (RESOLUTION ** 2) / 1e6
        outdir = os.path.join(workdir, "out")
        if os.path.isdir(outdir):
            shutil.rmtree(outdir)
        os.makedirs(outdir)
        t0 = time.time()
        commands.create_dem(pc, "dsm", output_type="max",
                            radiuses=["0.047", "0.066", "0.094"], gapfill=True,
                            outdir=outdir, resolution=RESOLUTION, max_workers=os.cpu_count(),
                            apply_smoothing=True, max_tiles=None)
        dt = time.time() - t0
        results[label] = (dx, dy, mp, dt)
        print(">> [%s] points=%d  extent=%.0f x %.0f m  raster=%.1f MP  create_dem=%.1f s\n"
              % (label, len(arr), dx, dy, mp, dt))

    a = results["WITH flyers"]; b = results["WITHOUT flyers"]
    print("=" * 66)
    print("SUMMARY  (%d flyers of %d points = %.3f%% of the cloud)"
          % (N_FLYERS, N_MAIN + N_FLYERS, 100.0 * N_FLYERS / (N_MAIN + N_FLYERS)))
    print("  raster : %.1f MP  ->  %.1f MP   (%.0fx smaller without flyers)"
          % (a[2], b[2], a[2] / max(b[2], 1e-9)))
    print("  time   : %.1f s   ->  %.1f s    (%.0fx faster without flyers)"
          % (a[3], b[3], a[3] / max(b[3], 1e-9)))


if __name__ == "__main__":
    main()