"""Enhanced building export with point cloud roof extraction.

Combines CityGML footprints with BDOM20RGBI/LPO point clouds for:
1. RANSAC-based roof plane extraction (where BDOM available)
2. LPO height refinement
3. DOP20 facade texturing
4. Fallback to standard CityGML pipeline
"""

from __future__ import annotations

import csv
import glob
import os
from dataclasses import dataclass
from typing import List, Optional, Tuple

import numpy as np
from osgeo import gdal

from .config import Config
from .gdal_utils import open_dataset
from .pointcloud import (
    PointCloud,
    find_laz_file,
    find_pointcloud_file,
    has_bdom_data,
    has_lpo_data,
    read_laz_file,
    read_pointcloud_file,
    sample_rgb_average,
)

try:
    from shapely.geometry import Polygon
    from shapely.ops import unary_union
    HAS_SHAPELY = True
except ImportError:
    HAS_SHAPELY = False


@dataclass
class RoofPlane:
    """A detected roof plane from RANSAC."""
    normal: np.ndarray  # Unit normal vector
    d: float            # Plane distance from origin
    inliers: np.ndarray  # Point indices
    centroid: np.ndarray  # Plane centroid


def ransac_plane_fit(
    points: np.ndarray,
    max_iterations: int = 1000,
    distance_threshold: float = 0.15,
    min_inliers: int = 50,
) -> Optional[RoofPlane]:
    """Fit a single plane using RANSAC.

    Args:
        points: Nx3 array of points
        max_iterations: Maximum RANSAC iterations
        distance_threshold: Inlier distance threshold in meters
        min_inliers: Minimum inliers for valid plane

    Returns:
        RoofPlane or None if no plane found
    """
    if len(points) < 3:
        return None

    best_inliers = None
    best_normal = None
    best_d = None

    n_points = len(points)

    for _ in range(max_iterations):
        # Sample 3 random points
        idx = np.random.choice(n_points, 3, replace=False)
        p1, p2, p3 = points[idx]

        # Compute plane normal
        v1 = p2 - p1
        v2 = p3 - p1
        normal = np.cross(v1, v2)

        norm_len = np.linalg.norm(normal)
        if norm_len < 1e-10:
            continue

        normal = normal / norm_len

        # Plane equation: normal . (x - p1) = 0
        # => normal . x = normal . p1 = d
        d = np.dot(normal, p1)

        # Compute distances to plane
        distances = np.abs(np.dot(points, normal) - d)

        # Find inliers
        inlier_mask = distances < distance_threshold
        inlier_count = np.sum(inlier_mask)

        if inlier_count >= min_inliers:
            if best_inliers is None or inlier_count > len(best_inliers):
                best_inliers = np.where(inlier_mask)[0]
                best_normal = normal
                best_d = d

    if best_inliers is None:
        return None

    # Refit plane using all inliers (least squares)
    inlier_points = points[best_inliers]
    centroid = np.mean(inlier_points, axis=0)

    # SVD to find best-fit plane
    centered = inlier_points - centroid
    _, _, vh = np.linalg.svd(centered)
    refined_normal = vh[2]  # Smallest singular value direction

    # Ensure normal points upward
    if refined_normal[2] < 0:
        refined_normal = -refined_normal

    refined_d = np.dot(refined_normal, centroid)

    return RoofPlane(
        normal=refined_normal,
        d=refined_d,
        inliers=best_inliers,
        centroid=centroid,
    )


