"""Prueba técnica: reconstrucción automática de una pieza 3D a partir de un plano DXF
con tres vistas (alzado, planta y perfil).

Método (descomposición en celdas):
  1. Leer las líneas visibles y ocultas; ignorar marco, cajetín, textos y cotas.
  2. Buscar la escala en los textos del plano ("ESCALA 1:2").
  3. Agrupar las líneas en vistas e identificar alzado, planta y perfil por su alineación.
  4. Para cada convención de proyección posible (europea/americana en planta y perfil):
       a. dividir el volumen en celdas con todos los planos que sugieren las vistas
          (incluidos los planos inclinados de las líneas diagonales);
       b. descartar celdas que se salen de la silueta de alguna vista;
       c. probar combinaciones de celdas cuya silueta coincide con las tres vistas;
       d. proyectar las aristas del sólido (visibles y ocultas) y compararlas con el plano.
  5. Quedarse con el sólido que reproduce el plano y exportarlo a STEP y STL.

Uso:  python reconstruir.py CUBE01.dxf
"""
import itertools
import math
import re
import sys
from pathlib import Path

import cadquery as cq
import ezdxf
import numpy as np
from OCP.BRepClass3d import BRepClass3d_SolidClassifier
from OCP.BRepIntCurveSurface import BRepIntCurveSurface_Inter
from OCP.gp import gp_Dir, gp_Lin, gp_Pnt
from OCP.TopAbs import TopAbs_IN
from shapely.geometry import LineString, MultiPoint
from shapely.ops import polygonize, unary_union

TOL = 1.0          # tolerancia en unidades de dibujo (las vistas no están perfectamente alineadas)
BIG = 1e4
MAX_CELDAS = 18    # límite de la búsqueda exhaustiva (2^18 combinaciones)


# ---------------------------------------------------------------- 1. lectura
def leer_dxf(ruta):
    doc = ezdxf.readfile(ruta)
    msp = doc.modelspace()
    avisos = []

    def es_oculta(e):
        lt = e.dxf.get("linetype", "ByLayer")
        if lt.lower() in ("bylayer", "byblock"):
            lt = doc.layers.get(e.dxf.layer).dxf.linetype
        nombre = (lt + " " + e.dxf.layer).upper()
        return any(k in nombre for k in ("HIDDEN", "DASH", "OCULT", "ISO02"))

    segs = []
    for e in msp:
        if e.dxftype() == "LINE":
            partes = [(e.dxf.start, e.dxf.end)]
        elif e.dxftype() == "LWPOLYLINE":
            pts = [p[:2] for p in e.get_points()]
            if e.closed:
                pts.append(pts[0])
            partes = list(zip(pts, pts[1:]))
        else:
            continue  # marco y cajetín (INSERT), textos, cotas, etc.
        for a, b in partes:
            segs.append((float(a[0]), float(a[1]), float(b[0]), float(b[1]), es_oculta(e)))

    textos = [t.dxf.text for t in msp.query("TEXT")] + [t.text for t in msp.query("MTEXT")]
    for ins in msp.query("INSERT"):
        textos += [a.dxf.text for a in ins.attribs]
    escala = None
    for t in textos:
        m = re.search(r"(?:ESCALA|SCALE)\s*:?\s*(\d+(?:[.,]\d+)?)\s*[:/]\s*(\d+(?:[.,]\d+)?)", t, re.I)
        if m:
            a, b = (float(x.replace(",", ".")) for x in m.groups())
            escala = (a, b)
            break
    if escala is None:
        avisos.append("No se ha encontrado la escala en el plano: se supone 1:1.")
        escala = (1.0, 1.0)

    n_cotas = len(msp.query("DIMENSION"))
    if n_cotas == 0:
        avisos.append("El plano no tiene cotas: las medidas se deducen de la geometría y la escala.")
    unidades = {0: "sin unidades", 1: "pulgadas", 4: "mm", 5: "cm", 6: "m"}.get(doc.header.get("$INSUNITS", 0), "otras")
    return segs, escala, n_cotas, unidades, avisos


