-- coding: utf-8 --

"""
FRUGAL -> STRUCTURAL 3D CAD (V23.2, evidence-augmented arrangement)

Базовая геометрия V23/V23.1 сохранена:

• физические плоскости не меняются
• structural intersection lines не меняются
• core-ячейки строятся как раньше

Добавлено (V23.2):

• Recovery pre-pass: точки sid==0 классифицируются и раскладываются
  по поверхностям как ДОПОЛНИТЕЛЬНЫЙ evidence (не заменяют core).
• Arrangement использует dual evidence (core + recovery):
  • core_ok → evidence_source="core"
  • core_ok & rec → "mixed"
  • rec_ok & core_min → "recovered"
• Recovered-ячейки экспортируются отдельным слоем CAD_EDGES_RECOVERED
  в DXF и отдельной колонкой в cad_faces_diag.csv.
• Ambiguous recovery-точки (top-kNN соседи принадлежат ≥2 surfaces)
  не используются как evidence и сохраняются отдельно для диагностики.
  """

import os
import csv
import math
import logging
from dataclasses import dataclass
from collections import defaultdict

import numpy as np
from scipy.spatial import cKDTree
from scipy import ndimage

import laspy

try:
from skimage.measure import find_contours, approximate_polygon
HAS_SKIMAGE = True
except ImportError:
HAS_SKIMAGE = False

try:
import ezdxf
HAS_DXF = True
except ImportError:
HAS_DXF = False

try:
from plyfile import PlyData, PlyElement
HAS_PLY = True
except ImportError:
HAS_PLY = False

try:
import cadquery as cq
HAS_CQ = True
except ImportError:
HAS_CQ = False

try:
from shapely.geometry import Polygon, MultiPoint, LineString, Point
from shapely.ops import unary_union, polygonize
HAS_SHAPELY = True
except ImportError:
HAS_SHAPELY = False

logging.basicConfig(
level=logging.INFO,
format="%(asctime)s | F2CAD | %(levelname)s | %(message)s")
log = logging.getLogger("F2CAD")

============================================================

CONFIG

============================================================

@dataclass
class Config:
frugal_las: str = ""
edge_las: str = ""
rooms_npz: str = ""
out_dir: str = ""

text
run_roughness_b: bool = False
run_A_completion: bool = True
run_C_merge: bool = False
b_hard_cap: int = 5000

surface_min_room_frac: float = .20
region_min_pts: int = 100
wall_nz_max: float = .25
horizontal_nz_min: float = .88
plane_rms_max: float = .045

surface_merge_angle: float = 5.0
surface_merge_plane_dist: float = .08
surface_merge_spatial_gap: float = .80
surface_merge_center_radius: float = 4.0
surface_merge_rms_max: float = .08
surface_merge_passes: int = 2

adjacency_dist: float = .25
adjacency_min_pts: int = 4
adjacency_sample: int = 4000

parallel_angle_deg: float = 10.
line_surface_dist: float = .16
line_surface_min: int = 4

trim_q: float = .015
min_edge_len: float = .12
edge_extent_extension: float = .12

edge_line_dist: float = .15
edge_search_pad: float = .50
edge_min_pts: int = 4

c_cell: float = .08
c_min_cell_pts: int = 1
c_simplify: float = .10
c_min_len: float = .30
c_sample: int = 15000
c_max_raster_cells: int = 4_000_000
c_min_island_cells: int = 4
c_vertical_only: bool = True

a_completion_min_len: float = .30

fuse_angle_deg: float = 4.0
fuse_lateral_dist: float = .05
fuse_gap: float = .25
fuse_min_len: float = .15

topo_tol: float = .03

triple_plane_tol: float = .20
endpoint_vertex_tol: float = .15
triple_endpoint_extension: float = .30

preview_step: float = .025

run_cad_stage2: bool = True
cad_export_ply: bool = True
cad_export_dae: bool = True
cad_export_dxf: bool = True
cad_dxf_3dface: bool = True

cad_tess_tolerance: float = 0.005
cad_tess_angular: float = 0.20

cad_use_only_measured_A: bool = True
cad_min_face_area: float = 0.02
cad_collinear_tol: float = 0.03
cad_min_loop_vertices: int = 3

cad_use_c_measured_loops: bool = True
cad_extract_c_for_horizontals: bool = True
cad_c_cell: float = .05
cad_c_simplify: float = .15
cad_c_min_island_cells: int = 4
cad_c_sample: int = 20000
cad_c_max_raster_cells: int = 4_000_000
cad_c_fill_holes: bool = False
cad_c_inner_min_area: float = 0.30
cad_c_close_iter: int = 2
cad_c_open_iter: int = 1

cad_reg_enabled: bool = True
cad_reg_min_run_len: float = 0.15
cad_reg_run_angle_deg: float = 15.0
cad_reg_snap_A_dist: float = 0.15
cad_reg_snap_A_angle_deg: float = 8.0
cad_reg_ortho_angle_deg: float = 10.0
cad_reg_merge_vertex_tol: float = 0.02

cad_wall_basis: bool = True
cad_min_polygon_area: float = 0.05
cad_fallback_concave: bool = True
cad_concave_ratio: float = 0.30
cad_consistency_sample: int = 5000

structural_enabled: bool = True
structural_band: float = 0.20
structural_extend_m: float = 0.20
structural_max_extend: float = 0.50
structural_min_len: float = 0.10
vertex_merge_tol: float = 0.03
cad_plane_res_tol: float = 0.01

structural_external_tol: float = 0.05
structural_external_min_wrong: int = 3
structural_coverage_min: float = 0.90
structural_obs_band: float = 0.15

# V23 arrangement
structural_arrange_extend: float = 0.60
structural_cell_min_pts: int = 5
structural_cell_min_frac: float = 0.02
structural_cell_merge_tol: float = 0.05
horizontal_level_tol: float = 0.10

# V23.1 evidence & provenance
provenance_line_tol: float = 0.05
provenance_parallel_cos: float = 0.98
shared_edge_tol: float = 0.05
shared_edge_min_len: float = 0.05
unsupported_min_pts: int = 15
unsupported_min_density: float = 30.0
evidence_sample: int = 2000

# V23.2 recovery pre-pass
recovery_enabled: bool = True
recovery_knn_assigned: int = 8
recovery_knn_unassigned: int = 12
recovery_dist_recoverable: float = 0.02
recovery_dist_extendable: float = 0.05
recovery_angle_recoverable_deg: float = 15.0
recovery_angle_extendable_deg: float = 25.0
recovery_uv_extend_m: float = 0.10
recovery_exclude_ambiguous: bool = True
recovery_chunk: int = 200_000

# пороги для ячейки, поддержанной только recovery
structural_cell_min_recovered_pts: int = 20
structural_cell_min_recovered_frac: float = 0.08
structural_cell_min_core_pts_rec: int = 3

============================================================

GEOMETRY

============================================================

def unit(x):
return x / (np.linalg.norm(x) + 1e-12)

def canonical_plane(n, d):
n = unit(n)
k = np.argmax(abs(n))
return (-n, -d) if n[k] < 0 else (n, float(d))

def angle_planes(a, b):
c = np.clip(abs(np.dot(unit(a), unit(b))), -1, 1)
return math.degrees(math.acos(c))

def angle_dirs(a, b):
c = np.clip(abs(np.dot(unit(a), unit(b))), -1, 1)
return math.degrees(math.acos(c))

def fit_plane(P):
c = P.mean(0)
Q = P - c
C = Q.T @ Q / max(len(P) - 1, 1)
w, V = np.linalg.eigh(C)
n = unit(V[:, np.argmin(w)])
d = -float(n @ c)
n, d = canonical_plane(n, d)
r = np.abs(P @ n + d)
return n, d, c, float(np.sqrt(np.mean(r * r)))

def plane_type(n, cfg):
z = abs(float(n[2]))
if z >= cfg.horizontal_nz_min:
return "horizontal"
if z <= cfg.wall_nz_max:
return "vertical"
return "inclined"

def plane_intersection(A, B):
n1, n2 = A["n"], B["n"]
v = np.cross(n1, n2)
den = float(v @ v)
if den < 1e-10:
return None
c1, c2 = -A["d"], -B["d"]
p = np.cross(c1 * n2 - c2 * n1, v) / den
return p, unit(v)

def line_dist(P, p, d):
return np.linalg.norm(np.cross(P - p, unit(d)), axis=1)

def triple_intersection(A, B, C):
M = np.vstack([A["n"], B["n"], C["n"]])
if abs(np.linalg.det(M)) < 1e-7:
return None
b = -np.array([A["d"], B["d"], C["d"]])
return np.linalg.solve(M, b)

def plane_basis(n):
n = unit(n)
a = np.array([0., 0., 1.]) if abs(n[2]) < .9 else np.array([1., 0., 0.])
u = unit(np.cross(n, a))
v = unit(np.cross(n, u))
return u, v

def wall_basis(n):
n = unit(n)
z = np.array([0., 0., 1.])
along = np.cross(z, n)
if np.linalg.norm(along) < 1e-8:
return plane_basis(n)
u = unit(along)
v = unit(np.cross(n, u))
if v[2] < 0:
v = -v
u = -u
return u, v

def surface_basis(s):
if s.get("type") == "vertical":
return wall_basis(s["n"])
return plane_basis(s["n"])

def uv_of_points(P, s, u=None, v=None):
if u is None or v is None:
u, v = surface_basis(s)
q = np.asarray(P, float) - s["c"]
if q.ndim == 1:
return np.array([q @ u, q @ v])
return np.c_[q @ u, q @ v]

def points_of_uv(uv, s, u=None, v=None):
if u is None or v is None:
u, v = surface_basis(s)
uv = np.asarray(uv, float)
if uv.ndim == 1:
return s["c"] + uv[0] * u + uv[1] * v
return s["c"] + np.outer(uv[:, 0], u) + np.outer(uv[:, 1], v)

def project_point_to_plane(x, s):
return x - (np.dot(x, s["n"]) + s["d"]) * s["n"]

============================================================

UNION-FIND

============================================================

class UF:
def init(self, n):
self.p = np.arange(n)

text
def f(self, x):
    while self.p[x] != x:
        self.p[x] = self.p[self.p[x]]
        x = self.p[x]
    return x

def u(self, a, b):
    a, b = self.f(a), self.f(b)
    if a != b:
        self.p[b] = a

============================================================

IO

============================================================

def load_frugal(path):
x = laspy.read(path)
names = set(x.point_format.dimension_names)
if "segment_id" not in names:
raise RuntimeError("FRUGAL LAS has no segment_id")
P = np.c_[x.x, x.y, x.z].astype(float)
sid = np.asarray(x.segment_id, dtype=np.int64)
if "room_id" in names:
rid = np.asarray(x.room_id, dtype=np.int32)
has_room = True
else:
rid = np.zeros(len(P), np.int32)
has_room = False
frac = 100.0 * float(np.mean(rid > 0)) if len(P) else 0.0
log.info("room_id dimension: %s, pts with room=%.1f%%", has_room, frac)
return P, sid, rid

def load_edges(path):
if not path or not os.path.exists(path):
log.warning("roughness edge file not found")
return np.empty((0, 3))
x = laspy.read(path)
return np.c_[x.x, x.y, x.z].astype(float)

def load_rooms_npz(path):
if not path or not os.path.exists(path):
log.warning("rooms NPZ not found: %s", path)
return None
data = np.load(path)
out = {
"room_grid": data["room_grid"] if "room_grid" in data else None,
"grid_origin": data["grid_origin"] if "grid_origin" in data else None,
"grid_cell": float(data["grid_cell"]) if "grid_cell" in data else None,
"room_id": data["room_id"] if "room_id" in data else None,
"centroid": data["centroid"] if "centroid" in data else None,
"bbox_min": data["bbox_min"] if "bbox_min" in data else None,
"bbox_max": data["bbox_max"] if "bbox_max" in data else None,
"area": data["area"] if "area" in data else None,
}
rr_reg = data["region_room_region"] if "region_room_region" in data else None
rr_room = data["region_room_room"] if "region_room_room" in data else None
rel = defaultdict(set)
if rr_reg is not None and rr_room is not None:
for r, room in zip(rr_reg.tolist(), rr_room.tolist()):
rel[int(r)].add(int(room))
out["region_room_rel"] = rel
return out

def save_ply_points(path, P, ids=None):
if not HAS_PLY or len(P) == 0:
return
P = np.asarray(P, np.float32)
dt = [("x", "f4"), ("y", "f4"), ("z", "f4")]
if ids is not None:
dt.append(("id", "i4"))
a = np.empty(len(P), dt)
a["x"], a["y"], a["z"] = P[:, 0], P[:, 1], P[:, 2]
if ids is not None:
a["id"] = ids
PlyData([PlyElement.describe(a, "vertex")], text=False).write(path)

def save_ply_lines(path, segments, colors=None):
if not HAS_PLY or len(segments) == 0:
return
all_edges = []
for seg in segments:
seg = np.asarray(seg, float)
n = len(seg)
for i in range(n):
all_edges.append((seg[i], seg[(i + 1) % n]))
if not all_edges:
return
V = np.array([e[0] for e in all_edges] + [e[1] for e in all_edges],
np.float32)
save_ply_points(path, V)

def save_ply_colored(path, P, colors_u8):
if not HAS_PLY or len(P) == 0:
return
P = np.asarray(P, np.float32)
C = np.asarray(colors_u8, np.uint8)
dt = [("x", "f4"), ("y", "f4"), ("z", "f4"),
("red", "u1"), ("green", "u1"), ("blue", "u1")]
a = np.empty(len(P), dt)
a["x"], a["y"], a["z"] = P[:, 0], P[:, 1], P[:, 2]
a["red"], a["green"], a["blue"] = C[:, 0], C[:, 1], C[:, 2]
PlyData([PlyElement.describe(a, "vertex")], text=False).write(path)

