424 lines
18 KiB
Python
424 lines
18 KiB
Python
"""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")
|