# ---------------------------------------------------------------- 2. vistas
def agrupar_vistas(segs):
    n = len(segs)
    padre = list(range(n))

    def raiz(i):
        while padre[i] != i:
            padre[i] = padre[padre[i]]
            i = padre[i]
        return i

    cajas = [(min(s[0], s[2]) - TOL, min(s[1], s[3]) - TOL, max(s[0], s[2]) + TOL, max(s[1], s[3]) + TOL) for s in segs]
    for i in range(n):
        for j in range(i + 1, n):
            a, b = cajas[i], cajas[j]
            if a[0] <= b[2] and b[0] <= a[2] and a[1] <= b[3] and b[1] <= a[3]:
                padre[raiz(i)] = raiz(j)
    grupos = {}
    for i in range(n):
        grupos.setdefault(raiz(i), []).append(segs[i])
    vistas = []
    for g in grupos.values():
        xs = [c for s in g for c in (s[0], s[2])]
        ys = [c for s in g for c in (s[1], s[3])]
        vistas.append({"segs": g, "u0": min(xs), "u1": max(xs), "v0": min(ys), "v1": max(ys)})
    return vistas


def identificar_vistas(vistas):
    def mismo_ancho(a, b):
        return abs(a["u0"] - b["u0"]) < 2 * TOL and abs(a["u1"] - b["u1"]) < 2 * TOL

    def mismo_alto(a, b):
        return abs(a["v0"] - b["v0"]) < 2 * TOL and abs(a["v1"] - b["v1"]) < 2 * TOL

    for alz in vistas:
        plantas = [v for v in vistas if v is not alz and mismo_ancho(alz, v)]
        perfiles = [v for v in vistas if v is not alz and mismo_alto(alz, v)]
        if plantas and perfiles:
            pl, pf = plantas[0], perfiles[0]
            return alz, pl, pf, ("debajo" if pl["v1"] < alz["v0"] else "encima"), ("derecha" if pf["u0"] > alz["u1"] else "izquierda")
    raise SystemExit("No se han podido identificar alzado, planta y perfil alineados.")


def a_coordenadas(alz, pl, pf, pos_pl, pos_pf, conv_pl, conv_pf):
    """Pasa cada vista a coordenadas de la pieza: alzado (X,Z), planta (X,Y), perfil (Y,Z).
    Y = 0 es la cara delantera. Devuelve también desde dónde mira el observador de cada vista."""
    D = pl["v1"] - pl["v0"]

    def map_pl(u, v):
        d = (pl["v1"] - v) if pos_pl == "debajo" else (v - pl["v0"])   # distancia al borde junto al alzado
        return (u - pl["u0"], D - d if conv_pl == "E" else d)

    def map_pf(u, v):
        d = (u - pf["u0"]) if pos_pf == "derecha" else (pf["u1"] - u)
        return (D - d if conv_pf == "E" else d, v - pf["v0"])

    def map_alz(u, v):
        return (u - alz["u0"], v - alz["v0"])

    out = {}
    for nombre, vista, f in (("alzado", alz, map_alz), ("planta", pl, map_pl), ("perfil", pf, map_pf)):
        out[nombre] = [(*f(s[0], s[1]), *f(s[2], s[3]), s[4]) for s in vista["segs"]]

    # Observador: en el sistema europeo la vista de debajo es la superior; en el americano, al revés.
    desde_arriba = (pos_pl == "debajo") == (conv_pl == "E")
    desde_izquierda = (pos_pf == "derecha") == (conv_pf == "E")
    observador = {
        "alzado": (0, -1, 0),
        "planta": (0, 0, 1) if desde_arriba else (0, 0, -1),
        "perfil": (-1, 0, 0) if desde_izquierda else (1, 0, 0),
    }
    return out, observador, D


PROY = {"alzado": (0, 2), "planta": (0, 1), "perfil": (1, 2)}   # ejes 3D que se ven en cada vista


# ---------------------------------------------------------------- 3. celdas
def agrupar_valores(vals):
    vals = sorted(vals)
    out = []
    for v in vals:
        if not out or v - out[-1] > TOL:
            out.append(v)
    return out


def region(segs):
    lineas = [LineString([(s[0], s[1]), (s[2], s[3])]) for s in segs if not s[4]]
    return unary_union(list(polygonize(unary_union(lineas))))


def semiespacio(p, n):
    n = cq.Vector(*n).normalized()
    aux = cq.Vector(1, 0, 0) if abs(n.x) < 0.9 else cq.Vector(0, 1, 0)
    plano = cq.Plane(origin=cq.Vector(*p), xDir=aux.cross(n).normalized(), normal=n)
    return cq.Workplane(plano).rect(BIG, BIG).extrude(BIG).val()