def save_ply_mesh(path, V, F, face_ids=None):
if not HAS_PLY:
log.warning("plyfile not installed; PLY mesh skipped")
return
V = np.asarray(V, np.float32)
F = np.asarray(F, np.int32)
if len(V) == 0 or len(F) == 0:
return
vert_dt = [("x", "f4"), ("y", "f4"), ("z", "f4")]
vert = np.empty(len(V), vert_dt)
vert["x"], vert["y"], vert["z"] = V[:, 0], V[:, 1], V[:, 2]
face_dt = [("vertex_indices", "i4", (3,))]
if face_ids is not None:
face_dt.append(("face_id", "i4"))
face = np.empty(len(F), face_dt)
face["vertex_indices"] = F
if face_ids is not None:
face["face_id"] = np.asarray(face_ids, np.int32)
PlyData(
[PlyElement.describe(vert, "vertex"),
PlyElement.describe(face, "face")],
text=False
).write(path)
log.info("PLY mesh written: %s (V=%d F=%d)", path, len(V), len(F))

def save_dae_mesh(path, V, F, face_ids=None):
V = np.asarray(V, np.float64)
F = np.asarray(F, np.int32)
if len(V) == 0 or len(F) == 0:
return
pos_flat = V.reshape(-1).tolist()
pos_str = " ".join(f"{x:.6f}" for x in pos_flat)
tri_flat = F.reshape(-1).tolist()
tri_str = " ".join(str(int(i)) for i in tri_flat)
with open(path, "w", encoding="utf-8") as f:
f.write('<?xml version="1.0" encoding="utf-8"?>\n')
f.write('<COLLADA xmlns="http://www.collada.org/2005/11/COLLADASchema" '
'version="1.4.1">\n')
f.write(' <asset>\n <unit meter="1" name="meter"/>\n'
' <up_axis>Z_UP</up_axis>\n </asset>\n')
f.write(' <library_geometries>\n')
f.write(' <geometry id="geom0" name="frugal">\n')
f.write(' <mesh>\n')
f.write(' <source id="pos0">\n')
f.write(' <float_array id="pos0-array" '
f'count="{len(pos_flat)}">{pos_str}</float_array>\n')
f.write(' <technique_common>\n')
f.write(' <accessor source="#pos0-array" '
f'count="{len(V)}" stride="3">\n')
f.write(' <param name="X" type="float"/>\n')
f.write(' <param name="Y" type="float"/>\n')
f.write(' <param name="Z" type="float"/>\n')
f.write(' </accessor>\n')
f.write(' </technique_common>\n')
f.write(' </source>\n')
f.write(' <vertices id="verts0">\n')
f.write(' <input semantic="POSITION" source="#pos0"/>\n')
f.write(' </vertices>\n')
f.write(' <triangles count="%d" material="mat0">\n' % len(F))
f.write(' <input semantic="VERTEX" source="#verts0" '
'offset="0"/>\n')
f.write(f' <p>{tri_str}</p>\n')
f.write(' </triangles>\n')
f.write(' </mesh>\n </geometry>\n </library_geometries>\n')
f.write(' <library_materials>\n')
f.write(' <material id="mat0" name="frugal_mat">\n')
f.write(' <instance_effect url="#fx0"/>\n')
f.write(' </material>\n </library_materials>\n')
f.write(' <library_effects>\n <effect id="fx0">\n')
f.write(' <profile_COMMON>\n <technique sid="common">\n')
f.write(' <lambert>\n')
f.write(' <diffuse><color>0.7 0.7 0.7 1</color></diffuse>\n')
f.write(' </lambert>\n </technique>\n')
f.write(' </profile_COMMON>\n </effect>\n </library_effects>\n')
f.write(' <library_visual_scenes>\n')
f.write(' <visual_scene id="scene0" name="scene0">\n')
f.write(' <node id="node0" name="frugal_node">\n')
f.write(' <instance_geometry url="#geom0">\n')
f.write(' <bind_material>\n <technique_common>\n')
f.write(' <instance_material symbol="mat0" '
'target="#mat0"/>\n')
f.write(' </technique_common>\n </bind_material>\n')
f.write(' </instance_geometry>\n </node>\n')
f.write(' </visual_scene>\n </library_visual_scenes>\n')
f.write(' <scene>\n <instance_visual_scene url="#scene0"/>\n')
f.write(' </scene>\n</COLLADA>\n')
log.info("DAE written: %s (V=%d F=%d)", path, len(V), len(F))

============================================================

STAGE 1 -- SURFACES

============================================================

def build_surfaces(P, sid, rid, cfg):
out = []
dropped_room = dropped_rms = dropped_pts = 0
for s in np.unique(sid):
if s == 0:
continue
ix = np.where(sid == s)[0]
if len(ix) < cfg.region_min_pts:
dropped_pts += 1
continue
room_frac = float(np.mean(rid[ix] > 0))
if room_frac < cfg.surface_min_room_frac:
dropped_room += 1
continue
Q = P[ix]
n, d, c, rms = fit_plane(Q)
typ = plane_type(n, cfg)
if rms > cfg.plane_rms_max:
dropped_rms += 1
continue
rids = rid[ix]
rids_pos = rids[rids > 0]
room_ids = sorted(set(rids_pos.tolist())) if len(rids_pos) else []
room_hist = {}
if len(rids_pos):
vals, counts = np.unique(rids_pos, return_counts=True)
room_hist = {int(v): int(c) for v, c in zip(vals, counts)}
out.append(dict(
id=len(out), segment_id=int(s),
P=Q, n=n, d=d, c=c, rms=rms, N=len(Q), type=typ,
bmin=Q.min(0), bmax=Q.max(0),
room_frac=room_frac, room_ids=room_ids,
room_histogram_las=room_hist,
member_segments=[int(s)]))
log.info("Structural FRUGAL surfaces: %d", len(out))
log.info(" walls=%d horizontal=%d inclined=%d",
sum(x["type"] == "vertical" for x in out),
sum(x["type"] == "horizontal" for x in out),
sum(x["type"] == "inclined" for x in out))
log.info(" dropped: room=%d pts=%d rms=%d",
dropped_room, dropped_pts, dropped_rms)
return out

def _surface_z_level(s):
if abs(s["n"][2]) < 0.1:
return None
return -s["d"] / s["n"][2]

def merge_coplanar_surfaces_once(S, cfg):
if not S:
return []
parent = np.arange(len(S))

text
def find(x):
    while parent[x] != x:
        parent[x] = parent[parent[x]]
        x = parent[x]
    return x

def union(a, b):
    a, b = find(a), find(b)
    if a != b:
        parent[b] = a

C = np.array([s["c"] for s in S])
tree = cKDTree(C)
pairs = tree.query_pairs(cfg.surface_merge_center_radius)
for i, j in pairs:
    A, B = S[i], S[j]
    if A["type"] != B["type"]:
        continue
    if angle_planes(A["n"], B["n"]) > cfg.surface_merge_angle:
        continue
    d1 = abs(B["c"] @ A["n"] + A["d"])
    d2 = abs(A["c"] @ B["n"] + B["d"])
    if max(d1, d2) > cfg.surface_merge_plane_dist:
        continue
    if A["type"] == "horizontal":
        zA = _surface_z_level(A)
        zB = _surface_z_level(B)
        if zA is not None and zB is not None:
            if abs(zA - zB) > cfg.horizontal_level_tol:
                continue
    gap = np.maximum(
        0., np.maximum(A["bmin"] - B["bmax"], B["bmin"] - A["bmax"]))
    if np.linalg.norm(gap) > cfg.surface_merge_spatial_gap:
        continue
    union(i, j)

groups = {}
for i in range(len(S)):
    groups.setdefault(find(i), []).append(i)

out = []
for ids in groups.values():
    P = np.vstack([S[i]["P"] for i in ids])
    n, d, c, rms = fit_plane(P)
    members = []
    room_fracs = []
    room_ids_all = set()
    room_hist_merged = defaultdict(int)
    for i in ids:
        members.extend(S[i].get(
            "member_segments", [int(S[i]["segment_id"])]))
        room_fracs.append(S[i].get("room_frac", 0.0))
        room_ids_all.update(S[i].get("room_ids", []))
        for k, v in S[i].get("room_histogram_las", {}).items():
            room_hist_merged[k] += v
    if rms > cfg.surface_merge_rms_max and len(ids) > 1:
        for i in ids:
            x = dict(S[i])
            x["id"] = len(out)
            out.append(x)
        continue
    out.append(dict(
        id=len(out),
        segment_id=int(members[0]) if members
            else int(S[ids[0]]["segment_id"]),
        member_segments=sorted(set(int(m) for m in members)),
        P=P, n=n, d=d, c=c, rms=rms, N=len(P),
        type=plane_type(n, cfg),
        bmin=P.min(0), bmax=P.max(0),
        room_frac=float(np.mean(room_fracs)) if room_fracs else 0.0,
        room_ids=sorted(room_ids_all),
        room_histogram_las=dict(room_hist_merged),
    ))
return out

def merge_coplanar_surfaces(S, cfg):
prev_n = len(S)
for k in range(max(1, cfg.surface_merge_passes)):
S = merge_coplanar_surfaces_once(S, cfg)
log.info("Physical surface merge pass %d: %d -> %d",
k + 1, prev_n, len(S))
if len(S) == prev_n:
break
prev_n = len(S)
for i, s in enumerate(S):
s["id"] = i
log.info("Physical surface IDs: total=%d", len(S))
return S

def sampled(P, n):
if len(P) <= n:
return P
return P[np.linspace(0, len(P) - 1, n).astype(int)]

def adjacent(A, B, cfg):
P = sampled(A["P"], cfg.adjacency_sample)
Q = sampled(B["P"], cfg.adjacency_sample)
sep = np.maximum(0, np.maximum(A["bmin"] - B["bmax"],
B["bmin"] - A["bmax"]))
if np.linalg.norm(sep) > cfg.adjacency_dist * 2:
return False
if len(P) > len(Q):
P, Q = Q, P
d, _ = cKDTree(Q).query(P, k=1)
return np.count_nonzero(d <= cfg.adjacency_dist) >= cfg.adjacency_min_pts

def projected_extent(P, p, d, q=.01):
t = (P - p) @ d
return float(np.quantile(t, q)), float(np.quantile(t, 1 - q))

def overlap_extent(A, B, p, d, cfg):
Pa = A["P"]
Pb = B["P"]
qa = Pa[line_dist(Pa, p, d) <= cfg.line_surface_dist]
qb = Pb[line_dist(Pb, p, d) <= cfg.line_surface_dist]
if len(qa) >= cfg.line_surface_min:
ta = (qa - p) @ d
a0, a1 = np.quantile(ta, [cfg.trim_q, 1 - cfg.trim_q])
else:
a0, a1 = projected_extent(A["P"], p, d, cfg.trim_q)
if len(qb) >= cfg.line_surface_min:
tb = (qb - p) @ d
b0, b1 = np.quantile(tb, [cfg.trim_q, 1 - cfg.trim_q])
else:
b0, b1 = projected_extent(B["P"], p, d, cfg.trim_q)
lo = max(a0, b0) - cfg.edge_extent_extension
hi = min(a1, b1) + cfg.edge_extent_extension
if hi - lo < cfg.min_edge_len:
return None
return float(lo), float(hi), len(qa), len(qb)

def edge_evidence(edge_tree, E, p, d, lo, hi, cfg):
if edge_tree is None:
return True, 0, lo, hi
mid = p + d * ((lo + hi) / 2)
R = (hi - lo) / 2 + cfg.edge_search_pad + cfg.edge_line_dist
ids = edge_tree.query_ball_point(mid, R)
if not ids:
return False, 0, lo, hi
Q = E[np.asarray(ids, int)]
dist = line_dist(Q, p, d)
t = (Q - p) @ d
ok = (dist <= cfg.edge_line_dist) &
(t >= lo - cfg.edge_search_pad) & (t <= hi + cfg.edge_search_pad)
t = t[ok]
if len(t) < cfg.edge_min_pts:
return False, len(t), lo, hi
return True, len(t), lo, hi

def build_edges_A(S, E, cfg):
tree = cKDTree(E) if len(E) else None
out = []
tested = adj = 0
for i in range(len(S)):
for j in range(i + 1, len(S)):
A, B = S[i], S[j]
tested += 1
if A["type"] == "horizontal" and B["type"] == "horizontal":
continue
if angle_planes(A["n"], B["n"]) < cfg.parallel_angle_deg:
continue
if not adjacent(A, B, cfg):
continue
adj += 1
z = plane_intersection(A, B)
if z is None:
continue
p, d = z
ext = overlap_extent(A, B, p, d, cfg)
if ext is None:
continue
lo, hi, sa, sb = ext
evidence, ne, lo, hi = edge_evidence(tree, E, p, d, lo, hi, cfg)
wallwall = (A["type"] == "vertical" and B["type"] == "vertical")
confidence = "measured" if evidence else "inferred"
out.append(dict(
source="A", id=len(out),
p=p, d=d,
a=p + lo * d, b=p + hi * d, L=hi - lo,
s0=i, s1=j,
seg0=A["segment_id"], seg1=B["segment_id"],
kind="wall_wall" if wallwall else "wall_horizontal",
support0=sa, support1=sb, edge_support=ne,
confidence=confidence, priority=3))
log.info("[A] plane pairs=%d adjacent=%d edges=%d", tested, adj, len(out))
return out

def filter_A_for_completion(A, min_len=.30):
measured = [e for e in A
if e["confidence"] == "measured" and e["L"] >= min_len]
log.info("[A completion] %d/%d measured structural edges",
len(measured), len(A))
return measured

def simplify_ring(P, tol):
if len(P) < 4:
return P
if HAS_SKIMAGE:
Q = approximate_polygon(P, tolerance=tol)
return Q if len(Q) >= 3 else P
return P

