#!/usr/bin/env python3
"""
Convert LAS lidar files to .mound binary format for Three.js rendering.

Usage:
    python las_to_mound.py input.las output.mound
"""

import sys
import struct
import numpy as np
from pathlib import Path

try:
    import laspy
except ImportError:
    print("Error: laspy not installed. Run: pip install laspy")
    sys.exit(1)

try:
    from scipy.spatial import Delaunay
except ImportError:
    print("Error: scipy not installed. Run: pip install scipy")
    sys.exit(1)

try:
    from pyproj import Transformer
except ImportError:
    print("Error: pyproj not installed. Run: pip install pyproj")
    sys.exit(1)


def transform_to_latlon(x, y, z):
    """Transform Ohio State Plane coordinates to lat/lon."""
    print("Transforming coordinates to lat/lon...")
    
    # Ohio State Plane South (EPSG:3735) in US Survey Feet to WGS84 (EPSG:4326)
    # Note: Ohio has two zones - North (3734) and South (3735)
    # Newark is in the North zone
    transformer = Transformer.from_crs("EPSG:3734", "EPSG:4326", always_xy=True)
    
    # Transform x,y (easting, northing) to lon, lat
    lon, lat = transformer.transform(x, y)
    
    # Convert elevation from US Survey Feet to meters
    z_meters = z * 0.3048006096012192
    
    print(f"Transformed to lat/lon bounds: lon[{lon.min():.6f}, {lon.max():.6f}] lat[{lat.min():.6f}, {lat.max():.6f}]")
    
    return lon, lat, z_meters


def read_las(filepath):
    """Read LAS file and extract ground points."""
    print(f"Reading {filepath}...")
    las = laspy.read(filepath)
    
    # Extract coordinates
    x = las.x
    y = las.y
    z = las.z
    
    # Filter for ground points (classification 2) if available
    if hasattr(las, 'classification'):
        ground_mask = las.classification == 2
        if ground_mask.any():
            print(f"Filtering {ground_mask.sum()} ground points from {len(x)} total points")
            x = x[ground_mask]
            y = y[ground_mask]
            z = z[ground_mask]
    
    print(f"Loaded {len(x)} points")
    return x, y, z


def triangulate_points(x, y, z):
    """Perform Delaunay triangulation on XY coordinates."""
    print("Performing Delaunay triangulation...")
    
    # Create 2D point array for triangulation
    points_2d = np.column_stack([x, y])
    
    # Triangulate
    tri = Delaunay(points_2d)
    
    print(f"Generated {len(tri.simplices)} triangles")
    
    return tri.simplices


def write_mound(filepath, x, y, z, indices):
    """Write .mound binary format."""
    print(f"Writing {filepath}...")
    
    # Convert to float32
    positions = np.column_stack([x, y, z]).astype(np.float32)
    indices = indices.astype(np.uint32)
    
    # Calculate bounds
    min_x, min_y, min_z = positions.min(axis=0)
    max_x, max_y, max_z = positions.max(axis=0)
    
    point_count = len(positions)
    triangle_count = len(indices)
    
    with open(filepath, 'wb') as f:
        # Header (64 bytes)
        f.write(b'LIDR')  # Magic number (4 bytes)
        f.write(struct.pack('I', 1))  # Version (4 bytes)
        f.write(struct.pack('I', point_count))  # Point count (4 bytes)
        f.write(struct.pack('I', triangle_count))  # Triangle count (4 bytes)
        f.write(struct.pack('f', min_x))  # Min X (4 bytes)
        f.write(struct.pack('f', min_y))  # Min Y (4 bytes)
        f.write(struct.pack('f', min_z))  # Min Z (4 bytes)
        f.write(struct.pack('f', max_x))  # Max X (4 bytes)
        f.write(struct.pack('f', max_y))  # Max Y (4 bytes)
        f.write(struct.pack('f', max_z))  # Max Z (4 bytes)
        f.write(b'\x00' * 24)  # Reserved (24 bytes)
        
        # Position data
        f.write(positions.tobytes())
        
        # Index data
        f.write(indices.tobytes())
    
    file_size = Path(filepath).stat().st_size / (1024 * 1024)
    print(f"Wrote {point_count} points, {triangle_count} triangles ({file_size:.2f} MB)")
    print(f"Bounds: X[{min_x:.2f}, {max_x:.2f}] Y[{min_y:.2f}, {max_y:.2f}] Z[{min_z:.2f}, {max_z:.2f}]")


def main():
    if len(sys.argv) != 3:
        print("Usage: python las_to_mound.py input.las output.mound")
        sys.exit(1)
    
    input_file = sys.argv[1]
    output_file = sys.argv[2]
    
    if not Path(input_file).exists():
        print(f"Error: Input file '{input_file}' not found")
        sys.exit(1)
    
    # Read LAS
    x, y, z = read_las(input_file)
    
    # Transform to lat/lon
    lon, lat, z_meters = transform_to_latlon(x, y, z)
    
    # Triangulate (using lon/lat as x/y)
    indices = triangulate_points(lon, lat, z_meters)
    
    # Write output (lon as x, lat as y, elevation as z)
    write_mound(output_file, lon, lat, z_meters, indices)
    
    print("Done!")


if __name__ == '__main__':
    main()