def crear_celdas(coords):
    rejilla = {0: [], 1: [], 2: []}
    inclinados = []
    for nombre, segs in coords.items():
        a, b = PROY[nombre]
        for x1, y1, x2, y2, _ in segs:
            rejilla[a] += [x1, x2]
            rejilla[b] += [y1, y2]
            if abs(x1 - x2) > TOL and abs(y1 - y2) > TOL:   # línea diagonal -> plano inclinado
                n = [0.0, 0.0, 0.0]
                n[a], n[b] = -(y2 - y1), (x2 - x1)
                p = [0.0, 0.0, 0.0]
                p[a], p[b] = x1, y1
                inclinados.append((tuple(p), tuple(n)))
    ejes = [agrupar_valores(rejilla[i]) for i in range(3)]

    celdas = []
    for i, j, k in itertools.product(range(len(ejes[0]) - 1), range(len(ejes[1]) - 1), range(len(ejes[2]) - 1)):
        x0, x1 = ejes[0][i], ejes[0][i + 1]
        y0, y1 = ejes[1][j], ejes[1][j + 1]
        z0, z1 = ejes[2][k], ejes[2][k + 1]
        trozos = [cq.Solid.makeBox(x1 - x0, y1 - y0, z1 - z0, pnt=cq.Vector(x0, y0, z0))]
        esquinas = np.array(list(itertools.product((x0, x1), (y0, y1), (z0, z1))))
        for p, n in inclinados:
            dist = (esquinas - np.array(p)) @ np.array(n)
            if dist.min() < -1e-6 and dist.max() > 1e-6:      # el plano atraviesa la caja
                hs = semiespacio(p, n)
                nuevos = []
                for t in trozos:
                    for parte in (t.intersect(hs), t.cut(hs)):
                        nuevos += [s for s in parte.Solids() if s.Volume() > 1e-6]
                trozos = nuevos
        celdas += trozos
    return celdas, ejes, inclinados


def proyeccion(solido, nombre):
    a, b = PROY[nombre]
    pts = [(v.toTuple()[a], v.toTuple()[b]) for v in solido.Vertices()]
    return MultiPoint(pts).convex_hull


# ---------------------------------------------------------------- 4. comparación
def oculto(forma, p, d):
    inter = BRepIntCurveSurface_Inter()
    inter.Init(forma.wrapped, gp_Lin(gp_Pnt(*p), gp_Dir(*d)), 1e-7)
    ws = []
    while inter.More():
        if inter.W() > 1e-5:
            ws.append(inter.W())
        inter.Next()
    params = [0.0] + sorted(set(round(w, 6) for w in ws))
    for t0, t1 in zip(params, params[1:]):
        q = np.array(p) + np.array(d) * (t0 + t1) / 2
        if BRepClass3d_SolidClassifier(forma.wrapped, gp_Pnt(*q), 1e-6).State() == TopAbs_IN:
            return True
    return False


def aristas_proyectadas(forma, nombre, obs):
    a, b = PROY[nombre]
    vis, ocu = [], []
    for e in forma.Edges():
        p0, p1 = np.array(e.startPoint().toTuple()), np.array(e.endPoint().toTuple())
        if math.dist((p0[a], p0[b]), (p1[a], p1[b])) < TOL / 2:
            continue                                   # arista perpendicular a la vista
        n = 5
        for i in range(n):
            q0, q1 = p0 + (p1 - p0) * i / n, p0 + (p1 - p0) * (i + 1) / n
            seg = (q0[a], q0[b], q1[a], q1[b])
            (ocu if oculto(forma, (q0 + q1) / 2, obs) else vis).append(seg)
    return vis, ocu


def muestras(segs, paso=TOL / 2):
    pts = []
    for x1, y1, x2, y2 in segs:
        n = max(1, int(math.dist((x1, y1), (x2, y2)) / paso))
        for i in range(n):
            t = (i + 0.5) / n
            pts.append((x1 + (x2 - x1) * t, y1 + (y2 - y1) * t))
    return np.array(pts).reshape(-1, 2)