def extract_roof_planes(
    points: PointCloud,
    max_planes: int = 6,
    ransac_threshold: float = 0.15,
    ransac_iterations: int = 1000,
    min_inliers: int = 50,
) -> List[RoofPlane]:
    """Extract multiple roof planes using iterative RANSAC.

    Args:
        points: Point cloud of roof points
        max_planes: Maximum number of planes to extract
        ransac_threshold: Distance threshold for inliers
        ransac_iterations: RANSAC iterations per plane
        min_inliers: Minimum points for valid plane

    Returns:
        List of extracted RoofPlane objects
    """
    if len(points) < min_inliers:
        return []

    xyz = points.xyz.copy()
    planes = []
    remaining_mask = np.ones(len(xyz), dtype=bool)

    for _ in range(max_planes):
        remaining_points = xyz[remaining_mask]

        if len(remaining_points) < min_inliers:
            break

        plane = ransac_plane_fit(
            remaining_points,
            max_iterations=ransac_iterations,
            distance_threshold=ransac_threshold,
            min_inliers=min_inliers,
        )

        if plane is None:
            break

        # Map inlier indices back to original array
        remaining_indices = np.where(remaining_mask)[0]
        original_inliers = remaining_indices[plane.inliers]
        plane.inliers = original_inliers

        planes.append(plane)

        # Remove inliers from remaining points
        remaining_mask[original_inliers] = False

    return planes


def classify_roof_type(planes: List[RoofPlane]) -> str:
    """Classify roof type from extracted planes.

    Args:
        planes: List of roof planes

    Returns:
        Roof type string: 'flat', 'gable', 'hip', 'complex'
    """
    if not planes:
        return "unknown"

    n_planes = len(planes)

    if n_planes == 1:
        # Check if plane is horizontal (flat roof)
        normal = planes[0].normal
        if abs(normal[2]) > 0.95:  # Nearly vertical normal
            return "flat"
        return "shed"

    if n_planes == 2:
        # Check for gable roof (two planes meeting at ridge)
        n1, n2 = planes[0].normal, planes[1].normal
        dot = abs(np.dot(n1, n2))
        if dot < 0.9:  # Planes at an angle
            return "gable"

    if n_planes == 4:
        return "hip"

    return "complex"


def refine_building_with_lpo(
    tile_id: str,
    footprint_bounds: Tuple[float, float, float, float],
    citygml_height: float,
    cfg: Config,
) -> Optional[float]:
    """Refine building height using LPO points.

    Args:
        tile_id: Tile identifier
        footprint_bounds: (xmin, ymin, xmax, ymax) of building
        citygml_height: Original height from CityGML
        cfg: Configuration

    Returns:
        Refined height or None if no LPO data
    """
    lpo_file = find_pointcloud_file(cfg.pointcloud.lpo_dir, tile_id)
    if not lpo_file:
        return None

    xmin, ymin, xmax, ymax = footprint_bounds

    # Load LPO points in building bounds
    lpo = read_pointcloud_file(lpo_file, bounds=(xmin, ymin, xmax, ymax))
    if len(lpo) < 10:
        return None

    # Get 95th percentile height
    max_z = float(np.percentile(lpo.z, 95))

    # Load DGM for ground elevation
    dgm_pattern = os.path.join(cfg.raw.dgm1_dir, f"*{tile_id}*.tif")
    dgm_files = glob.glob(dgm_pattern)
    if not dgm_files:
        return None

    dgm_ds = open_dataset(dgm_files[0], required=False)
    if dgm_ds is None:
        return None

    dgm = dgm_ds.GetRasterBand(1).ReadAsArray()
    gt = dgm_ds.GetGeoTransform()

    # Sample ground elevation at building center
    center_x = (xmin + xmax) / 2
    center_y = (ymin + ymax) / 2
    col = int((center_x - gt[0]) / gt[1])
    row = int((center_y - gt[3]) / gt[5])

    if 0 <= row < dgm.shape[0] and 0 <= col < dgm.shape[1]:
        ground_z = float(dgm[row, col])
        refined_height = max_z - ground_z
        if refined_height > 0:
            return refined_height

    return None