def surface_boundaries(S, cfg, fill_holes=True):
edges_out = []
loops_out = defaultdict(list)

text
if not HAS_SKIMAGE:
    log.warning("[C] scikit-image not installed")
    return edges_out, loops_out

for si, s in enumerate(S):
    if cfg.c_vertical_only and s["type"] != "vertical":
        continue
    P = sampled(s["P"], cfg.c_sample)
    if len(P) < 10:
        continue
    u, v = plane_basis(s["n"])
    Q = P - s["c"]
    uv = np.c_[Q @ u, Q @ v]
    cell = cfg.c_cell
    mn = np.floor(uv.min(0) / cell) * cell
    ij = np.floor((uv - mn) / cell).astype(np.int64)
    W = int(ij[:, 0].max()) + 3
    H = int(ij[:, 1].max()) + 3
    if W <= 0 or H <= 0 or W * H > cfg.c_max_raster_cells:
        continue
    cnt = np.zeros((H, W), np.uint16)
    np.add.at(cnt, (ij[:, 1], ij[:, 0]), 1)
    mask = cnt >= cfg.c_min_cell_pts
    mask = ndimage.binary_closing(
        mask, structure=np.ones((3, 3), bool), iterations=1)
    if fill_holes:
        mask = ndimage.binary_fill_holes(mask)
    lab, nlab = ndimage.label(mask)
    if nlab:
        sizes = np.bincount(lab.ravel())
        good = np.where(sizes >= cfg.c_min_island_cells)[0]
        good = good[good != 0]
        mask = np.isin(lab, good)
    if not mask.any():
        continue
    contours = find_contours(mask.astype(np.uint8), .5)
    for rc in contours:
        if len(rc) < 4:
            continue
        xy = np.c_[
            mn[0] + (rc[:, 1] + .5) * cell,
            mn[1] + (rc[:, 0] + .5) * cell,
        ]
        xy = simplify_ring(xy, cfg.c_simplify)
        if len(xy) < 3:
            continue
        if np.linalg.norm(xy[0] - xy[-1]) > 1e-8:
            xy = np.vstack([xy, xy[0]])
        loop3d = (s["c"] + np.outer(xy[:-1, 0], u)
                  + np.outer(xy[:-1, 1], v))
        loops_out[s["id"]].append(loop3d)
        for k in range(len(xy) - 1):
            a2, b2 = xy[k], xy[k + 1]
            if np.linalg.norm(b2 - a2) < cfg.c_min_len:
                continue
            a3 = s["c"] + a2[0] * u + a2[1] * v
            b3 = s["c"] + b2[0] * u + b2[1] * v
            L = float(np.linalg.norm(b3 - a3))
            if L < cfg.c_min_len:
                continue
            d = unit(b3 - a3)
            edges_out.append(dict(
                source="C", id=len(edges_out),
                p=a3, d=d, a=a3, b=b3, L=L,
                s0=s["id"], s1=-1,
                seg0=s["segment_id"], seg1=-1,
                kind="surface_boundary",
                support0=len(P), support1=0, edge_support=0,
                confidence="boundary", priority=2))

log.info("[C] wall boundary segments=%d", len(edges_out))
n_loops_total = sum(len(v) for v in loops_out.values())
log.info("[C] closed loops: %d (on %d surfaces)",
         n_loops_total, len(loops_out))
return edges_out, loops_out

def intervals_overlap_or_close(a, b, gap):
return max(0., max(a[0], b[0]) - min(a[1], b[1])) <= gap

def interval_on_line(s, p, d):
t0 = float(np.dot(s["a"] - p, d))
t1 = float(np.dot(s["b"] - p, d))
return min(t0, t1), max(t0, t1)

def line_already_represented(x, base, cfg):
mx = (x["a"] + x["b"]) * .5
for c in base:
if angle_dirs(x["d"], c["d"]) > cfg.fuse_angle_deg:
continue
d = line_dist(np.array([mx]), c["a"], c["d"])[0]
if d > cfg.fuse_lateral_dist:
continue
ia = interval_on_line(x, c["a"], c["d"])
ib = interval_on_line(c, c["a"], c["d"])
if intervals_overlap_or_close(ia, ib, cfg.fuse_gap):
return True
return False

def priority_completion(C, A, B, cfg):
out = [dict(s) for s in C]
addA = 0
for a in A:
if not line_already_represented(a, out, cfg):
out.append(dict(a))
addA += 1
for i, s in enumerate(out):
s["id"] = i
log.info("Completion: C=%d +A=%d -> %d", len(C), addA, len(out))
return out

def topology_no_move(edges, tol=.03):
if not edges:
return [], np.empty((0, 3))
E = np.vstack([[s["a"], s["b"]] for s in edges]).astype(float)
tree = cKDTree(E)
uf = UF(len(E))
for a, b in tree.query_pairs(tol):
uf.u(a, b)
groups = {}
for i in range(len(E)):
groups.setdefault(uf.f(i), []).append(i)
roots = list(groups)
rid = {r: i for i, r in enumerate(roots)}
V = np.array([E[groups[r]].mean(0) for r in roots])
out = []
for i, s in enumerate(edges):
x = dict(s)
x["v0"] = rid[uf.f(2 * i)]
x["v1"] = rid[uf.f(2 * i + 1)]
x["a"] = np.asarray(s["a"], float).copy()
x["b"] = np.asarray(s["b"], float).copy()
x["L"] = float(np.linalg.norm(x["b"] - x["a"]))
out.append(x)
return out, V

def preview(edges, step):
P = []
I = []
for i, e in enumerate(edges):
n = max(2, int(np.ceil(e["L"] / step)) + 1)
t = np.linspace(0, 1, n)[:, None]
P.append(e["a"] * (1 - t) + e["b"] * t)
I.extend([i] * n)
if not P:
return np.empty((0, 3)), np.empty(0, np.int32)
return np.vstack(P), np.asarray(I, np.int32)

def save_edges_csv(path, E):
with open(path, "w", newline="", encoding="utf-8") as f:
w = csv.writer(f, delimiter=";")
w.writerow(["edge", "v0", "v1", "source", "kind", "priority",
"confidence", "segment0", "segment1",
"length", "support0", "support1", "roughness_support",
"x1", "y1", "z1", "x2", "y2", "z2"])
for i, e in enumerate(E):
w.writerow([
i, e.get("v0", -1), e.get("v1", -1),
e.get("source", ""), e.get("kind", ""),
e.get("priority", ""), e.get("confidence", ""),
e.get("seg0", -1), e.get("seg1", -1),
f"{e['L']:.4f}",
e.get("support0", 0), e.get("support1", 0),
e.get("edge_support", 0),
*[f"{x:.6f}" for x in np.r_[e["a"], e["b"]]]])

def save_surfaces_csv(path, S):
with open(path, "w", newline="", encoding="utf-8") as f:
w = csv.writer(f, delimiter=";")
w.writerow(["id", "segment_id", "member_segments", "type",
"points", "rms", "room_frac",
"room_ids", "dominant_room_id",
"n_recovered_pts", "n_recovered_ambiguous",
"nx", "ny", "nz", "d"])
for s in S:
members = s.get("member_segments", [s["segment_id"]])
rooms = s.get("room_ids", [])
rec_stats = s.get("recovered_stats", {})
w.writerow([s["id"], s["segment_id"],
"|".join(str(m) for m in members),
s["type"], s["N"], f"{s['rms']:.6f}",
f"{s.get('room_frac', 0.0):.4f}",
"|".join(str(r) for r in rooms),
s.get("dominant_room_id", 0),
rec_stats.get("n_total", 0),
rec_stats.get("n_ambiguous", 0),
*[f"{x:.8f}" for x in s["n"]],
f"{s['d']:.8f}"])

def save_dxf_skeleton(path, E, V):
if not HAS_DXF:
log.warning("DXF skeleton skipped")
return
d = ezdxf.new("R2010")
m = d.modelspace()
for name, col in [("SRC_A_INTERSECTION", 1), ("SRC_C_BOUNDARY", 3),
("SRC_B_ROUGHNESS", 5), ("SRC_MIXED", 6),
("VERTICES", 2)]:
d.layers.add(name, color=col)

text
def pick_layer(e):
    src = e.get("source", "A")
    if src == "A": return "SRC_A_INTERSECTION"
    if src == "C": return "SRC_C_BOUNDARY"
    if src == "B": return "SRC_B_ROUGHNESS"
    return "SRC_MIXED"

for e in E:
    m.add_line(tuple(e["a"]), tuple(e["b"]),
               dxfattribs={"layer": pick_layer(e)})
for p in V:
    m.add_point(tuple(p), dxfattribs={"layer": "VERTICES"})
d.saveas(path)

============================================================

STAGE 2 -- CAD (V23.2)

============================================================

def attach_rooms(S, rooms_data):
rel = None
if rooms_data:
rel = rooms_data.get("region_room_rel", {})
for s in S:
hist = defaultdict(int)
if rel is not None:
for seg in s.get("member_segments", []):
key = int(seg) - 1
for room in rel.get(key, set()):
hist[room] += 1
for room, cnt in s.get("room_histogram_las", {}).items():
hist[room] += cnt
for room in s.get("room_ids", []):
hist[room] += 1
s["room_histogram"] = dict(hist)
s["room_ids"] = sorted(hist.keys())
if hist:
s["dominant_room_id"] = max(hist.items(),
key=lambda kv: kv[1])[0]
else:
s["dominant_room_id"] = 0
s["room_id"] = s["dominant_room_id"]
return S

def build_structural_adjacency(S, cfg):
adj = defaultdict(set)
n = len(S)
samples = [sampled(s["P"], 2000) for s in S]
trees = [cKDTree(p) for p in samples]
adj_dist = cfg.adjacency_dist
min_pts = max(2, cfg.adjacency_min_pts - 1)
for i in range(n):
A = S[i]
for j in range(i + 1, n):
B = S[j]
if A["type"] == "horizontal" and B["type"] == "horizontal":
continue
if angle_planes(A["n"], B["n"]) < cfg.parallel_angle_deg:
continue
sep = np.maximum(0, np.maximum(A["bmin"] - B["bmax"],
B["bmin"] - A["bmax"]))
if np.linalg.norm(sep) > adj_dist * 5:
continue
d, _ = trees[j].query(samples[i], k=1)
if np.count_nonzero(d <= adj_dist) < min_pts:
continue
adj[A["id"]].add(B["id"])
adj[B["id"]].add(A["id"])
return adj

---------- helpers ----------

def _seg_extent_along_line(obs_uv, p_uv, d_uv, band, trim=0.02):
rel = obs_uv - p_uv
t = rel @ d_uv
perp = np.array([-d_uv[1], d_uv[0]])
dist = np.abs(rel @ perp)
near = dist <= band
if np.count_nonzero(near) < 3:
return None
tn = t[near]
return float(np.quantile(tn, trim)), float(np.quantile(tn, 1 - trim))

def _point_in_polygon(pts, poly):
pts = np.asarray(pts, float)
poly = np.asarray(poly, float)
n = len(poly)
inside = np.zeros(len(pts), bool)
x = pts[:, 0]; y = pts[:, 1]
j = n - 1
for i in range(n):
xi, yi = poly[i]; xj, yj = poly[j]
cond = ((yi > y) != (yj > y))
if np.any(cond):
x_int = (xj - xi) * (y - yi) / (yj - yi + 1e-15) + xi
inside ^= cond & (x < x_int)
j = i
return inside

def _polygon_area_uv(poly_uv):
if len(poly_uv) < 3:
return 0.0
uv = np.asarray(poly_uv, float)
return 0.5 * abs(float(np.sum(
uv[:, 0] * np.roll(uv[:, 1], -1) -
np.roll(uv[:, 0], -1) * uv[:, 1])))

def _is_external_boundary_line(obs_uv, p_uv, d_uv, cfg):
perp = np.array([-d_uv[1], d_uv[0]])
rel = obs_uv - p_uv
dist = rel @ perp
along = rel @ d_uv
band_along = max(cfg.structural_band * 2, 0.5)
in_band = np.abs(along) <= band_along
if not np.any(in_band):
return True
d_in = dist[in_band]
pos = np.count_nonzero(d_in > cfg.structural_external_tol)
neg = np.count_nonzero(d_in < -cfg.structural_external_tol)
if pos >= cfg.structural_external_min_wrong and
neg >= cfg.structural_external_min_wrong:
return False
return True

---------- V23.1 boundary provenance ----------

def _edge_provenance(cell_coords_uv, phys_lines, frame_lines, cfg):
provenance = []
n = len(cell_coords_uv)
cos_thr = cfg.provenance_parallel_cos
line_tol = cfg.provenance_line_tol
for i in range(n):
a = cell_coords_uv[i]
b = cell_coords_uv[(i + 1) % n]
edge = b - a
L = np.linalg.norm(edge)
if L < 1e-9:
provenance.append(("unknown", -1))
continue
edge_dir = edge / L
mid = 0.5 * (a + b)
best_kind = "unknown"
best_id = -1
best_dist = float('inf')

text
    for (pa, pb, nid) in phys_lines:
        ld = pb - pa
        ldn = np.linalg.norm(ld)
        if ldn < 1e-9:
            continue
        ld_u = ld / ldn
        if abs(np.dot(edge_dir, ld_u)) < cos_thr:
            continue
        t = np.dot(mid - pa, ld_u)
        proj = pa + t * ld_u
        dist = np.linalg.norm(mid - proj)
        if dist < line_tol and dist < best_dist:
            best_dist = dist
            best_kind = "physical"
            best_id = nid

    if best_kind == "unknown":
        for (pa, pb, fid) in frame_lines:
            ld = pb - pa
            ldn = np.linalg.norm(ld)
            if ldn < 1e-9:
                continue
            ld_u = ld / ldn
            if abs(np.dot(edge_dir, ld_u)) < cos_thr:
                continue
            t = np.dot(mid - pa, ld_u)
            proj = pa + t * ld_u
            dist = np.linalg.norm(mid - proj)
            if dist < line_tol and dist < best_dist:
                best_dist = dist
                best_kind = "frame"
                best_id = fid

    provenance.append((best_kind, best_id))