def distancia(pts, segs):
    if len(pts) == 0:
        return np.zeros(0)
    if not segs:
        return np.full(len(pts), np.inf)
    s = np.array(segs, dtype=float)
    a, b = s[:, :2], s[:, 2:]
    ab = b - a
    t = np.clip(((pts[:, None, :] - a) * ab).sum(-1) / np.maximum((ab * ab).sum(-1), 1e-12), 0, 1)
    proj = a + t[..., None] * ab
    return np.linalg.norm(pts[:, None, :] - proj, axis=-1).min(axis=1)


def errores_vista(forma, nombre, obs, dibujo):
    cvis, cocu = aristas_proyectadas(forma, nombre, obs)
    dvis = [s[:4] for s in dibujo if not s[4]]
    docu = [s[:4] for s in dibujo if s[4]]
    err = 0
    err += int((distancia(muestras(cvis), dvis) > TOL).sum())                 # aristas que el plano no tiene
    m = muestras(cocu)
    m = m[distancia(m, cvis) > TOL]                                           # lo visible tapa lo oculto
    err += int((distancia(m, docu) > TOL).sum())
    err += int((distancia(muestras(dvis), cvis) > TOL).sum())                 # líneas del plano que faltan
    err += int((distancia(muestras(docu), cocu) > TOL).sum())
    return err


# ---------------------------------------------------------------- 5. búsqueda
def reconstruir(coords, observador):
    regiones = {n: region(s) for n, s in coords.items()}
    celdas, ejes, inclinados = crear_celdas(coords)
    proy = [{n: proyeccion(c, n) for n in coords} for c in celdas]

    # Una celda que se sale de la silueta de alguna vista no puede formar parte de la pieza.
    validas = [i for i in range(len(celdas))
               if all(proy[i][n].difference(regiones[n].buffer(TOL / 2)).area < 1e-3 for n in coords)]
    if len(validas) > MAX_CELDAS:
        raise SystemExit(f"Demasiadas celdas ({len(validas)}) para la búsqueda exhaustiva de la prueba.")

    candidatos = []
    for r in range(1, len(validas) + 1):
        for combo in itertools.combinations(validas, r):
            ok = all(unary_union([proy[i][n] for i in combo]).symmetric_difference(regiones[n]).area
                     < 0.005 * regiones[n].area for n in coords)
            if ok:
                candidatos.append(combo)

    resultados = []
    for combo in candidatos:
        forma = celdas[combo[0]]
        if len(combo) > 1:
            forma = forma.fuse(*[celdas[i] for i in combo[1:]])
        forma = forma.clean()
        if len(forma.Solids()) != 1:
            continue                                   # piezas sueltas: no es una pieza única
        forma = forma.Solids()[0]
        err = sum(errores_vista(forma, n, observador[n], coords[n]) for n in coords)
        resultados.append((err, combo, forma))
        if err == 0 and len([r for r in resultados if r[0] == 0]) > 3:
            break
    resultados.sort(key=lambda r: r[0])
    return resultados, len(celdas), len(validas), len(candidatos), ejes, inclinados


# ---------------------------------------------------------------- salida
def vista_previa(forma, ruta_png):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection

    verts, tris = forma.tessellate(0.1)
    v = np.array([p.toTuple() for p in verts])
    caras = v[np.array(tris)]
    normales = np.cross(caras[:, 1] - caras[:, 0], caras[:, 2] - caras[:, 0])
    normales /= np.linalg.norm(normales, axis=1, keepdims=True) + 1e-12
    luz = np.clip(normales @ np.array([0.4, -0.6, 0.7]), 0.15, 1)
    fig = plt.figure(figsize=(6, 6))
    ax = fig.add_subplot(projection="3d")
    ax.add_collection3d(Poly3DCollection(caras, facecolors=plt.cm.Blues(0.35 + 0.55 * luz), edgecolor="none"))
    for e in forma.Edges():
        a, b = e.startPoint().toTuple(), e.endPoint().toTuple()
        ax.plot(*zip(a, b), color="k", lw=0.8)
    lo, hi = v.min(0), v.max(0)
    ax.set_xlim(lo[0], hi[0]); ax.set_ylim(lo[1], hi[1]); ax.set_zlim(lo[2], hi[2])
    ax.set_box_aspect(hi - lo)
    ax.set_xlabel("X (mm)"); ax.set_ylabel("Y (mm, 0 = frente)"); ax.set_zlabel("Z (mm)")
    ax.view_init(elev=25, azim=-60)
    plt.tight_layout()
    plt.savefig(ruta_png, dpi=110)


