点云-俯视图

"""
Point Cloud Ground Plan → Top-Down Image
=========================================
Extract a thin ground slice from a point cloud, voxel-downsample to balance
station density, and render an XY plan-view image with per-point RGB colours.

Supports: .e57, .las, .laz, .ply, .xyz, .txt, .asc, .csv

Usage:
    python bottom-generator.py --input scan.e57
    python bottom-generator.py --input warehouse.e57 --resolution 0.02
    python bottom-generator.py --input scan.e57 --mode band --z-min 0 --z-max 2

Ground mode (default):
    Extracts a thin slice at the cloud floor: Z = [floor, floor + slice_thickness].
    Default slice_thickness = 5 cm — keeps only floor returns, drops shelf/object noise.

Canvas sizing:
    Width/height are derived from XY extent ÷ --resolution (metres per pixel).
    E.g. 147 m span at 2 cm/px → ~7350 px.  Override with --width if needed.
"""

from __future__ import annotations

import argparse
import json
import math
import pathlib
import sys

import numpy as np

# ---------------------------------------------------------------------------
# PARAMETERS  (edit defaults or override via CLI)
# ---------------------------------------------------------------------------
INPUT_PATH        = ""
OUTPUT_PATH       = ""
MODE              = "ground"       # ground | band
Z_MIN             = 0.0            # metres above floor (band mode)
Z_MAX             = 0.05           # metres above floor (band mode / slice thickness)
SLICE_THICKNESS   = 0.05           # ground mode: metres above floor
RESOLUTION        = 0.02           # metres per pixel (2 cm)
MAX_DIMENSION     = 16384          # cap longest canvas side (0 = no cap)
VOXEL_SIZE        = 0.02           # voxel grid cell (m); 0 = disabled
IMAGE_WIDTH       = 0              # 0 = auto from resolution
POINT_SIZE        = 1
MAX_POINTS        = 0              # random subsample after voxel (0 = all)
BACKGROUND        = (255, 255, 255)
MARGIN            = 0.02           # fraction of XY extent added as border
# ---------------------------------------------------------------------------

try:
    import laspy
    from PIL import Image
except ImportError as e:
    sys.exit(
        f"Missing dependency: {e}\n"
        "Install with:  pip install 'laspy[lazrs]' pye57 numpy Pillow"
    )

try:
    import pye57

    _HAS_PYE57 = True
except ImportError:
    _HAS_PYE57 = False

try:
    import open3d as o3d

    _HAS_OPEN3D = True
except ImportError:
    _HAS_OPEN3D = False