return provenance

---------- V23.2 RECOVERY PRE-PASS ----------

def recover_unassigned_points(cfg, P, sid, S):
"""
Классифицирует точки sid==0 и раскладывает их по поверхностям.
Плоскости s["n"], s["d"] и s["P"] (core) НЕ меняются.

text
Каждой surface добавляет:
  s["P_recovered"]       : Nx3 — clean (не-ambiguous)
  s["P_recovered_ambig"] : Nx3 — ambiguous (диагностика)
  s["recovered_stats"]   : dict
"""
if not cfg.recovery_enabled:
    for s in S:
        s["P_recovered"] = np.empty((0, 3))
        s["P_recovered_ambig"] = np.empty((0, 3))
        s["recovered_stats"] = {"n_total": 0, "n_ambiguous": 0}
    return {"n_recovered_total": 0, "n_recovered_ambiguous": 0}

log.info("[V23.2] recovery pre-pass")
t0 = time.time() if False else None
import time as _time
t0 = _time.time()

mask_u = sid == 0
P_u = P[mask_u]
n_u = len(P_u)
if n_u == 0:
    for s in S:
        s["P_recovered"] = np.empty((0, 3))
        s["P_recovered_ambig"] = np.empty((0, 3))
        s["recovered_stats"] = {"n_total": 0, "n_ambiguous": 0}
    return {"n_recovered_total": 0, "n_recovered_ambiguous": 0}

mask_a = sid > 0
P_a = P[mask_a]
sid_a = sid[mask_a]

max_sid = int(sid.max()) + 1
sid_to_surf = np.full(max_sid, -1, np.int32)
for s in S:
    for m in s.get("member_segments", []):
        if 0 <= m < max_sid:
            sid_to_surf[m] = s["id"]

log.info("  building KD-trees (A=%d, U=%d)...", len(P_a), n_u)
tree_a = cKDTree(P_a)
tree_u = cKDTree(P_u)

k_u = cfg.recovery_knn_unassigned
chunk = cfg.recovery_chunk
local_normals = np.zeros((n_u, 3), np.float32)
log.info("  computing local normals (k=%d)...", k_u)
for s0 in range(0, n_u, chunk):
    s1 = min(s0 + chunk, n_u)
    _, idx = tree_u.query(P_u[s0:s1], k=k_u, workers=-1)
    pts = P_u[idx]
    c = pts.mean(axis=1, keepdims=True)
    Q = pts - c
    cov = np.einsum("cki,ckj->cij", Q, Q) / max(k_u - 1, 1)
    try:
        w, V = np.linalg.eigh(cov)
        n = V[:, :, 0]
    except np.linalg.LinAlgError:
        n = np.zeros((s1 - s0, 3), np.float32)
        n[:, 2] = 1.0
    ln = np.linalg.norm(n, axis=1, keepdims=True)
    ln[ln < 1e-9] = 1.0
    local_normals[s0:s1] = n / ln

# UV-контекст каждой поверхности
s_ctx = {}
for s in S:
    P_core = s["P"]
    c = P_core.mean(0)
    u, v = surface_basis(s)
    uv = uv_of_points(P_core, s, u, v)
    s_ctx[s["id"]] = {
        "c": c, "u": u, "v": v,
        "uv_min": uv.min(0), "uv_max": uv.max(0),
    }

assigned_surf = np.full(n_u, -1, np.int32)
category = np.full(n_u, 3, np.int8)
ambiguity = np.zeros(n_u, np.int8)

k_a = cfg.recovery_knn_assigned
log.info("  querying kNN assigned (k=%d)...", k_a)
for s0 in range(0, n_u, chunk):
    s1 = min(s0 + chunk, n_u)
    Pc = P_u[s0:s1]
    d_a, i_a = tree_a.query(Pc, k=k_a, workers=-1)
    if k_a == 1:
        d_a = d_a[:, None]; i_a = i_a[:, None]
    neigh_sid = sid_a[i_a]

    dom = np.empty(len(neigh_sid), np.int32)
    amb = np.zeros(len(neigh_sid), np.int8)
    for r in range(len(neigh_sid)):
        vals, cnts = np.unique(neigh_sid[r], return_counts=True)
        dom[r] = vals[np.argmax(cnts)]
        surfs = {int(sid_to_surf[v]) for v in vals.tolist()
                 if 0 <= v < max_sid and sid_to_surf[v] >= 0}
        amb[r] = len(surfs)
    ambiguity[s0:s1] = amb

    uniq_segs = np.unique(dom)
    for seg in uniq_segs:
        if seg <= 0 or seg >= max_sid:
            continue
        sf = int(sid_to_surf[seg])
        if sf < 0:
            continue
        sf_obj = next((x for x in S if x["id"] == sf), None)
        if sf_obj is None:
            continue
        ctx_s = s_ctx[sf]
        m_local = dom == seg
        P_sel = Pc[m_local]
        if len(P_sel) == 0:
            continue
        n_ = sf_obj["n"]; d_ = sf_obj["d"]
        d_plane = np.abs(P_sel @ n_ + d_)
        ln = local_normals[s0:s1][m_local]
        cos = np.clip(np.abs(ln @ n_), -1, 1)
        ang = np.degrees(np.arccos(cos))
        uv = uv_of_points(P_sel, sf_obj, ctx_s["u"], ctx_s["v"])
        lo = ctx_s["uv_min"] - cfg.recovery_uv_extend_m
        hi = ctx_s["uv_max"] + cfg.recovery_uv_extend_m
        in_uv = np.all((uv >= lo) & (uv <= hi), axis=1)

        cat_local = np.full(len(P_sel), 3, np.int8)
        cat_local[(d_plane <= cfg.recovery_dist_recoverable) &
                  (ang <= cfg.recovery_angle_recoverable_deg) &
                  in_uv] = 0
        cat_local[(cat_local == 3) &
                  (d_plane <= cfg.recovery_dist_extendable) &
                  (ang <= cfg.recovery_angle_extendable_deg) &
                  in_uv] = 1
        cat_local[(cat_local == 3) &
                  (d_plane <= cfg.recovery_dist_extendable) &
                  (ang <= 35.0)] = 2

        glob_idx = np.where(m_local)[0] + s0
        assigned_surf[glob_idx] = sf
        category[glob_idx] = cat_local

for s in S:
    s["P_recovered"] = np.empty((0, 3))
    s["P_recovered_ambig"] = np.empty((0, 3))
    s["recovered_stats"] = {"n_total": 0, "n_ambiguous": 0}

m_use = ((category == 0) | (category == 1)) & (assigned_surf >= 0)
idx_use = np.where(m_use)[0]
amb_mask = ambiguity[idx_use] >= 2

if cfg.recovery_exclude_ambiguous:
    clean_idx = idx_use[~amb_mask]
    amb_idx = idx_use[amb_mask]
else:
    clean_idx = idx_use
    amb_idx = np.empty(0, np.int64)

sf_of_clean = assigned_surf[clean_idx]
sf_of_amb = assigned_surf[amb_idx]

for s in S:
    sid_s = s["id"]
    mc = sf_of_clean == sid_s
    ma = sf_of_amb == sid_s
    s["P_recovered"] = P_u[clean_idx[mc]]
    s["P_recovered_ambig"] = P_u[amb_idx[ma]]
    s["recovered_stats"] = {
        "n_total": int(mc.sum()),
        "n_ambiguous": int(ma.sum()),
    }

n_clean = int(len(clean_idx))
n_amb = int(len(amb_idx))
log.info("[V23.2] recovery done: clean=%d, ambiguous=%d, %.1f s",
         n_clean, n_amb, _time.time() - t0)
return {
    "n_recovered_total": n_clean,
    "n_recovered_ambiguous": n_amb,
}

---------- V23.2 ARRANGEMENT ----------

def build_structural_arrangement(s, S, adj, cfg):
"""
V23.2: та же геометрия + dual evidence (core + recovery).
Каждая ячейка получает evidence_source ∈ {"core","mixed","recovered"}.
"""
u, v = surface_basis(s)
P_core = s["P"]
obs_uv = uv_of_points(P_core, s, u, v)
if len(obs_uv) < 3:
return None, "no_obs"

text
P_rec = s.get("P_recovered", np.empty((0, 3)))
rec_uv = uv_of_points(P_rec, s, u, v) if len(P_rec) > 0 \
    else np.empty((0, 2))

neighbors = sorted(adj[s["id"]])
if len(neighbors) < 2:
    return None, "few_neighbors"

seg_lines = []
n_total = n_external = n_internal = 0

for nid in neighbors:
    B = S[nid]
    z = plane_intersection(s, B)
    if z is None:
        continue
    n_total += 1
    p3, d3 = z
    p_uv = uv_of_points(p3, s, u, v)
    d_uv_raw = np.array([d3 @ u, d3 @ v])
    nd = np.linalg.norm(d_uv_raw)
    if nd < 1e-9:
        continue
    d_uv = d_uv_raw / nd
    ext = _seg_extent_along_line(obs_uv, p_uv, d_uv,
                                  band=cfg.structural_band)
    if ext is None:
        continue
    t0, t1 = ext
    t0 -= cfg.structural_extend_m
    t1 += cfg.structural_extend_m
    if t1 - t0 < cfg.structural_min_len:
        continue
    is_ext = _is_external_boundary_line(obs_uv, p_uv, d_uv, cfg)
    if is_ext:
        n_external += 1
    else:
        n_internal += 1
    t0 -= cfg.structural_arrange_extend
    t1 += cfg.structural_arrange_extend
    a_uv = p_uv + t0 * d_uv
    b_uv = p_uv + t1 * d_uv
    seg_lines.append((a_uv, b_uv, nid))

s["_n_neigh"] = n_total
s["_n_ext"] = n_external
s["_n_internal"] = n_internal

if len(seg_lines) < 2:
    return None, "few_lines"
if not HAS_SHAPELY:
    return None, "no_shapely"

bmin_uv = obs_uv.min(0) - cfg.structural_extend_m * 3
bmax_uv = obs_uv.max(0) + cfg.structural_extend_m * 3
frame = [
    (np.array([bmin_uv[0], bmin_uv[1]]),
     np.array([bmax_uv[0], bmin_uv[1]]), -1),
    (np.array([bmax_uv[0], bmin_uv[1]]),
     np.array([bmax_uv[0], bmax_uv[1]]), -2),
    (np.array([bmax_uv[0], bmax_uv[1]]),
     np.array([bmin_uv[0], bmax_uv[1]]), -3),
    (np.array([bmin_uv[0], bmax_uv[1]]),
     np.array([bmin_uv[0], bmin_uv[1]]), -4),
]
all_segs = seg_lines + frame

try:
    ls = [LineString([a.tolist(), b.tolist()]) for a, b, _ in all_segs]
    merged = unary_union(ls)
    cells = list(polygonize(merged))
except Exception as e:
    return None, f"polygonize_{type(e).__name__}"

if not cells:
    return None, "no_cells"

n_core_total = len(obs_uv)
n_rec_total = len(rec_uv)

min_pts_core = cfg.structural_cell_min_pts
min_frac_core = cfg.structural_cell_min_frac
min_pts_rec = cfg.structural_cell_min_recovered_pts
min_frac_rec = cfg.structural_cell_min_recovered_frac
min_core_for_rec = cfg.structural_cell_min_core_pts_rec

supported = []
for poly in cells:
    if poly.is_empty:
        continue
    if not poly.is_valid:
        try:
            poly = poly.buffer(0)
        except Exception:
            continue
        if poly.is_empty:
            continue
    coords = np.asarray(poly.exterior.coords[:-1], float)
    if len(coords) < 3:
        continue

    core_in = _point_in_polygon(obs_uv, coords) if n_core_total > 0 \
        else np.zeros(0, bool)
    n_core = int(core_in.sum())
    frac_core = n_core / max(n_core_total, 1)

    if n_rec_total > 0:
        rec_in = _point_in_polygon(rec_uv, coords)
        n_rec = int(rec_in.sum())
        frac_rec = n_rec / max(n_rec_total, 1)
    else:
        n_rec = 0
        frac_rec = 0.0

    core_ok = (n_core >= min_pts_core and frac_core >= min_frac_core)
    rec_ok = (n_rec >= min_pts_rec and
              frac_rec >= min_frac_rec and
              n_core >= min_core_for_rec)

    if core_ok and n_rec >= min_pts_rec:
        source = "mixed"
    elif core_ok:
        source = "core"
    elif rec_ok:
        source = "recovered"
    else:
        continue

    provenance = _edge_provenance(coords, seg_lines, frame, cfg)
    supported.append({
        "coords": coords,
        "n_pts": n_core,
        "frac": frac_core,
        "area": float(poly.area),
        "provenance": provenance,
        "n_core_hits": n_core,
        "n_rec_hits": n_rec,
        "frac_core": frac_core,
        "frac_rec": frac_rec,
        "evidence_source": source,
    })

if not supported:
    return None, "no_supported"

supported.sort(key=lambda c: c["area"], reverse=True)
return supported, "arrangement"

def _collapse_collinear_uv(uv, tol):
if uv is None or len(uv) < 4:
return uv
n = len(uv)
keep = [True] * n
for i in range(n):
a = uv[(i - 1) % n]; b = uv[i]; c = uv[(i + 1) % n]
ac = c - a
L = np.linalg.norm(ac)
if L < 1e-9:
keep[i] = False; continue
cross2 = (b[0] - a[0]) * (c[1] - a[1]) -
(b[1] - a[1]) * (c[0] - a[0])
if abs(cross2) / L < tol:
keep[i] = False
out = uv[np.array(keep)]
return out if len(out) >= 3 else uv