def main(ruta):
    ruta = Path(ruta)
    segs, (ea, eb), n_cotas, unidades, avisos = leer_dxf(ruta)
    factor = eb / ea
    vistas = agrupar_vistas(segs)
    alz, pl, pf, pos_pl, pos_pf = identificar_vistas(vistas)

    inf = [f"# Informe de reconstrucción: {ruta.name}", "",
           f"- Líneas leídas: {len(segs)} ({sum(s[4] for s in segs)} ocultas); unidades del archivo: {unidades}",
           f"- Escala: {ea:g}:{eb:g} -> factor x{factor:g}; cotas en el plano: {n_cotas}",
           f"- Vistas encontradas: {len(vistas)}; planta {pos_pl} del alzado, perfil a la {pos_pf}",
           f"- Desalineación planta/alzado: {abs(pl['u0'] - alz['u0']):.2f}; perfil/alzado: {abs(pf['v0'] - alz['v0']):.2f} (unidades de dibujo)"]
    D_pf = pf["u1"] - pf["u0"]
    D_pl = pl["v1"] - pl["v0"]
    if abs(D_pf - D_pl) > 2 * TOL:
        avisos.append(f"La profundidad no coincide entre planta ({D_pl:.1f}) y perfil ({D_pf:.1f}).")

    nombres = {"E": "europeo (ISO-E)", "A": "americano (ISO-A)"}
    mejor = None
    inf += ["", "## Hipótesis de proyección", "", "| Planta | Perfil | Celdas (válidas) | Candidatos por silueta | Error mínimo |", "|---|---|---|---|---|"]
    for conv_pl, conv_pf in itertools.product("EA", repeat=2):
        coords, obs, D = a_coordenadas(alz, pl, pf, pos_pl, pos_pf, conv_pl, conv_pf)
        res, n_c, n_v, n_cand, ejes, incl = reconstruir(coords, obs)
        e_min = res[0][0] if res else "-"
        inf.append(f"| {nombres[conv_pl]} | {nombres[conv_pf]} | {n_c} ({n_v}) | {n_cand} | {e_min} |")
        print(f"planta {conv_pl} / perfil {conv_pf}: celdas {n_c}, validas {n_v}, candidatos {n_cand}, error {e_min}")
        if res and (mejor is None or res[0][0] < mejor[0][0][0]):
            mejor = (res, conv_pl, conv_pf)

    if mejor is None:
        raise SystemExit("No se ha encontrado ninguna pieza compatible con las vistas.")
    res, conv_pl, conv_pf = mejor
    err, combo, forma = res[0]
    exactas = [r for r in res if r[0] == 0]
    if conv_pl != conv_pf:
        avisos.append(f"Sistemas de proyección mezclados: la planta sigue el sistema {nombres[conv_pl]} "
                      f"y el perfil el {nombres[conv_pf]}.")
    if err > 0:
        avisos.append(f"Ninguna pieza reproduce exactamente el plano (error mínimo {err}): revisar manualmente.")
    if len(exactas) > 1:
        avisos.append(f"Hay {len(exactas)} piezas distintas que encajan con las vistas: el plano es ambiguo.")

    forma = forma.scale(factor)
    bb = forma.BoundingBox()
    inf += ["", "## Resultado", "",
            f"- Convención elegida: planta {nombres[conv_pl]}, perfil {nombres[conv_pf]}",
            f"- Celdas usadas: {len(combo)}; error de reproducción del plano: {err}",
            f"- Tamaño real: {bb.xlen:.1f} x {bb.ylen:.1f} x {bb.zlen:.1f} mm; volumen {forma.Volume():.0f} mm³",
            f"- Caras: {len(forma.Faces())}, aristas: {len(forma.Edges())}",
            "", "## Avisos", ""] + [f"- {a}" for a in avisos]

    base = ruta.with_suffix("")
    cq.exporters.export(cq.Workplane().add(forma), str(base) + ".step")
    cq.exporters.export(cq.Workplane().add(forma), str(base) + ".stl")
    vista_previa(forma, str(base) + "_3d.png")
    Path(str(base) + "_informe.md").write_text("\n".join(inf) + "\n", encoding="utf-8")
    print("\n".join(inf))


if __name__ == "__main__":
    main(sys.argv[1] if len(sys.argv) > 1 else "CUBE01.dxf")