def sample_facade_color(
    tile_id: str,
    world_x: float,
    world_y: float,
    cfg: Config,
) -> Tuple[int, int, int]:
    """Sample facade color from orthophoto.

    Args:
        tile_id: Tile identifier
        world_x, world_y: World coordinates
        cfg: Configuration

    Returns:
        RGB tuple (0-255)
    """
    ortho_path = os.path.join(cfg.export.ortho_dir, f"{tile_id}.jpg")
    if not os.path.exists(ortho_path):
        return (180, 170, 160)  # Default neutral

    ds = open_dataset(ortho_path, required=False)
    if ds is None:
        return (180, 170, 160)

    gt = ds.GetGeoTransform()
    col = int((world_x - gt[0]) / gt[1])
    row = int((world_y - gt[3]) / gt[5])

    if 0 <= row < ds.RasterYSize and 0 <= col < ds.RasterXSize:
        try:
            r = ds.GetRasterBand(1).ReadAsArray(col, row, 1, 1)[0, 0]
            g = ds.GetRasterBand(2).ReadAsArray(col, row, 1, 1)[0, 0]
            b = ds.GetRasterBand(3).ReadAsArray(col, row, 1, 1)[0, 0]
            return (int(r), int(g), int(b))
        except Exception:
            pass

    return (180, 170, 160)


def process_building_with_bdom(
    tile_id: str,
    footprint_bounds: Tuple[float, float, float, float],
    cfg: Config,
) -> Optional[dict]:
    """Process a building using BDOM point cloud for roof geometry.

    Args:
        tile_id: Tile identifier
        footprint_bounds: (xmin, ymin, xmax, ymax) of building
        cfg: Configuration

    Returns:
        Dict with roof planes and metadata, or None if processing fails
    """
    eb_cfg = cfg.buildings_enhanced

    bdom_file = find_laz_file(cfg.pointcloud.bdom_dir, tile_id)
    if not bdom_file:
        return None

    xmin, ymin, xmax, ymax = footprint_bounds

    # Load BDOM points in building bounds (with small buffer)
    buffer = 0.5
    bdom = read_laz_file(bdom_file, bounds=(
        xmin - buffer, ymin - buffer,
        xmax + buffer, ymax + buffer
    ))

    if len(bdom) < eb_cfg.min_plane_inliers:
        return None

    # Load DGM for ground reference
    dgm_pattern = os.path.join(cfg.raw.dgm1_dir, f"*{tile_id}*.tif")
    dgm_files = glob.glob(dgm_pattern)
    if not dgm_files:
        return None

    dgm_ds = open_dataset(dgm_files[0], required=False)
    if dgm_ds is None:
        return None

    dgm = dgm_ds.GetRasterBand(1).ReadAsArray()
    gt = dgm_ds.GetGeoTransform()

    # Get ground elevation at building center
    center_x = (xmin + xmax) / 2
    center_y = (ymin + ymax) / 2
    col = int((center_x - gt[0]) / gt[1])
    row = int((center_y - gt[3]) / gt[5])

    if 0 <= row < dgm.shape[0] and 0 <= col < dgm.shape[1]:
        ground_z = float(dgm[row, col])
    else:
        ground_z = float(np.min(bdom.z))

    # Filter to points above ground (roof points)
    roof_threshold = ground_z + 2.0  # At least 2m above ground
    roof_mask = bdom.z > roof_threshold
    roof_points = PointCloud(
        x=bdom.x[roof_mask],
        y=bdom.y[roof_mask],
        z=bdom.z[roof_mask],
        r=bdom.r[roof_mask] if bdom.r is not None else None,
        g=bdom.g[roof_mask] if bdom.g is not None else None,
        b=bdom.b[roof_mask] if bdom.b is not None else None,
    )

    if len(roof_points) < eb_cfg.min_plane_inliers:
        return None

    # Extract roof planes
    planes = extract_roof_planes(
        roof_points,
        max_planes=eb_cfg.max_roof_planes,
        ransac_threshold=eb_cfg.ransac_threshold_m,
        ransac_iterations=eb_cfg.ransac_iterations,
        min_inliers=eb_cfg.min_plane_inliers,
    )

    if not planes:
        return None

    # Get roof RGB if available
    roof_rgb = None
    if roof_points.has_rgb:
        r_avg = int(np.mean(roof_points.r))
        g_avg = int(np.mean(roof_points.g))
        b_avg = int(np.mean(roof_points.b))
        roof_rgb = (r_avg, g_avg, b_avg)

    # Compute building height
    max_z = float(np.max(roof_points.z))
    height = max_z - ground_z

    return {
        "planes": planes,
        "roof_type": classify_roof_type(planes),
        "height": height,
        "ground_z": ground_z,
        "roof_rgb": roof_rgb,
        "source": "bdom",
    }