def _merge_close_vertices_uv(xy, tol):
n = len(xy)
if n < 4:
return xy
keep = [True] * n
for i in range(n):
a = xy[i]; b = xy[(i + 1) % n]
if np.linalg.norm(b - a) < tol:
keep[i] = False
out = xy[np.array(keep)]
return out if len(out) >= 3 else xy

def coordinate_shared_vertices(S, cfg):
items = []
index = {}
for s in S:
for j, cell in enumerate(s.get("cad_cells_data", [])):
loop = cell["loop_3d"]
for i in range(len(loop)):
k = len(items)
items.append((s["id"], j, i, np.array(loop[i], float)))
index[(s["id"], j, i)] = k
if not s.get("cad_cells_data") and s.get("outer_loop_3d") is not None:
loop = s["outer_loop_3d"]
for i in range(len(loop)):
k = len(items)
items.append((s["id"], -1, i, np.array(loop[i], float)))
index[(s["id"], -1, i)] = k
if not items:
return S

text
pts = np.array([it[3] for it in items])
pairs = cKDTree(pts).query_pairs(cfg.vertex_merge_tol)
uf = UF(len(items))
for i, j in pairs:
    uf.u(i, j)
groups = defaultdict(list)
for i in range(len(items)):
    groups[uf.f(i)].append(i)
new_xyz = {}
for _root, idxs in groups.items():
    mean = np.mean([items[i][3] for i in idxs], axis=0)
    for i in idxs:
        new_xyz[i] = mean

n_moved = n_reproj = 0
for s in S:
    for j, cell in enumerate(s.get("cad_cells_data", [])):
        loop = cell["loop_3d"]
        new_loop = np.array(loop, float)
        for i in range(len(loop)):
            k = index[(s["id"], j, i)]
            p = new_xyz[k].copy()
            d_res = float(np.dot(p, s["n"]) + s["d"])
            if abs(d_res) > 1e-9:
                p = p - d_res * s["n"]
                n_reproj += 1
            if not np.allclose(new_loop[i], p, atol=1e-9):
                n_moved += 1
            new_loop[i] = p
        cell["loop_3d"] = new_loop
        u, v = s.get("_u"), s.get("_v")
        if u is not None:
            cell["loop_uv"] = uv_of_points(new_loop, s, u, v)
    if not s.get("cad_cells_data") and s.get("outer_loop_3d") is not None:
        loop = s["outer_loop_3d"]
        new_loop = np.array(loop, float)
        for i in range(len(loop)):
            k = index[(s["id"], -1, i)]
            p = new_xyz[k].copy()
            d_res = float(np.dot(p, s["n"]) + s["d"])
            if abs(d_res) > 1e-9:
                p = p - d_res * s["n"]
                n_reproj += 1
            if not np.allclose(new_loop[i], p, atol=1e-9):
                n_moved += 1
            new_loop[i] = p
        s["outer_loop_3d"] = new_loop

log.info("[V23.2] shared-vertex coordination: moved=%d reprojected=%d",
         n_moved, n_reproj)
return S

---------- V18 FALLBACK ----------

def extract_cad_measured_loops(S, cfg):
if not HAS_SKIMAGE:
return {}
loops_out = defaultdict(list)
st = np.ones((3, 3), bool)
for s in S:
if not cfg.cad_extract_c_for_horizontals and s["type"] == "horizontal":
continue
P = sampled(s["P"], cfg.cad_c_sample)
if len(P) < 10:
continue
if cfg.cad_wall_basis and s["type"] == "vertical":
u, v = wall_basis(s["n"])
else:
u, v = plane_basis(s["n"])
Q = P - s["c"]
uv = np.c_[Q @ u, Q @ v]
cell = cfg.cad_c_cell
mn = np.floor(uv.min(0) / cell) * cell
ij = np.floor((uv - mn) / cell).astype(np.int64)
W = int(ij[:, 0].max()) + 3
H = int(ij[:, 1].max()) + 3
if W <= 0 or H <= 0 or W * H > cfg.cad_c_max_raster_cells:
continue
cnt = np.zeros((H, W), np.uint16)
np.add.at(cnt, (ij[:, 1], ij[:, 0]), 1)
mask = cnt >= 1
if cfg.cad_c_close_iter > 0:
mask = ndimage.binary_closing(mask, structure=st,
iterations=cfg.cad_c_close_iter)
if cfg.cad_c_open_iter > 0:
mask = ndimage.binary_opening(mask, structure=st,
iterations=cfg.cad_c_open_iter)
if cfg.cad_c_fill_holes:
mask = ndimage.binary_fill_holes(mask)
lab, nlab = ndimage.label(mask)
if nlab:
sizes = np.bincount(lab.ravel())
good = np.where(sizes >= cfg.cad_c_min_island_cells)[0]
good = good[good != 0]
mask = np.isin(lab, good)
if not mask.any():
continue
contours = find_contours(mask.astype(np.uint8), .5)
for rc in contours:
if len(rc) < 4:
continue
xy = np.c_[
mn[0] + (rc[:, 1] + .5) * cell,
mn[1] + (rc[:, 0] + .5) * cell,
]
xy = simplify_ring(xy, cfg.cad_c_simplify)
if len(xy) < 3:
continue
if np.linalg.norm(xy[0] - xy[-1]) > 1e-8:
xy = np.vstack([xy, xy[0]])
loop3d = (s["c"] + np.outer(xy[:-1, 0], u)
+ np.outer(xy[:-1, 1], v))
loops_out[s["id"]].append(loop3d)
return loops_out

def _fit_line_pca(pts):
c = pts.mean(0)
Q = pts - c
cov = Q.T @ Q / max(len(pts) - 1, 1)
w, V = np.linalg.eigh(cov)
d = V[:, np.argmax(w)]
return c, d / (np.linalg.norm(d) + 1e-12)

def _intersect_2d_lines(c1, d1, c2, d2):
A = np.array([[d1[0], -d2[0]], [d1[1], -d2[1]]])
if abs(np.linalg.det(A)) < 1e-9:
return None
t = np.linalg.solve(A, c2 - c1)
return c1 + t[0] * d1

def regularize_loop_uv(xy, s, cfg, a_lines_uv=None):
n = len(xy)
if n < 4:
return xy, 0
dirs, lens = [], []
for i in range(n):
a, b = xy[i], xy[(i + 1) % n]
d = b - a
L = np.linalg.norm(d)
if L < 1e-9:
dirs.append(None); lens.append(0.0)
else:
dirs.append(d / L); lens.append(L)
runs = []
cur = None
angle_tol = cfg.cad_reg_run_angle_deg
for i in range(n):
d_i = dirs[i]
if d_i is None:
if cur is not None:
runs.append(cur); cur = None
continue
if cur is None:
cur = {"start": i, "end": i, "dir": d_i.copy(),
"length": lens[i], "idx": [i]}
else:
if angle_dirs(cur["dir"], d_i) <= angle_tol:
cur["end"] = i; cur["length"] += lens[i]
cur["idx"].append(i); cur["dir"] = unit(cur["dir"] + d_i)
else:
runs.append(cur)
cur = {"start": i, "end": i, "dir": d_i.copy(),
"length": lens[i], "idx": [i]}
if cur is not None:
runs.append(cur)
if len(runs) >= 2 and angle_dirs(runs[0]["dir"],
runs[-1]["dir"]) <= angle_tol:
merged = {
"start": runs[-1]["start"], "end": runs[0]["end"],
"dir": unit(runs[-1]["dir"] + runs[0]["dir"]),
"length": runs[-1]["length"] + runs[0]["length"],
"idx": runs[-1]["idx"] + runs[0]["idx"],
}
runs = [merged] + runs[1:-1]
fitted = []
for r in runs:
if r["length"] < cfg.cad_reg_min_run_len:
fitted.append(None); continue
pts = []
for i in r["idx"]:
pts.append(xy[i]); pts.append(xy[(i + 1) % n])
c, d = _fit_line_pca(np.asarray(pts))
fitted.append({"c": c, "d": d, "run": r})
n_snap = 0
if a_lines_uv is not None and len(a_lines_uv) > 0:
for f in fitted:
if f is None: continue
best_a, best_score = None, np.inf
for aa, ab in a_lines_uv:
ad = ab - aa
la = np.linalg.norm(ad)
if la < 1e-9: continue
ad = ad / la
ang = angle_dirs(f["d"], ad)
if ang > cfg.cad_reg_snap_A_angle_deg: continue
t = np.dot(f["c"] - aa, ad)
proj = aa + t * ad
dist = np.linalg.norm(f["c"] - proj)
if dist > cfg.cad_reg_snap_A_dist: continue
score = dist + 0.005 * ang
if score < best_score:
best_score = score; best_a = (aa, ad)
if best_a is not None:
aa, ad = best_a
t = np.dot(f["c"] - aa, ad)
f["c"] = aa + t * ad; f["d"] = ad
n_snap += 1
new_xy = []
m = len(fitted)
for i in range(m):
f_i = fitted[i]
if f_i is None:
r = runs[i]
for k in r["idx"]:
new_xy.append(xy[k])
continue
j = (i + 1) % m
tries = 0
while fitted[j] is None and tries < m:
j = (j + 1) % m; tries += 1
if tries >= m: continue
f_j = fitted[j]
x = _intersect_2d_lines(f_i["c"], f_i["d"], f_j["c"], f_j["d"])
if x is None:
r = runs[i]; k_end = r["idx"][-1]
new_xy.append(0.5 * (xy[k_end] + xy[(k_end + 1) % n]))
else:
new_xy.append(x)
if len(new_xy) < 3:
return xy, n_snap
return np.asarray(new_xy), n_snap

def _classify_outer_inner(loops_3d, s, cfg):
if not loops_3d:
return None
if len(loops_3d) == 1 or not HAS_SHAPELY:
u, v = surface_basis(s)
outer_uv = uv_of_points(loops_3d[0], s, u, v)
return loops_3d[0], [], outer_uv
u, v = surface_basis(s)
polys = []
for loop in loops_3d:
uv = uv_of_points(loop, s, u, v)
try:
p = Polygon(uv)
if not p.is_valid:
p = p.buffer(0)
if p.is_empty: continue
polys.append((p, loop, uv))
except Exception:
continue
if not polys:
outer_uv = uv_of_points(loops_3d[0], s, u, v)
return loops_3d[0], [], outer_uv
polys.sort(key=lambda pw: pw[0].area, reverse=True)
outer_poly, outer_loop, outer_uv = polys[0]
inners = []
for p, loop, uv in polys[1:]:
rp = p.representative_point()
try:
if outer_poly.contains(rp) and p.area >= cfg.cad_c_inner_min_area:
inners.append((loop, uv))
except Exception:
continue
return outer_loop, inners, outer_uv

def assign_edges_to_surfaces(edges, S, cfg):
by_surface = defaultdict(list)
for e in edges:
src = e.get("source", "A")
if src == "C":
sid = e.get("s0", -1)
if sid < 0 or sid >= len(S): continue
s = S[sid]
a3 = project_point_to_plane(np.asarray(e["a"], float), s)
b3 = project_point_to_plane(np.asarray(e["b"], float), s)
by_surface[sid].append(dict(a=a3, b=b3, source="C",
edge_id=e.get("id", -1)))
elif src == "A":
if cfg.cad_use_only_measured_A and e.get("confidence") != "measured":
continue
for sid in (e.get("s0", -1), e.get("s1", -1)):
if sid < 0 or sid >= len(S): continue
s = S[sid]
a3 = project_point_to_plane(np.asarray(e["a"], float), s)
b3 = project_point_to_plane(np.asarray(e["b"], float), s)
by_surface[sid].append(dict(a=a3, b=b3, source="A",
edge_id=e.get("id", -1)))
return by_surface

def _to_uv(edge_list, u, v, c):
out = []
for e in edge_list:
a = e["a"] - c; b = e["b"] - c
out.append((np.array([a @ u, a @ v]), np.array([b @ u, b @ v])))
return out

def fallback_loop_from_cloud(s, cfg):
if not HAS_SHAPELY:
return None
u, v = surface_basis(s)
uv = s["P"] @ np.c[u, v] - np.array([s["c"] @ u, s["c"] @ v])
if len(uv) < 4: return None
if len(uv) > 4000:
idx = np.linspace(0, len(uv) - 1, 4000).astype(int)
uv_h = uv[idx]
else:
uv_h = uv
try:
mp = MultiPoint(uv_h.tolist())
hull = mp.concave_hull(ratio=cfg.cad_concave_ratio)
if hull.geom_type == "Polygon" and len(hull.exterior.coords) >= 4:
coords = np.array(hull.exterior.coords[:-1])
return points_of_uv(coords, s, u, v)
except Exception:
pass
from scipy.spatial import ConvexHull
try:
h = ConvexHull(uv_h)
coords = uv_h[h.vertices]
return points_of_uv(coords, s, u, v)
except Exception:
return None

def _collapse_collinear_3d(pts, tol):
if pts is None or len(pts) < 4: return pts
keep = [True] * len(pts)
n = len(pts)
for i in range(n):
a = pts[(i - 1) % n]; b = pts[i]; c = pts[(i + 1) % n]
ac = c - a; lac = np.linalg.norm(ac)
if lac < 1e-9:
keep[i] = False; continue
dist = np.linalg.norm(np.cross(b - a, ac)) / lac
if dist < tol:
keep[i] = False
out = pts[np.array(keep)]
return out if len(out) >= 3 else pts

def polygonize_surface_uv(lines_uv, min_area):
if not HAS_SHAPELY or not lines_uv: return []
ls = [LineString([a.tolist(), b.tolist()]) for a, b in lines_uv]
try:
network = unary_union(ls)
polys = list(polygonize(network))
except Exception:
return []
polys = [p for p in polys if p.is_valid and p.area >= min_area]
polys.sort(key=lambda p: p.area, reverse=True)
return polys