# ───────────────────────────────────────────────────────────────────────────
# Load point cloud
# ───────────────────────────────────────────────────────────────────────────
def _apply_pose(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    translation: np.ndarray,
    rotation_q: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Apply E57 scan pose (quaternion [w,x,y,z] + translation)."""
    trans = np.asarray(translation, dtype=np.float64)
    rot_q = np.asarray(rotation_q, dtype=np.float64)
    is_identity_rot = np.allclose(rot_q, [1, 0, 0, 0], atol=1e-6)
    is_zero_trans = np.allclose(trans, [0, 0, 0], atol=1e-6)
    if is_identity_rot and is_zero_trans:
        return x, y, z

    w, qx, qy, qz = rot_q
    r = np.array([
        [1 - 2 * (qy * qy + qz * qz), 2 * (qx * qy - qz * w), 2 * (qx * qz + qy * w)],
        [2 * (qx * qy + qz * w), 1 - 2 * (qx * qx + qz * qz), 2 * (qy * qz - qx * w)],
        [2 * (qx * qz - qy * w), 2 * (qy * qz + qx * w), 1 - 2 * (qx * qx + qy * qy)],
    ])
    pts = r @ np.vstack([x, y, z]) + trans[:, None]
    return pts[0], pts[1], pts[2]


def _should_apply_pose(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    translation: np.ndarray,
    rotation_q: np.ndarray,
) -> bool:
    """
    Skip pose when cartesian* already stores global coordinates.
    Cyclone-registered multi-scan E57 often has points aligned with
    pose.translation — re-applying rotation would double-transform.
    """
    trans = np.asarray(translation, dtype=np.float64)
    rot_q = np.asarray(rotation_q, dtype=np.float64)
    if np.allclose(rot_q, [1, 0, 0, 0], atol=1e-6) and np.allclose(trans, [0, 0, 0], atol=1e-6):
        return False

    centroid_xy = np.array([float(x.mean()), float(y.mean())])
    if np.linalg.norm(centroid_xy - trans[:2]) < 2.0:
        return False
    return True


def _normalize_rgb(r: np.ndarray, g: np.ndarray, b: np.ndarray) -> np.ndarray:
    """Return (N, 3) uint8 RGB array."""
    rgb = np.stack([r, g, b], axis=1).astype(np.float64)
    if rgb.size == 0:
        return rgb.astype(np.uint8)
    peak = float(rgb.max())
    if peak <= 1.0:
        rgb *= 255.0
    return np.clip(rgb, 0, 255).astype(np.uint8)


def load_las(path: str) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]:
    las = laspy.read(path)
    x = np.asarray(las.x, dtype=np.float64)
    y = np.asarray(las.y, dtype=np.float64)
    z = np.asarray(las.z, dtype=np.float64)
    rgb = None
    dims = {d.name.lower(): d.name for d in las.point_format.dimensions}
    if all(k in dims for k in ("red", "green", "blue")):
        rgb = _normalize_rgb(
            np.asarray(las[dims["red"]]),
            np.asarray(las[dims["green"]]),
            np.asarray(las[dims["blue"]]),
        )
    return x, y, z, rgb


def load_xyz(path: str) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]:
    data = np.loadtxt(path, comments=["#", "//"], usecols=(0, 1, 2), dtype=np.float64)
    if data.ndim == 1:
        data = data[np.newaxis, :]
    return data[:, 0], data[:, 1], data[:, 2], None


def load_e57(path: str) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]:
    if not _HAS_PYE57:
        sys.exit("pye57 is required for E57 files.  pip install pye57")

    e = pye57.E57(path)
    xs_all, ys_all, zs_all, rgb_all = [], [], [], []
    has_color = False

    for scan_idx in range(e.scan_count):
        h = e.get_header(scan_idx)
        data = e.read_scan(
            scan_idx,
            intensity=False,
            colors=True,
            row_column=False,
            ignore_missing_fields=True,
        )
        x = np.asarray(data["cartesianX"], dtype=np.float64)
        y = np.asarray(data["cartesianY"], dtype=np.float64)
        z = np.asarray(data["cartesianZ"], dtype=np.float64)
        if _should_apply_pose(x, y, z, h.translation, h.rotation):
            x, y, z = _apply_pose(x, y, z, h.translation, h.rotation)
        else:
            print(f"      Scan {scan_idx}: global coords (pose skipped)", flush=True)

        scan_rgb = None
        if all(k in data for k in ("colorRed", "colorGreen", "colorBlue")):
            has_color = True
            scan_rgb = _normalize_rgb(
                np.asarray(data["colorRed"]),
                np.asarray(data["colorGreen"]),
                np.asarray(data["colorBlue"]),
            )

        xs_all.append(x)
        ys_all.append(y)
        zs_all.append(z)
        if scan_rgb is not None:
            rgb_all.append(scan_rgb)
        print(f"      Scan {scan_idx}: {len(x):,} points", flush=True)

    rgb = np.concatenate(rgb_all) if has_color else None
    return np.concatenate(xs_all), np.concatenate(ys_all), np.concatenate(zs_all), rgb


def load_ply(path: str) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]:
    if not _HAS_OPEN3D:
        sys.exit("open3d is required for PLY files.  pip install open3d")

    pcd = o3d.io.read_point_cloud(path)
    pts = np.asarray(pcd.points, dtype=np.float64)
    if len(pts) == 0:
        sys.exit(f"No points found in {path}")
    rgb = None
    if pcd.has_colors():
        rgb = _normalize_rgb(
            np.asarray(pcd.colors[:, 0]),
            np.asarray(pcd.colors[:, 1]),
            np.asarray(pcd.colors[:, 2]),
        )
    return pts[:, 0], pts[:, 1], pts[:, 2], rgb


def load_pointcloud(path: str) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray | None]:
    print(f"\nLoading {path} …", flush=True)
    ext = pathlib.Path(path).suffix.lower()
    if ext == ".e57":
        x, y, z, rgb = load_e57(path)
    elif ext in (".las", ".laz"):
        x, y, z, rgb = load_las(path)
    elif ext == ".ply":
        x, y, z, rgb = load_ply(path)
    elif ext in (".xyz", ".txt", ".asc", ".csv"):
        x, y, z, rgb = load_xyz(path)
    else:
        sys.exit(
            f"Unsupported file format: {ext}\n"
            "Expected: .las, .laz, .e57, .ply, .xyz, .txt, .asc, .csv"
        )

    print(
        f"  total  : {len(x):,} points\n"
        f"  X range: [{x.min():.3f}, {x.max():.3f}]  span {x.max()-x.min():.2f} m\n"
        f"  Y range: [{y.min():.3f}, {y.max():.3f}]  span {y.max()-y.min():.2f} m\n"
        f"  Z range: [{z.min():.3f}, {z.max():.3f}]\n"
        f"  colour : {'yes' if rgb is not None else 'no (height colormap will be used)'}"
    )
    return x, y, z, rgb


# ───────────────────────────────────────────────────────────────────────────
# Filter + downsample + render
# ───────────────────────────────────────────────────────────────────────────
def height_colormap(z: np.ndarray, z_lo: float, z_hi: float) -> np.ndarray:
    """Blue→green→red gradient by height within the band."""
    span = max(z_hi - z_lo, 1e-9)
    t = np.clip((z - z_lo) / span, 0.0, 1.0)
    r = (np.clip(2 * t - 0.5, 0, 1) * 255).astype(np.uint8)
    g = (np.clip(1 - 2 * np.abs(t - 0.5), 0, 1) * 255).astype(np.uint8)
    b = (np.clip(1.5 - 2 * t, 0, 1) * 255).astype(np.uint8)
    return np.stack([r, g, b], axis=1)


def estimate_floor_z(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    grid_size: float | None = None,
    min_peak_frac: float = 0.05,
) -> float:
    """
    Estimate floor height from per-cell minimum Z histogram.
    Finds the lowest significant peak — robust against outliers below the
    true floor and against ceiling/shelf peaks above it.
    """
    x_span = float(x.max() - x.min())
    y_span = float(y.max() - y.min())
    if grid_size is None:
        grid_size = max(0.25, min(x_span, y_span) * 0.01)

    origin = np.array([x.min(), y.min()], dtype=np.float64)
    ij = (np.floor((np.column_stack([x, y]) - origin) / grid_size)).astype(np.int64)
    ny = int(ij[:, 1].max()) + 1
    keys = ij[:, 0].astype(np.int64) * ny + ij[:, 1]
    n_cells = int(keys.max()) + 1
    cell_min_z = np.full(n_cells, np.inf, dtype=np.float64)
    np.minimum.at(cell_min_z, keys, z)
    valid = cell_min_z[np.isfinite(cell_min_z)]
    if len(valid) < 10:
        return float(z.min())

    z_lo = float(valid.min())
    z_hi = float(np.percentile(valid, 95))
    n_bins = max(50, min(300, int((z_hi - z_lo) / 0.02)))
    hist, edges = np.histogram(valid, bins=n_bins, range=(z_lo, z_hi))
    min_count = max(5, int(hist.max() * min_peak_frac))

    peaks: list[int] = []
    for i in range(1, len(hist) - 1):
        if hist[i] >= hist[i - 1] and hist[i] > hist[i + 1] and hist[i] >= min_count:
            peaks.append(i)

    if not peaks:
        floor = float(np.median(valid))
    else:
        floor = 0.5 * (edges[peaks[0]] + edges[peaks[0] + 1])

    print(
        f"  floor   : {floor:.4f}m  (grid={grid_size:.2f}m, "
        f"{len(valid):,} cells, z.min={z.min():.4f}m)",
        flush=True,
    )
    return floor


def resolve_height_band(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    mode: str,
    z_min: float,
    z_max: float,
    slice_thickness: float,
    relative: bool,
    floor_grid: float,
) -> tuple[float, float]:
    """Return absolute (z_lo, z_hi) for the extraction band."""
    if mode == "ground":
        if relative:
            floor = estimate_floor_z(x, y, z, grid_size=floor_grid if floor_grid > 0 else None)
            abs_lo = floor + z_min
            abs_hi = floor + z_min + slice_thickness
            print(
                f"  ground slice: floor + [{z_min}, {z_min + slice_thickness}]m  "
                f"→ Z=[{abs_lo:.4f}, {abs_hi:.4f}]m  (thickness={slice_thickness*100:.0f} cm)"
            )
        else:
            abs_lo = z_min
            abs_hi = z_min + slice_thickness
            print(
                f"  ground slice (absolute): Z=[{abs_lo:.4f}, {abs_hi:.4f}]m  "
                f"(thickness={slice_thickness*100:.0f} cm)"
            )
    else:
        if relative:
            floor = estimate_floor_z(x, y, z, grid_size=floor_grid if floor_grid > 0 else None)
            abs_lo = floor + z_min
            abs_hi = floor + z_max
            print(
                f"  height band (relative): detected floor({floor:.3f}m) + "
                f"[{z_min}, {z_max}]m → Z=[{abs_lo:.4f}, {abs_hi:.4f}]m"
            )
        else:
            abs_lo, abs_hi = z_min, z_max
            print(f"  height band (absolute): Z=[{abs_lo:.4f}, {abs_hi:.4f}]m")

    if abs_lo > abs_hi:
        abs_lo, abs_hi = abs_hi, abs_lo
    return abs_lo, abs_hi


def filter_height_band(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    rgb: np.ndarray | None,
    abs_lo: float,
    abs_hi: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    mask = (z >= abs_lo) & (z <= abs_hi)
    xs, ys, zs = x[mask], y[mask], z[mask]
    if len(xs) == 0:
        sys.exit("No points found in the requested height band.")

    if rgb is not None:
        colours = rgb[mask]
    else:
        colours = height_colormap(zs, abs_lo, abs_hi)

    print(f"  filtered: {len(xs):,} points  ({100 * mask.mean():.2f}% of cloud)")
    return xs, ys, zs, colours


def voxel_downsample(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    rgb: np.ndarray,
    voxel_size: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Voxel-grid filter: one representative point per cell with averaged colour.
    Balances density between dense scanner stations and sparse distant areas.
    """
    n = len(x)
    if voxel_size <= 0 or n == 0:
        return x, y, z, rgb

    origin = np.array([x.min(), y.min(), z.min()], dtype=np.float64)
    ijk = (np.floor((np.column_stack([x, y, z]) - origin) / voxel_size)).astype(np.int64)

    nx = int(ijk[:, 0].max()) + 1
    ny = int(ijk[:, 1].max()) + 1
    nz = int(ijk[:, 2].max()) + 1
    keys = (
        ijk[:, 0].astype(np.int64) * (ny * nz)
        + ijk[:, 1].astype(np.int64) * nz
        + ijk[:, 2]
    )

    _, inverse = np.unique(keys, return_inverse=True)
    counts = np.bincount(inverse, minlength=inverse.max() + 1).astype(np.float64)

    x_out = np.bincount(inverse, weights=x) / counts
    y_out = np.bincount(inverse, weights=y) / counts
    z_out = np.bincount(inverse, weights=z) / counts
    r_out = np.bincount(inverse, weights=rgb[:, 0].astype(np.float64)) / counts
    g_out = np.bincount(inverse, weights=rgb[:, 1].astype(np.float64)) / counts
    b_out = np.bincount(inverse, weights=rgb[:, 2].astype(np.float64)) / counts
    rgb_out = np.stack([r_out, g_out, b_out], axis=1).clip(0, 255).astype(np.uint8)

    print(
        f"  voxel   : {n:,} → {len(x_out):,} pts  "
        f"(cell={voxel_size*100:.1f} cm, avg {n/len(x_out):.1f} pts/cell)",
        flush=True,
    )
    return x_out, y_out, z_out, rgb_out


def subsample_points(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray,
    rgb: np.ndarray,
    max_points: int,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    if max_points <= 0 or len(x) <= max_points:
        return x, y, z, rgb
    rng = np.random.default_rng(42)
    idx = rng.choice(len(x), size=max_points, replace=False)
    print(f"  subsample: {max_points:,} / {len(x):,} points", flush=True)
    return x[idx], y[idx], z[idx], rgb[idx]


def compute_canvas_size(
    x: np.ndarray,
    y: np.ndarray,
    resolution: float,
    margin: float,
    max_dimension: int,
    fixed_width: int,
    bounds: tuple[float, float, float, float] | None = None,
) -> tuple[int, int, float, float, float, float]:
    """
    Compute image dimensions from physical extent.
    Returns (width, height, xmin, xmax, ymin, ymax) with margin applied.
    """
    if bounds is not None:
        xmin, xmax, ymin, ymax = bounds
    else:
        xmin, xmax = float(x.min()), float(x.max())
        ymin, ymax = float(y.min()), float(y.max())
    dx = max(xmax - xmin, 1e-9)
    dy = max(ymax - ymin, 1e-9)
    pad_x = dx * margin
    pad_y = dy * margin
    xmin -= pad_x
    xmax += pad_x
    ymin -= pad_y
    ymax += pad_y
    dx = xmax - xmin
    dy = ymax - ymin

    if fixed_width > 0:
        width = fixed_width
        height = max(1, int(round(width * dy / dx)))
        effective_res = dx / max(width - 1, 1)
        print(
            f"  canvas : {width}×{height} px  (fixed width, ~{effective_res*100:.2f} cm/px)",
            flush=True,
        )
    else:
        width = max(64, int(math.ceil(dx / resolution)))
        height = max(64, int(math.ceil(dy / resolution)))
        if max_dimension > 0:
            longest = max(width, height)
            if longest > max_dimension:
                scale = max_dimension / longest
                width = max(64, int(round(width * scale)))
                height = max(64, int(round(height * scale)))
                effective = dx / max(width - 1, 1)
                print(
                    f"  canvas : {width}×{height} px  "
                    f"(capped at {max_dimension}, effective ~{effective*100:.2f} cm/px)",
                    flush=True,
                )
            else:
                print(
                    f"  canvas : {width}×{height} px  "
                    f"({resolution*100:.1f} cm/px, extent {dx:.1f}×{dy:.1f} m)",
                    flush=True,
                )
        else:
            print(
                f"  canvas : {width}×{height} px  "
                f"({resolution*100:.1f} cm/px, extent {dx:.1f}×{dy:.1f} m)",
                flush=True,
            )

    return width, height, xmin, xmax, ymin, ymax


def render_plan_image(
    x: np.ndarray,
    y: np.ndarray,
    rgb: np.ndarray,
    width: int,
    height: int,
    xmin: float,
    xmax: float,
    ymin: float,
    ymax: float,
    point_size: int,
    background: tuple[int, int, int] | None,
) -> Image.Image:
    dx = max(xmax - xmin, 1e-9)
    dy = max(ymax - ymin, 1e-9)

    px = np.floor((x - xmin) / dx * (width - 1)).astype(np.int32)
    py = np.floor((ymax - y) / dy * (height - 1)).astype(np.int32)
    px = np.clip(px, 0, width - 1)
    py = np.clip(py, 0, height - 1)

    transparent = background is None
    if transparent:
        img = np.zeros((height, width, 4), dtype=np.uint8)
    else:
        img = np.full((height, width, 4), (*background, 255), dtype=np.uint8)

    filled = np.zeros((height, width), dtype=bool)

    if point_size <= 1:
        img[py, px, :3] = rgb
        if transparent:
            img[py, px, 3] = 255
        filled[py, px] = True
    else:
        half = point_size // 2
        for dy_off in range(-half, half + 1):
            row = np.clip(py + dy_off, 0, height - 1)
            for dx_off in range(-half, half + 1):
                col = np.clip(px + dx_off, 0, width - 1)
                img[row, col, :3] = rgb
                if transparent:
                    img[row, col, 3] = 255
                filled[row, col] = True

    n_filled = int(filled.sum())
    coverage = 100.0 * n_filled / (width * height)
    print(f"  coverage: {n_filled:,} / {width*height:,} px ({coverage:.1f}%)", flush=True)

    return Image.fromarray(img, mode="RGBA")


def parse_background(value: str) -> tuple[int, int, int] | None:
    if value.strip().lower() in ("transparent", "none", "alpha"):
        return None
    parts = value.split(",")
    if len(parts) != 3:
        sys.exit("--background must be transparent, or R,G,B with three integers 0-255")
    try:
        rgb = tuple(int(p.strip()) for p in parts)
    except ValueError:
        sys.exit("--background must be transparent, or R,G,B with three integers 0-255")
    if any(c < 0 or c > 255 for c in rgb):
        sys.exit("--background values must be in [0, 255]")
    return rgb  # type: ignore[return-value]


def default_output_path(
    input_path: str,
    mode: str,
    z_min: float,
    z_max: float,
    slice_thickness: float,
    resolution: float,
) -> pathlib.Path:
    stem = pathlib.Path(input_path).stem
    parent = pathlib.Path(input_path).parent
    if mode == "ground":
        tag = f"ground_{slice_thickness*100:.0f}cm_{resolution*100:.0f}cmpx"
    else:
        tag = f"z{z_min:g}-{z_max:g}m_{resolution*100:.0f}cmpx".replace(".", "p")
    return parent / f"{stem}_plan_{tag}.png"


# ───────────────────────────────────────────────────────────────────────────
# CLI
# ───────────────────────────────────────────────────────────────────────────
def build_parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(
        description="Render a top-down ground plan from a point cloud.",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=(
            "Examples:\n"
            "  python bottom-generator.py --input warehouse.e57\n"
            "  python bottom-generator.py --input scan.e57 --resolution 0.01\n"
            "  python bottom-generator.py --input scan.e57 --mode band --z-min 0 --z-max 2 --transparent\n"
            "  python bottom-generator.py --input scan.e57 --slice-thickness 0.03 --voxel 0.02"
        ),
    )
    parser.add_argument(
        "--input", required=True,
        help="Input point cloud (.las/.laz/.e57/.ply/.xyz/.txt/.asc/.csv)",
    )
    parser.add_argument(
        "--output", default=None,
        help="Output PNG path (default: auto-named)",
    )
    parser.add_argument(
        "--mode", choices=("ground", "band"), default=MODE,
        help="ground = thin floor slice (default); band = arbitrary height range",
    )
    parser.add_argument(
        "--z-min", type=float, default=Z_MIN,
        help="Lower bound above floor (m); ground mode uses as floor offset",
    )
    parser.add_argument(
        "--z-max", type=float, default=Z_MAX,
        help="Upper bound above floor (m); band mode only",
    )
    parser.add_argument(
        "--slice-thickness", type=float, default=SLICE_THICKNESS,
        help=f"Ground slice thickness in metres (default: {SLICE_THICKNESS} = 5 cm)",
    )
    parser.add_argument(
        "--absolute-z", action="store_true",
        help="Treat height bounds as absolute survey Z",
    )
    parser.add_argument(
        "--resolution", type=float, default=RESOLUTION,
        help=f"Metres per pixel for auto canvas sizing (default: {RESOLUTION} = 2 cm/px)",
    )
    parser.add_argument(
        "--width", type=int, default=IMAGE_WIDTH,
        help="Fixed image width in pixels (0 = auto from --resolution)",
    )
    parser.add_argument(
        "--max-dimension", type=int, default=MAX_DIMENSION,
        help=f"Cap longest canvas side when auto-sizing (default: {MAX_DIMENSION}, 0=no cap)",
    )
    parser.add_argument(
        "--voxel", type=float, default=VOXEL_SIZE,
        help=f"Voxel grid size in metres before render (default: {VOXEL_SIZE}, 0=disabled)",
    )
    parser.add_argument(
        "--no-voxel", action="store_true",
        help="Disable voxel downsampling",
    )
    parser.add_argument(
        "--point-size", type=int, default=POINT_SIZE,
        help=f"Draw each point as an N×N pixel block (default: {POINT_SIZE})",
    )
    parser.add_argument(
        "--max-points", type=int, default=MAX_POINTS,
        help="Random subsample cap after voxel filter (0 = all)",
    )
    parser.add_argument(
        "--background", default="255,255,255",
        help="Background as R,G,B, or 'transparent' (default: white)",
    )
    parser.add_argument(
        "--transparent", action="store_true",
        help="Shorthand for --background transparent",
    )
    parser.add_argument(
        "--floor-grid", type=float, default=0.0,
        help="XY grid size (m) for floor detection (0 = auto from scene extent)",
    )
    parser.add_argument(
        "--margin", type=float, default=MARGIN,
        help=f"XY border margin as fraction of extent (default: {MARGIN})",
    )
    return parser


def main() -> None:
    args = build_parser().parse_args()
    input_path = pathlib.Path(args.input)
    if not input_path.exists():
        sys.exit(f"Input file not found: {input_path}")

    if args.resolution <= 0 and args.width <= 0:
        sys.exit("Set --resolution > 0 or --width > 0")
    if args.width > 0 and args.width < 64:
        sys.exit("--width must be >= 64")
    if args.point_size < 1:
        sys.exit("--point-size must be >= 1")
    if args.slice_thickness <= 0:
        sys.exit("--slice-thickness must be > 0")

    voxel_size = 0.0 if args.no_voxel else args.voxel

    x, y, z, rgb = load_pointcloud(str(input_path))
    xy_bounds = (float(x.min()), float(x.max()), float(y.min()), float(y.max()))
    abs_lo, abs_hi = resolve_height_band(
        x, y, z, args.mode, args.z_min, args.z_max,
        args.slice_thickness, relative=not args.absolute_z,
        floor_grid=args.floor_grid,
    )
    xs, ys, zs, colours = filter_height_band(x, y, z, rgb, abs_lo, abs_hi)

    if voxel_size > 0:
        xs, ys, zs, colours = voxel_downsample(xs, ys, zs, colours, voxel_size)

    xs, ys, zs, colours = subsample_points(xs, ys, zs, colours, args.max_points)

    width, height, xmin, xmax, ymin, ymax = compute_canvas_size(
        xs, ys,
        resolution=args.resolution,
        margin=args.margin,
        max_dimension=args.max_dimension,
        fixed_width=args.width,
        bounds=xy_bounds,
    )

    image = render_plan_image(
        xs, ys, colours,
        width=width,
        height=height,
        xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
        point_size=args.point_size,
        background=parse_background("transparent" if args.transparent else args.background),
    )

    out_path = pathlib.Path(args.output) if args.output else default_output_path(
        str(input_path), args.mode, args.z_min, args.z_max,
        args.slice_thickness, args.resolution,
    )
    out_path.parent.mkdir(parents=True, exist_ok=True)
    image.save(out_path, format="PNG", optimize=True)

    bounds_meta = {
        "minX": xmin,
        "maxX": xmax,
        "minY": ymin,
        "maxY": ymax,
        "widthPx": width,
        "heightPx": height,
        "resolution": args.resolution if args.width <= 0 else None,
        "margin": args.margin,
        "renderBounds": True,
        "rawBounds": {
            "minX": xy_bounds[0],
            "maxX": xy_bounds[1],
            "minY": xy_bounds[2],
            "maxY": xy_bounds[3],
        },
    }
    bounds_path = out_path.with_suffix(".bounds.json")
    with open(bounds_path, "w", encoding="utf-8") as f:
        json.dump(bounds_meta, f, indent=2)
    print(f"  bounds  → {bounds_path}", flush=True)
    print(f"\nSaved → {out_path}  ({out_path.stat().st_size / 1024 / 1024:.2f} MB)")


if __name__ == "__main__":
    main()