def export_buildings_enhanced(cfg: Config) -> int:
    """Export enhanced buildings with point cloud fusion.

    For tiles with BDOM20RGBI:
    - Extract roof planes using RANSAC
    - Use RGB from point cloud for facade colors

    For tiles with only LPO:
    - Refine building heights from LPO
    - Use DOP20 for facade colors

    For tiles without point cloud:
    - Fall back to standard CityGML pipeline

    Args:
        cfg: Configuration

    Returns:
        Exit code (0 = success)
    """
    eb_cfg = cfg.buildings_enhanced

    # Check for manifest
    if not os.path.exists(cfg.export.manifest_path):
        print(f"[buildings_enhanced] ERROR: Manifest not found: {cfg.export.manifest_path}")
        print("[buildings_enhanced] Run heightmap export first.")
        return 1

    # Load tiles
    tiles = []
    with open(cfg.export.manifest_path, "r") as f:
        reader = csv.DictReader(f)
        for row in reader:
            tiles.append({
                "tile_id": row["tile_id"],
                "xmin": float(row["xmin"]),
                "ymin": float(row["ymin"]),
                "xmax": float(row["xmax"]),
                "ymax": float(row["ymax"]),
            })

    if not tiles:
        print("[buildings_enhanced] No tiles in manifest.")
        return 1

    os.makedirs(eb_cfg.out_dir, exist_ok=True)

    stats = {"bdom": 0, "lpo": 0, "citygml": 0, "failed": 0}

    for tile in tiles:
        tile_id = tile["tile_id"]
        print(f"[buildings_enhanced] Processing {tile_id}...")

        # Check data availability
        has_bdom = has_bdom_data(tile_id, cfg.pointcloud.bdom_dir)
        has_lpo = has_lpo_data(tile_id, cfg.pointcloud.lpo_dir)

        if has_bdom and eb_cfg.use_bdom_roof:
            print(f"[buildings_enhanced] {tile_id}: Using BDOM for roof extraction")
            stats["bdom"] += 1
            # BDOM processing would create enhanced GLB here
            # For now, we log that BDOM is available

        elif has_lpo and eb_cfg.use_lpo_refinement:
            print(f"[buildings_enhanced] {tile_id}: Using LPO for height refinement")
            stats["lpo"] += 1
            # LPO refinement would adjust CityGML heights here

        elif eb_cfg.fallback_to_citygml:
            print(f"[buildings_enhanced] {tile_id}: Falling back to CityGML")
            stats["citygml"] += 1
            # Standard CityGML export

        else:
            print(f"[buildings_enhanced] {tile_id}: No data available, skipping")
            stats["failed"] += 1

    print(f"[buildings_enhanced] DONE. "
          f"BDOM={stats['bdom']}, LPO={stats['lpo']}, "
          f"CityGML={stats['citygml']}, Failed={stats['failed']}")

    # Note: Full GLB export would integrate with existing buildings.py
    # This module provides the data fusion logic; actual mesh export
    # requires integration with trimesh and GLB generation from buildings.py

    if stats["bdom"] > 0 or stats["lpo"] > 0:
        print("[buildings_enhanced] Note: Full GLB export requires integration with buildings.py")
        print("[buildings_enhanced] Use --export buildings for standard CityGML export")

    return 0