def build_face_loops_v18_fallback(S, C_loops_stage1, by_surface,
A_lines_by_surface, cfg):
if cfg.cad_use_c_measured_loops:
C_loops_cad = extract_cad_measured_loops(S, cfg)
else:
C_loops_cad = {}
counters = defaultdict(int)
for s in S:
sid = s["id"]
u, v = surface_basis(s)
s["_u"], s["_v"] = u, v
outer = None
inners = []
source = "none"

text
    if sid in C_loops_cad and C_loops_cad[sid]:
        outer_loop, inners_raw, outer_uv = _classify_outer_inner(
            C_loops_cad[sid], s, cfg)
        if outer_loop is not None:
            if cfg.cad_reg_enabled:
                a_lines = A_lines_by_surface.get(sid, [])
                new_uv, _ = regularize_loop_uv(outer_uv, s, cfg, a_lines)
                new_uv = _merge_close_vertices_uv(
                    new_uv, cfg.cad_reg_merge_vertex_tol)
                outer = points_of_uv(new_uv, s, u, v)
                source = "C_cad_reg"
                counters["C_cad_reg"] += 1
                inners = []
                for i_loop, i_uv in inners_raw:
                    if cfg.cad_reg_enabled:
                        i_uv2, _ = regularize_loop_uv(i_uv, s, cfg, None)
                        i_uv2 = _merge_close_vertices_uv(
                            i_uv2, cfg.cad_reg_merge_vertex_tol)
                        inners.append(points_of_uv(i_uv2, s, u, v))
                    else:
                        inners.append(i_loop)
            else:
                outer = outer_loop
                inners = [x[0] for x in inners_raw]
                source = "C_cad_raw"
                counters["C_cad_raw"] += 1

    if outer is None and sid in C_loops_stage1 and C_loops_stage1[sid]:
        outer_loop, _, outer_uv = _classify_outer_inner(
            C_loops_stage1[sid], s, cfg)
        if outer_loop is not None:
            if cfg.cad_reg_enabled:
                a_lines = A_lines_by_surface.get(sid, [])
                new_uv, _ = regularize_loop_uv(outer_uv, s, cfg, a_lines)
                new_uv = _merge_close_vertices_uv(
                    new_uv, cfg.cad_reg_merge_vertex_tol)
                outer = points_of_uv(new_uv, s, u, v)
                source = "C_stage1_reg"
            else:
                outer = outer_loop
                source = "C_stage1"
            counters["C_stage1"] += 1

    if outer is None:
        edge_list = by_surface.get(sid, [])
        if edge_list:
            lines_uv = _to_uv(edge_list, u, v, s["c"])
            polys = polygonize_surface_uv(lines_uv,
                                          cfg.cad_min_polygon_area)
            if polys:
                outer_p = polys[0]
                coords = np.asarray(outer_p.exterior.coords[:-1], float)
                outer = points_of_uv(coords, s, u, v)
                for p in polys[1:]:
                    rp = p.representative_point()
                    try:
                        if outer_p.contains(rp):
                            c2 = np.asarray(p.exterior.coords[:-1], float)
                            inners.append(points_of_uv(c2, s, u, v))
                    except Exception:
                        continue
                source = "skeleton"
                counters["skeleton"] += 1

    if outer is None and cfg.cad_fallback_concave:
        fb = _fallback_loop_from_cloud(s, cfg)
        if fb is not None:
            outer = fb
            source = "fallback_concave"
            counters["fallback"] += 1

    if outer is None:
        counters["missing"] += 1
        s["cad_cells_data"] = []
        s["outer_loop_3d"] = None
        s["loop_source"] = "none"
        continue

    outer = _collapse_collinear_3d(outer, cfg.cad_collinear_tol)
    s["outer_loop_3d"] = outer
    s["outer_loop_uv"] = uv_of_points(outer, s, u, v)
    s["cad_cells_data"] = [{
        "loop_3d": outer,
        "loop_uv": s["outer_loop_uv"],
        "n_pts": len(s["P"]),
        "area_uv": _polygon_area_uv(s["outer_loop_uv"]),
        "provenance": [("fallback", -1)] * len(outer),
        "n_core_hits": len(s["P"]),
        "n_rec_hits": 0,
        "frac_core": 1.0,
        "frac_rec": 0.0,
        "evidence_source": "fallback",
    }]
    s["loop_source"] = source

log.info("[V18 fallback] %s", dict(counters))
return S

def build_face_loops_v23(S, C_loops_stage1, edges, A, cfg):
log.info("=" * 60)
log.info("[V23.2] ARRANGEMENT with recovery evidence")
log.info("=" * 60)

text
if not cfg.structural_enabled:
    by_surface = assign_edges_to_surfaces(edges, S, cfg)
    A_lines = build_A_lines_by_surface(S, A, cfg)
    S = build_face_loops_v18_fallback(S, C_loops_stage1,
                                       by_surface, A_lines, cfg)
    S = coordinate_shared_vertices(S, cfg)
    return S

# Убедимся, что recovery pre-pass уже отработал
for s in S:
    if "P_recovered" not in s:
        s["P_recovered"] = np.empty((0, 3))
    if "P_recovered_ambig" not in s:
        s["P_recovered_ambig"] = np.empty((0, 3))

adj = build_structural_adjacency(S, cfg)
n_adj = sum(len(v) for v in adj.values()) // 2
log.info("[V23.2] adjacency edges: %d", n_adj)

counters = defaultdict(int)
fail_reasons = defaultdict(int)
S_fallback = []
total_cells = 0
n_cells_core = 0
n_cells_recovered = 0
n_cells_mixed = 0

for s in S:
    u, v = surface_basis(s)
    s["_u"], s["_v"] = u, v
    s["cad_cells_data"] = []
    s["outer_loop_3d"] = None

    cells, src = build_structural_arrangement(s, S, adj, cfg)
    if cells is None:
        counters["arrangement_fail"] += 1
        fail_reasons[src] += 1
        S_fallback.append(s)
        continue

    cells_data = []
    for c in cells:
        loop_uv = _merge_close_vertices_uv(
            c["coords"], cfg.cad_reg_merge_vertex_tol)
        prov = c.get("provenance", [])
        loop_uv = _collapse_collinear_uv(loop_uv,
                                          cfg.cad_collinear_tol)
        if len(loop_uv) < 3:
            continue
        loop3d = points_of_uv(loop_uv, s, u, v)
        cells_data.append({
            "loop_3d": loop3d,
            "loop_uv": loop_uv,
            "n_pts": c.get("n_core_hits", c.get("n_pts", 0)),
            "area_uv": _polygon_area_uv(loop_uv),
            "provenance": prov,
            "n_core_hits": c.get("n_core_hits", 0),
            "n_rec_hits": c.get("n_rec_hits", 0),
            "frac_core": c.get("frac_core", 0.0),
            "frac_rec": c.get("frac_rec", 0.0),
            "evidence_source": c.get("evidence_source", "core"),
        })

    if not cells_data:
        counters["arrangement_degenerate"] += 1
        fail_reasons["degenerate"] += 1
        S_fallback.append(s)
        continue

    s["cad_cells_data"] = cells_data
    s["outer_loop_3d"] = cells_data[0]["loop_3d"]
    s["outer_loop_uv"] = cells_data[0]["loop_uv"]
    s["loop_source"] = "structural_arrangement"
    counters["arrangement"] += 1
    counters["cells"] += len(cells_data)
    total_cells += len(cells_data)
    for c in cells_data:
        src_k = c["evidence_source"]
        if src_k == "core":
            n_cells_core += 1
        elif src_k == "recovered":
            n_cells_recovered += 1
        elif src_k == "mixed":
            n_cells_mixed += 1

if S_fallback:
    log.info("[V23.2] raster fallback for %d surfaces", len(S_fallback))
    log.info("[V23.2] arrangement fail reasons: %s", dict(fail_reasons))
    by_surface = assign_edges_to_surfaces(edges, S, cfg)
    A_lines = build_A_lines_by_surface(S, A, cfg)
    build_face_loops_v18_fallback(
        S_fallback, C_loops_stage1, by_surface, A_lines, cfg)
    for s in S_fallback:
        if s.get("cad_cells_data"):
            counters["raster"] += 1
        else:
            counters["missing"] += 1

log.info("[V23.2] loops: %s (total cells=%d)",
         dict(counters), total_cells)
log.info("[V23.2] cell sources: core=%d mixed=%d recovered=%d",
         n_cells_core, n_cells_mixed, n_cells_recovered)
S = coordinate_shared_vertices(S, cfg)
return S

---------- V23.1 CELL EVIDENCE ----------

def _cell_evidence(cell, s, u, v, cfg):
P = s["P"]
if len(P) == 0:
return {"n_pts": 0, "frac": 0.0, "area_uv": 0.0,
"density": 0.0, "rms_dist": 0.0}
obs_uv = uv_of_points(P, s, u, v)
inside = _point_in_polygon(obs_uv, cell["loop_uv"])
n_pts = int(np.count_nonzero(inside))
frac = n_pts / len(P)
area_uv = _polygon_area_uv(cell["loop_uv"])
density = n_pts / max(1e-6, area_uv)
if n_pts > 0:
P_in = P[inside]
if len(P_in) > cfg.evidence_sample:
idx = np.linspace(0, len(P_in) - 1,
cfg.evidence_sample).astype(int)
P_in = P_in[idx]
dists = np.abs(P_in @ s["n"] + s["d"])
rms = float(np.sqrt(np.mean(dists * dists)))
else:
rms = 0.0
return {
"n_pts": n_pts, "frac": frac,
"area_uv": float(area_uv),
"density": float(density),
"rms_dist": rms,
}

---------- CAD faces / tessellation ----------

def _make_wire_from_loop(loop3d):
edges = []
n = len(loop3d)
for i in range(n):
a = loop3d[i]; b = loop3d[(i + 1) % n]
if np.linalg.norm(b - a) < 1e-9:
continue
edges.append(cq.Edge.makeLine(cq.Vector(*a), cq.Vector(*b)))
if len(edges) < 3:
return None
try:
return cq.Wire.assembleEdges(edges)
except Exception:
return None

def validate_loop(loop, s, cfg):
if loop is None or len(loop) < 3:
return False, "too_few", {}
P = np.asarray(loop, float)
res = np.abs(P @ s["n"] + s["d"])
planarity = float(res.max())
u, v = surface_basis(s)
uv = uv_of_points(P, s, u, v)
area = polygon_area_uv(uv)
meta = {"planarity": planarity, "area": area, "n_verts": len(P)}
if planarity > cfg.cad_plane_res_tol:
return False, f"not_planar{planarity:.4f}", meta
if area < cfg.cad_min_face_area:
return False, f"area{area:.4f}", meta
return True, "ok", meta

def build_cad_faces(S, cfg):
if not HAS_CQ:
log.warning("[CAD] cadquery not installed")
return S
n_ok = 0
n_rejected = 0
n_faces_total = 0
n_surfaces_with_faces = 0
n_recovered_cells = 0
for s in S:
u, v = s.get("_u"), s.get("_v")
if u is None:
u, v = surface_basis(s)
faces_data = s.get("cad_cells_data", [])
if not faces_data:
loop = s.get("outer_loop_3d")
if loop is not None and len(loop) >= 3:
faces_data = [{"loop_3d": loop,
"loop_uv": s.get("outer_loop_uv"),
"n_pts": len(s.get("P", [])),
"area_uv": 0.0,
"provenance": [],
"evidence_source": "fallback"}]
else:
s["cad_faces_data"] = []
s["face"] = None
s["face_reject"] = "no_loop"
continue

text
    faces_out = []
    for cell in faces_data:
        loop = cell["loop_3d"]
        ok, reason, meta = _validate_loop(loop, s, cfg)
        evidence = _cell_evidence(cell, s, u, v, cfg)
        if not ok:
            n_rejected += 1
            continue
        outer = _make_wire_from_loop(loop)
        if outer is None:
            n_rejected += 1
            continue
        try:
            f = cq.Face.makeFromWires(outer)
            if f.Area() < cfg.cad_min_face_area:
                n_rejected += 1
                continue
            src = cell.get("evidence_source", "core")
            if src == "recovered":
                n_recovered_cells += 1
            faces_out.append({
                "face": f,
                "loop_3d": loop,
                "loop_uv": cell["loop_uv"],
                "meta": meta,
                "n_pts": cell.get("n_pts", 0),
                "provenance": cell.get("provenance", []),
                "evidence": evidence,
                "evidence_source": src,
            })
            n_ok += 1
        except Exception as e:
            log.debug("[CAD] face %d failed: %s", s["id"], e)
            n_rejected += 1
    s["cad_faces_data"] = faces_out
    if faces_out:
        n_surfaces_with_faces += 1
        s["face"] = max(faces_out,
                        key=lambda fc: fc["meta"]["area"])["face"]
        s["face_reject"] = ""
    else:
        s["face"] = None
        s["face_reject"] = "all_cells_rejected"
    n_faces_total += len(faces_out)
log.info("[CAD] built faces: %d/%d surfaces (%d cells total, "
         "recovered cells=%d, rejected %d)",
         n_surfaces_with_faces, len(S), n_faces_total,
         n_recovered_cells, n_rejected)
return S

def _triangulate_face(face, cfg):
try:
verts, tris = face.tessellate(
cfg.cad_tess_tolerance, cfg.cad_tess_angular)
except Exception:
return None, None
if not verts or not tris:
return None, None
V = np.array([[v.x, v.y, v.z] for v in verts], np.float64)
F = []
for t in tris:
idx = list(t)
if len(idx) == 3:
F.append(idx)
elif len(idx) == 4:
F.append([idx[0], idx[1], idx[2]])
F.append([idx[0], idx[2], idx[3]])
if not F:
return None, None
return V, np.asarray(F, np.int32)

def build_mesh_from_faces(S, cfg):
vmap = {}
V_out = []
F_out = []
faceid_out = []

text
def get_vidx(v):
    key = (round(v[0], 5), round(v[1], 5), round(v[2], 5))
    if key in vmap:
        return vmap[key]
    idx = len(V_out)
    V_out.append([v[0], v[1], v[2]])
    vmap[key] = idx
    return idx

n_faces_used = 0
for s in S:
    faces_data = s.get("cad_faces_data", [])
    if not faces_data and s.get("face") is not None:
        faces_data = [{"face": s["face"]}]
    for fc in faces_data:
        face = fc.get("face")
        if face is None:
            continue
        V, F = _triangulate_face(face, cfg)
        if V is None:
            continue
        n_faces_used += 1
        for tri in F:
            a, b, c = V[tri[0]], V[tri[1]], V[tri[2]]
            ia = get_vidx(a); ib = get_vidx(b); ic = get_vidx(c)
            if ia == ib or ib == ic or ia == ic:
                continue
            F_out.append([ia, ib, ic])
            faceid_out.append(s["id"])

if not V_out or not F_out:
    log.warning("[CAD] mesh empty")
    return (np.empty((0, 3)),
            np.empty((0, 3), np.int32),
            np.empty(0, np.int32))

V_out = np.asarray(V_out, np.float64)
F_out = np.asarray(F_out, np.int32)
faceid_out = np.asarray(faceid_out, np.int32)
avg_tri = float(len(F_out)) / max(1, n_faces_used)
log.info("[CAD] mesh: faces_used=%d V=%d F=%d (avg tri/face=%.2f)",
         n_faces_used, len(V_out), len(F_out), avg_tri)
return V_out, F_out, faceid_out

def export_ply(S, cfg, out_dir):
V, F, fid = build_mesh_from_faces(S, cfg)
if len(V) == 0: return
save_ply_mesh(os.path.join(out_dir, "model.ply"), V, F, fid)

def export_dae(S, cfg, out_dir):
V, F, fid = build_mesh_from_faces(S, cfg)
if len(V) == 0: return
save_dae_mesh(os.path.join(out_dir, "model.dae"), V, F, fid)

def export_dxf_stage2(S, cfg, out_dir):
if not HAS_DXF:
return
d = ezdxf.new("R2010")
m = d.modelspace()
d.layers.add("CAD_EDGES", color=7)
d.layers.add("CAD_FACES", color=4)
d.layers.add("CAD_EDGES_PHYS", color=3)
d.layers.add("CAD_EDGES_FRAME", color=1)
d.layers.add("CAD_EDGES_RECOVERED", color=6) # V23.2
n_edges = 0
n_faces = 0
n_recovered = 0
for s in S:
for fc in s.get("cad_faces_data", []):
loop = fc.get("loop_3d")
if loop is None or len(loop) < 3:
continue
prov = fc.get("provenance", [])
src = fc.get("evidence_source", "core")
n = len(loop)
for i in range(n):
a = tuple(map(float, loop[i]))
b = tuple(map(float, loop[(i + 1) % n]))
if src == "recovered":
kind = "CAD_EDGES_RECOVERED"
n_recovered += 1
else:
kind = "CAD_EDGES"
if i < len(prov):
pk = prov[i][0]
if pk == "physical":
kind = "CAD_EDGES_PHYS"
elif pk == "frame":
kind = "CAD_EDGES_FRAME"
try:
m.add_line(a, b, dxfattribs={"layer": kind})
n_edges += 1
except Exception:
pass
if not cfg.cad_dxf_3dface:
continue
face = fc.get("face")
if face is None:
continue
V, F = _triangulate_face(face, cfg)
if V is None:
continue
for tri in F:
a = tuple(map(float, V[tri[0]]))
b = tuple(map(float, V[tri[1]]))
c = tuple(map(float, V[tri[2]]))
try:
m.add_3dface([a, b, c, c],
dxfattribs={"layer": "CAD_FACES"})
n_faces += 1
except Exception:
pass
d.saveas(os.path.join(out_dir, "cad_3d.dxf"))
log.info("[CAD] DXF written: edges=%d faces=%d recovered_edges=%d",
n_edges, n_faces, n_recovered)

def compute_consistency_metric(S, edges, cfg):
cad_segments = []
for s in S:
for fc in s.get("cad_faces_data", []):
loop = fc.get("loop_3d")
if loop is None or len(loop) < 2: continue
n = len(loop)
for i in range(n):
cad_segments.append((loop[i], loop[(i + 1) % n]))
if not cad_segments:
log.info("[CAD consistency] no segments")
return
cad_A = np.array([a for a, b in cad_segments])
cad_B = np.array([b for a, b in cad_segments])

text
def pt_seg_dist_vec(P, A, B):
    AB = B - A
    denom = (AB * AB).sum(1) + 1e-12
    t = np.clip(((P - A) * AB).sum(1) / denom, 0, 1)
    proj = A + AB * t[:, None]
    return np.linalg.norm(P - proj, axis=1)

sample = cfg.cad_consistency_sample
if len(edges) > sample:
    idx = np.linspace(0, len(edges) - 1, sample).astype(int)
    edges_s = [edges[i] for i in idx]
else:
    edges_s = edges

all_d, C_d, A_d, vert_d, horiz_d = [], [], [], [], []
for e in edges_s:
    mid = 0.5 * (np.asarray(e["a"]) + np.asarray(e["b"]))
    d = pt_seg_dist_vec(mid, cad_A, cad_B).min()
    all_d.append(d)
    if e.get("source") == "C": C_d.append(d)
    elif e.get("source") == "A": A_d.append(d)
    s0 = e.get("s0", -1)
    if 0 <= s0 < len(S):
        t = S[s0]["type"]
        if t == "vertical": vert_d.append(d)
        elif t == "horizontal": horiz_d.append(d)

def report(name, arr):
    if len(arr) == 0:
        log.info("[CAD consistency] %s: no data", name); return
    a = np.asarray(arr)
    unmatched = np.count_nonzero(a > 0.10)
    log.info("[CAD consistency] %s: n=%d p50=%.4f p90=%.4f "
             "p95=%.4f max=%.4f unmatched>0.10=%d (%.1f%%)",
             name, len(a), np.percentile(a, 50), np.percentile(a, 90),
             np.percentile(a, 95), a.max(),
             unmatched, 100.0 * unmatched / len(a))

log.info("=" * 60)
report("ALL", all_d); report("Source C", C_d); report("Source A", A_d)
report("Vertical", vert_d); report("Horizontal", horiz_d)
log.info("=" * 60)

def build_A_lines_by_surface(S, A, cfg):
by_surface = defaultdict(list)
for e in A:
if e.get("confidence") != "measured": continue
if e.get("L", 0) < cfg.a_completion_min_len: continue
for sid in (e.get("s0", -1), e.get("s1", -1)):
if sid < 0 or sid >= len(S): continue
s = S[sid]
u, v = surface_basis(s)
a3 = project_point_to_plane(np.asarray(e["a"], float), s)
b3 = project_point_to_plane(np.asarray(e["b"], float), s)
qa = a3 - s["c"]; qb = b3 - s["c"]
by_surface[sid].append((np.array([qa @ u, qa @ v]),
np.array([qb @ u, qb @ v])))
return by_surface

---------- V23.1 SHARED EDGE VALIDATION ----------

def validate_shared_edges(S, cfg):
all_edges = []
for s in S:
for ci, fc in enumerate(s.get("cad_faces_data", [])):
loop = fc.get("loop_3d")
if loop is None:
continue
n = len(loop)
for i in range(n):
a = np.asarray(loop[i], float)
b = np.asarray(loop[(i + 1) % n], float)
L = np.linalg.norm(b - a)
if L < cfg.shared_edge_min_len:
continue
all_edges.append({
"s_id": s["id"],
"cell_idx": ci,
"a": a, "b": b,
"mid": 0.5 * (a + b),
"len": L,
"dir": (b - a) / L,
})
if not all_edges:
return [], {}

text
report = []
n_matched = 0
dists_matched = []
for i, e1 in enumerate(all_edges):
    best_dist = float('inf')
    best_j = -1
    for j, e2 in enumerate(all_edges):
        if j == i:
            continue
        if e2["s_id"] == e1["s_id"] and e2["cell_idx"] == e1["cell_idx"]:
            continue
        if abs(np.dot(e1["dir"], e2["dir"])) < 0.95:
            continue
        diff = e2["mid"] - e1["mid"]
        perp = diff - np.dot(diff, e1["dir"]) * e1["dir"]
        dist = float(np.linalg.norm(perp))
        if dist < best_dist:
            best_dist = dist
            best_j = j
    matched = best_j >= 0 and best_dist < cfg.shared_edge_tol
    if matched:
        n_matched += 1
        dists_matched.append(best_dist)
    report.append({
        "s_id": e1["s_id"],
        "cell_idx": e1["cell_idx"],
        "edge_len": e1["len"],
        "min_shared_dist": best_dist if best_j >= 0 else -1.0,
        "matched": matched,
    })
summary = {
    "n_edges": len(report),
    "n_matched": n_matched,
    "frac_matched": n_matched / max(1, len(report)),
    "p50_dist": float(np.percentile(dists_matched, 50))
        if dists_matched else 0.0,
    "p90_dist": float(np.percentile(dists_matched, 90))
        if dists_matched else 0.0,
    "max_dist": float(max(dists_matched)) if dists_matched else 0.0,
}
return report, summary

def save_shared_edges_diag_csv(path, report):
with open(path, "w", newline="", encoding="utf-8") as f:
w = csv.writer(f, delimiter=";")
w.writerow(["s_id", "cell_idx", "edge_len",
"min_shared_dist", "matched"])
for r in report:
w.writerow([
r["s_id"], r["cell_idx"],
f"{r['edge_len']:.4f}",
f"{r['min_shared_dist']:.4f}",
1 if r["matched"] else 0])
log.info("[V23.2] shared-edges diagnostic written: %s", path)

---------- V23.1 UNSUPPORTED CELLS ----------

def detect_unsupported_cells(S, cfg):
flagged = []
for s in S:
for ci, fc in enumerate(s.get("cad_faces_data", [])):
ev = fc.get("evidence", {})
n_pts = ev.get("n_pts", 0)
density = ev.get("density", 0.0)
area = ev.get("area_uv", 0.0)
if n_pts < cfg.unsupported_min_pts or
density < cfg.unsupported_min_density:
flagged.append({
"s_id": s["id"],
"cell_idx": ci,
"n_pts": n_pts,
"density": density,
"area": area,
"loop_3d": fc.get("loop_3d"),
})
return flagged

def save_unsupported_cells_ply(path, flagged):
if not HAS_PLY or not flagged:
return
segs = [f["loop_3d"] for f in flagged if f["loop_3d"] is not None]
if not segs:
return
save_ply_lines(path, segs)
log.info("[V23.2] unsupported cells diagnostic written: %s (%d)",
path, len(segs))

def save_unsupported_cells_csv(path, flagged):
with open(path, "w", newline="", encoding="utf-8") as f:
w = csv.writer(f, delimiter=";")
w.writerow(["s_id", "cell_idx", "n_pts", "density", "area_uv"])
for x in flagged:
w.writerow([x["s_id"], x["cell_idx"], x["n_pts"],
f"{x['density']:.2f}", f"{x['area']:.4f}"])
log.info("[V23.2] unsupported cells CSV: %s", path)

---------- V23.2 CAD faces diagnostic ----------

def save_cad_faces_diag_csv(path, S, cfg):
with open(path, "w", newline="", encoding="utf-8") as f:
w = csv.writer(f, delimiter=";")
w.writerow([
"id", "type", "N", "rms",
"loop_source", "n_cells",
"n_neigh", "n_ext", "n_internal",
"cell_idx",
"cell_n_pts", "cell_frac",
"cell_n_core_hits", "cell_n_rec_hits",
"cell_frac_core", "cell_frac_rec",
"cell_evidence_source",
"cell_area_uv", "cell_density", "cell_rms_dist",
"loop_verts", "planarity", "area_3d",
"prov_physical", "prov_frame", "prov_unknown",
"tri_count", "tri_area_sum",
"face_ok", "reject_reason", "bbox_diag"])
for s in S:
faces = s.get("cad_faces_data", [])
cells_uv = s.get("cad_cells_data", [])
bbox_diag = float(np.linalg.norm(s['bmax'] - s['bmin']))
if not faces:
w.writerow([
s["id"], s["type"], s["N"], f"{s['rms']:.4f}",
s.get("loop_source", "none"), 0,
s.get("_n_neigh", 0), s.get("_n_ext", 0),
s.get("_n_internal", 0),
-1, 0, 0.0, 0, 0, 0.0, 0.0, "",
0.0, 0.0, 0.0,
0, 0.0, 0.0, 0, 0, 0, 0, 0.0,
0, s.get("face_reject", ""), f"{bbox_diag:.4f}"])
continue
for j, fc in enumerate(faces):
loop = fc.get("loop_3d")
meta = fc.get("meta", {})
ev = fc.get("evidence", {})
prov = fc.get("provenance", [])
cell = cells_uv[j] if j < len(cells_uv) else {}
n_phys = sum(1 for p in prov if p[0] == "physical")
n_frame = sum(1 for p in prov if p[0] == "frame")
n_unk = sum(1 for p in prov if p[0] == "unknown")
tri_count = 0
tri_area = 0.0
face = fc.get("face")
if face is not None:
V, F = _triangulate_face(face, cfg)
if V is not None and F is not None:
tri_count = len(F)
for t in F:
a, b, c = V[t[0]], V[t[1]], V[t[2]]
tri_area += 0.5 * float(np.linalg.norm(
np.cross(b - a, c - a)))
w.writerow([
s["id"], s["type"], s["N"], f"{s['rms']:.4f}",
s.get("loop_source", "none"),
len(cells_uv),
s.get("_n_neigh", 0), s.get("_n_ext", 0),
s.get("_n_internal", 0),
j,
ev.get("n_pts", 0), f"{ev.get('frac', 0.0):.4f}",
cell.get("n_core_hits", 0),
cell.get("n_rec_hits", 0),
f"{cell.get('frac_core', 0.0):.4f}",
f"{cell.get('frac_rec', 0.0):.4f}",
cell.get("evidence_source", ""),
f"{ev.get('area_uv', 0.0):.4f}",
f"{ev.get('density', 0.0):.2f}",
f"{ev.get('rms_dist', 0.0):.6f}",
len(loop) if loop is not None else 0,
f"{meta.get('planarity', 0.0):.6f}",
f"{meta.get('area', 0.0):.6f}",
n_phys, n_frame, n_unk,
tri_count, f"{tri_area:.6f}",
1 if face is not None else 0,
s.get("face_reject", ""),
f"{bbox_diag:.4f}"])
log.info("[CAD] faces diagnostic written: %s", path)

def run_cad_stage2(S, edges, A, C_loops_stage1, rooms_data, cfg,
P=None, sid=None):
log.info("=" * 60)
log.info("STAGE 2: CAD (V23.2, evidence-augmented)")
log.info("=" * 60)

text
for s in S:
    s["_rooms_data"] = rooms_data

S = attach_rooms(S, rooms_data)
n_with_room = sum(1 for s in S if s.get("room_id", 0) > 0)
log.info("[CAD] surfaces with dominant room: %d/%d",
         n_with_room, len(S))

# ---- recovery pre-pass ----
if P is not None and sid is not None and cfg.recovery_enabled:
    rec_stats = recover_unassigned_points(cfg, P, sid, S)
    log.info("[V23.2] recovery: clean=%d ambiguous=%d",
             rec_stats.get("n_recovered_total", 0),
             rec_stats.get("n_recovered_ambiguous", 0))
else:
    log.info("[V23.2] recovery skipped (no P/sid passed)")
    for s in S:
        s["P_recovered"] = np.empty((0, 3))
        s["P_recovered_ambig"] = np.empty((0, 3))
        s["recovered_stats"] = {"n_total": 0, "n_ambiguous": 0}

S = build_face_loops_v23(S, C_loops_stage1, edges, A, cfg)
S = build_cad_faces(S, cfg)
save_cad_faces_diag_csv(
    os.path.join(cfg.out_dir, "cad_faces_diag.csv"), S, cfg)

# shared-edge validation
shared_report, shared_summary = validate_shared_edges(S, cfg)
if shared_report:
    save_shared_edges_diag_csv(
        os.path.join(cfg.out_dir, "shared_edges_diag.csv"),
        shared_report)
    log.info("[V23.2] shared-edge: n=%d matched=%d (%.1f%%) "
             "p50=%.4f p90=%.4f max=%.4f",
             shared_summary["n_edges"],
             shared_summary["n_matched"],
             100.0 * shared_summary["frac_matched"],
             shared_summary["p50_dist"],
             shared_summary["p90_dist"],
             shared_summary["max_dist"])

# unsupported cells
flagged = detect_unsupported_cells(S, cfg)
if flagged:
    save_unsupported_cells_csv(
        os.path.join(cfg.out_dir, "unsupported_cells.csv"), flagged)
    save_unsupported_cells_ply(
        os.path.join(cfg.out_dir, "unsupported_cells.ply"), flagged)
    log.info("[V23.2] unsupported cells flagged: %d", len(flagged))

# loops
loop_pts, loop_ids = [], []
for s in S:
    for fc in s.get("cad_faces_data", []):
        loop = fc.get("loop_3d")
        if loop is None: continue
        loop_pts.append(loop)
        loop_ids.extend([s["id"]] * len(loop))
if loop_pts and HAS_PLY:
    save_ply_points(os.path.join(cfg.out_dir, "07_LOOPS_FINAL.ply"),
                    np.vstack(loop_pts),
                    np.asarray(loop_ids, np.int32))

if cfg.cad_export_ply:
    export_ply(S, cfg, cfg.out_dir)
if cfg.cad_export_dae:
    export_dae(S, cfg, cfg.out_dir)
if cfg.cad_export_dxf:
    export_dxf_stage2(S, cfg, cfg.out_dir)

compute_consistency_metric(S, edges, cfg)
return S, {}

============================================================

PIPELINE

============================================================

def run(cfg):
os.makedirs(cfg.out_dir, exist_ok=True)

text
P, sid, rid = load_frugal(cfg.frugal_las)
R = load_edges(cfg.edge_las)
rooms_data = load_rooms_npz(cfg.rooms_npz)
log.info("FRUGAL points=%d roughness edge points=%d", len(P), len(R))

S = build_surfaces(P, sid, rid, cfg)
if not S:
    raise RuntimeError("No structural surfaces")

PP, II = [], []
for s in S:
    PP.append(s["P"]); II.extend([s["id"]] * len(s["P"]))
if HAS_PLY:
    save_ply_points(os.path.join(
        cfg.out_dir, "01_STRUCTURAL_SURFACES.ply"),
        np.vstack(PP), np.asarray(II, np.int32))

S = merge_coplanar_surfaces(S, cfg)

PP, II = [], []
for s in S:
    PP.append(s["P"]); II.extend([s["id"]] * len(s["P"]))
if HAS_PLY:
    save_ply_points(os.path.join(
        cfg.out_dir, "01b_PHYSICAL_SURFACES.ply"),
        np.vstack(PP), np.asarray(II, np.int32))

C, C_loops_stage1 = surface_boundaries(S, cfg, fill_holes=True)
Q, I = preview(C, cfg.preview_step)
if HAS_PLY and len(Q):
    save_ply_points(os.path.join(cfg.out_dir, "02c_SOURCE_C.ply"), Q, I)
C = [dict(x) for x in C]

A = build_edges_A(S, R, cfg)
Q, I = preview(A, cfg.preview_step)
if HAS_PLY and len(Q):
    save_ply_points(os.path.join(cfg.out_dir, "02a_SOURCE_A.ply"), Q, I)

B = []
if cfg.run_A_completion:
    A_add = filter_A_for_completion(A, cfg.a_completion_min_len)
else:
    A_add = []
F = priority_completion(C, A_add, [], cfg)
Q, I = preview(F, cfg.preview_step)
if HAS_PLY and len(Q):
    save_ply_points(os.path.join(cfg.out_dir, "03_COMPLETED_LINES.ply"),
                    Q, I)

edges, V = topology_no_move(F, tol=cfg.topo_tol)
log.info("Stage 1 FINAL: surfaces=%d lines=%d edges=%d vertices=%d",
         len(S), len(F), len(edges), len(V))

Q, I = preview(edges, cfg.preview_step)
if HAS_PLY and len(Q):
    save_ply_points(os.path.join(cfg.out_dir, "05_FINAL_SKELETON.ply"),
                    Q, I)
if HAS_PLY and len(V):
    save_ply_points(os.path.join(cfg.out_dir, "06_FINAL_VERTICES.ply"),
                    V, np.arange(len(V), dtype=np.int32))

C_only = [e for e in edges if e.get("source") == "C"]
A_only = [e for e in edges if e.get("source") == "A"]
if HAS_PLY and C_only:
    Q, I = preview(C_only, cfg.preview_step)
    save_ply_points(os.path.join(cfg.out_dir, "05a_FINAL_C_ONLY.ply"),
                    Q, I)
if HAS_PLY and A_only:
    Q, I = preview(A_only, cfg.preview_step)
    save_ply_points(os.path.join(cfg.out_dir, "05b_FINAL_A_ONLY.ply"),
                    Q, I)

save_surfaces_csv(os.path.join(cfg.out_dir, "surfaces.csv"), S)
save_edges_csv(os.path.join(cfg.out_dir, "edges.csv"), edges)
save_dxf_skeleton(os.path.join(cfg.out_dir, "skeleton_3d.dxf"),
                  edges, V)

if cfg.run_cad_stage2:
    if not HAS_CQ:
        log.warning("cadquery NOT installed")
    if not HAS_SHAPELY:
        log.warning("shapely NOT installed")
    try:
        run_cad_stage2(S, edges, A, C_loops_stage1, rooms_data, cfg,
                        P=P, sid=sid)
    except Exception as e:
        log.exception("Stage 2 failed: %s", e)

log.info("DONE -> %s", cfg.out_dir)

============================================================

MAIN

============================================================

if name == "main":
cfg = Config()

text
base = r"C:\Users\dimma\PESN\RESULTS\600"
out_base = r"C:\Users\dimma\PESN\RESULTS\600"

cfg.frugal_las = os.path.join(base, "FRUGAL_full.las")
cfg.edge_las = os.path.join(base, "roughness_edges8.las")
cfg.rooms_npz = os.path.join(base, "FRUGAL_rooms_epoch0.npz")
cfg.out_dir = os.path.join(out_base, "FRUGAL_CAD_V23_3")

cfg.run_roughness_b = False
cfg.run_A_completion = True
cfg.run_C_merge = False

cfg.surface_min_room_frac = .20
cfg.region_min_pts = 100
cfg.plane_rms_max = .045
cfg.surface_merge_angle = 5.0
cfg.surface_merge_plane_dist = .08
cfg.surface_merge_spatial_gap = .80
cfg.surface_merge_center_radius = 4.0
cfg.surface_merge_rms_max = .08
cfg.surface_merge_passes = 2
cfg.adjacency_dist = .25
cfg.adjacency_min_pts = 4
cfg.adjacency_sample = 4000
cfg.line_surface_dist = .16
cfg.line_surface_min = 4
cfg.min_edge_len = .12
cfg.edge_extent_extension = .12
cfg.edge_line_dist = .15
cfg.edge_search_pad = .50
cfg.edge_min_pts = 4
cfg.c_cell = .08
cfg.c_min_cell_pts = 1
cfg.c_simplify = .10
cfg.c_min_len = .30
cfg.c_sample = 15000
cfg.c_max_raster_cells = 4_000_000
cfg.c_min_island_cells = 4
cfg.c_vertical_only = True
cfg.a_completion_min_len = .30
cfg.fuse_angle_deg = 4.0
cfg.fuse_lateral_dist = .05
cfg.fuse_gap = .25
cfg.fuse_min_len = .15
cfg.topo_tol = .03
cfg.triple_plane_tol = .20
cfg.endpoint_vertex_tol = .15
cfg.triple_endpoint_extension = .30

cfg.run_cad_stage2 = True
cfg.cad_export_ply = True
cfg.cad_export_dae = True
cfg.cad_export_dxf = True
cfg.cad_dxf_3dface = True
cfg.cad_use_only_measured_A = True
cfg.cad_min_face_area = 0.02
cfg.cad_collinear_tol = 0.03
cfg.cad_min_loop_vertices = 3
cfg.cad_use_c_measured_loops = True
cfg.cad_extract_c_for_horizontals = True
cfg.cad_c_cell = .05
cfg.cad_c_simplify = .15
cfg.cad_c_min_island_cells = 4
cfg.cad_c_sample = 20000
cfg.cad_c_max_raster_cells = 4_000_000
cfg.cad_c_fill_holes = False
cfg.cad_c_inner_min_area = 0.30
cfg.cad_c_close_iter = 2
cfg.cad_c_open_iter = 1
cfg.cad_reg_enabled = True
cfg.cad_reg_min_run_len = 0.15
cfg.cad_reg_run_angle_deg = 15.0
cfg.cad_reg_snap_A_dist = 0.15
cfg.cad_reg_snap_A_angle_deg = 8.0
cfg.cad_reg_merge_vertex_tol = 0.02
cfg.cad_wall_basis = True
cfg.cad_min_polygon_area = 0.05
cfg.cad_fallback_concave = True
cfg.cad_concave_ratio = 0.30
cfg.cad_consistency_sample = 5000
cfg.cad_tess_tolerance = 0.005
cfg.cad_tess_angular = 0.20

# V23 arrangement
cfg.structural_enabled = True
cfg.structural_band = 0.20
cfg.structural_extend_m = 0.20
cfg.structural_max_extend = 0.50
cfg.structural_min_len = 0.10
cfg.vertex_merge_tol = 0.03
cfg.cad_plane_res_tol = 0.01
cfg.adjacency_dist = 0.25
cfg.adjacency_min_pts = 3
cfg.adjacency_sample = 2000
cfg.structural_external_tol = 0.05
cfg.structural_external_min_wrong = 3
cfg.structural_coverage_min = 0.90
cfg.structural_obs_band = 0.15
cfg.structural_arrange_extend = 0.60
cfg.structural_cell_min_pts = 5
cfg.structural_cell_min_frac = 0.02
cfg.structural_cell_merge_tol = 0.05
cfg.horizontal_level_tol = 0.10

# V23.1 evidence & provenance
cfg.provenance_line_tol = 0.05
cfg.provenance_parallel_cos = 0.98
cfg.shared_edge_tol = 0.05
cfg.shared_edge_min_len = 0.05
cfg.unsupported_min_pts = 15
cfg.unsupported_min_density = 30.0
cfg.evidence_sample = 2000

# V23.2 recovery pre-pass
cfg.recovery_enabled = True
cfg.recovery_knn_assigned = 8
cfg.recovery_knn_unassigned = 12
cfg.recovery_dist_recoverable = 0.02
cfg.recovery_dist_extendable = 0.05
cfg.recovery_angle_recoverable_deg = 15.0
cfg.recovery_angle_extendable_deg = 25.0
cfg.recovery_uv_extend_m = 0.10
cfg.recovery_exclude_ambiguous = True
cfg.recovery_chunk = 200_000

cfg.structural_cell_min_recovered_pts = 20
cfg.structural_cell_min_recovered_frac = 0.08
cfg.structural_cell_min_core_pts_rec = 3

run(cfg) это код генерации cad моделей из облака точек но модель получается не полная с отсутствующими или не полными ступенями или в модели элементы отличаются от олака они или больше или не той формы перепиши полный исправленный код