-- 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 моделей из облака точек но модель получается не полная с отсутствующими или не полными ступенями или в модели элементы отличаются от олака они или больше или не той формы перепиши полный исправленный код

Скачать frugal_cad_v24_robust.py где код или ссылка

Проверил V24 ROBUST и сопоставил его с твоим FRUGAL_CAD V23 Structural Arrangement, результатами аудитов и подходами Poux.

Мой вывод: запускать V24 ROBUST как замену V23 пока не рекомендую. В нём есть полезные идеи, особенно поиск новых плоскостей среди segment_id=0, но обнаружены серьёзные алгоритмические ошибки. Кроме того, код возвращается к подходу, от которого ты уже отказался из-за рваных границ.

Разберу конкретно.

1. Главная архитектурная проблема: V24 отказывается от Structural Arrangement

В V23 геометрия формируется следующим образом:

text
Physical planes
      ↓
Structural intersections
      ↓
Planar arrangement
      ↓
Supported CAD cells
      ↓
CAD faces

В V24 ROBUST:

text
Physical planes
      ↓
Raster footprint
      ↓
Marching squares
      ↓
Douglas-Peucker
      ↓
Soft snapping
      ↓
CAD faces

Это фактически возврат к V18, но с дополнительным восстановлением точек.

Ключевой участок:

python
def build_surface_loops(S, cfg):
    ...
    loops = extract_surface_footprints(s, cfg)

То есть структурные пересечения больше не определяют форму CAD-граней.

Они используются только здесь:

python
uv = snap_loop_uv(
    q["loop_uv"],
    inter.get(s["id"], []),
    cfg
)

И только при близости до 45 мм.

Это противоречит твоему основному требованию: получить чёткую конструктивную геометрию с согласованными стенами, потолками, колоннами и лестницами.

V24 может сохранить больше видимых фрагментов, но с большой вероятностью снова создаст обрезанные поверхности и пробелы на стыках.

2. Критическая ошибка: recovery выполняется дважды

В run():

python
used_recovery = recover_unassigned_to_existing(
    P, sid, S, cfg
)

После этого:

python
S, used_discovery = discover_missing_planes(
    P, sid, used_recovery, S, cfg
)

Затем:

python
used_recovery2 = recover_unassigned_to_existing(
    P, sid, S, cfg
)

Но внутри recover_unassigned_to_existing():

python
for s in S:
    s["P_recovered"] = np.empty((0, 3), float)
    s["P_recovered_ambig"] = np.empty((0, 3), float)

Второй вызов полностью уничтожает результаты первого recovery.

Более того, discover_missing_planes() может присоединить найденные точки к существующим поверхностям:

python
S[i]["P_recovered"] = np.vstack(...)

Эти точки тоже стираются вторым вызовом.

Это не просто неэффективность.

Это логическая ошибка, из-за которой часть найденной геометрии теряется до построения CAD.

Исправление

Recovery должен выполняться один раз.

После discovery нужно сохранить уже найденные точки и отдельно обработать новые поверхности, не сбрасывая P_recovered.

3. Ещё одна серьёзная ошибка: найденные плоскости могут использовать одни и те же точки

В discover_missing_planes():

python
candidates.sort(
    key=lambda s: len(s["P"]),
    reverse=True
)

Затем:

python
if float(np.mean(used_local[loc])) > 0.60:
    continue

То есть новый кандидат отклоняется только тогда, когда более 60% его точек уже использованы.

Следовательно, два принятых кандидата могут иметь значительное пересечение по исходным точкам.

Например:

text
Candidate A: 1000 points
Candidate B: 800 points

Overlap: 400 points

Оба кандидата будут приняты.

Тогда одна и та же геометрия может породить две разные CAD-плоскости.

Для лестниц и подиумов это особенно опасно.

Каждая восстановленная точка должна иметь единственного владельца либо явно храниться как неоднозначное evidence без включения в геометрию нескольких физических поверхностей.

4. Discovery использует нестабильное квантование плоскостей

Сейчас:

python
def quantized_plane_key(p, n, cfg):

    qn = tuple(
        np.round(
            n / cfg.discover_normal_quant
        ).astype(np.int16)
    )

    d = -float(n @ p)

    qd = int(
        round(d / cfg.discover_d_quant)
    )

    return qn + (qd,)

Параметры:

python
discover_normal_quant = 0.16
discover_d_quant = 0.055

Проблема в том, что даже точки одной физической плоскости могут попасть в разные группы из-за небольших изменений локальных нормалей.

И наоборот, близкие, но различные плоскости могут попасть в одну группу.

Особенно опасна независимая компонентная квантизация нормали: она не эквивалентна ограничению угла между плоскостями.

Это допустимо как способ генерации кандидатов, но не как окончательное решение о принадлежности точек физической плоскости.

Что лучше

Использовать quantized keys только для быстрого поиска кандидатов.

После этого обязательно проверять:

• угловую согласованность нормалей;
• расстояния до общей плоскости;
• пространственную связность;
• RMS и p95;
• распределение точек по площади.

5. Ошибка в восстановлении исходных точек обнаруженной плоскости

В discover_missing_planes():

python
res = np.abs(Prem @ s["n"] + s["d"])

near = res <= cfg.discover_rms_max * 1.5

Здесь проверяется расстояние до бесконечной плоскости.

Затем применяется UV bounding box.

Но UV bounding box может содержать пустые области между разными конструктивными элементами.

Это позволяет кандидату захватывать точки другой поверхности, если они лежат на близкой плоскости.

Особенно при ступенчатых конструкциях.

Для discovery необходима локальная пространственная связность после присоединения исходных точек, а не только до него.

6. Ошибка в обработке отверстий

В extract_surface_footprints():

python
# suppress contours that are holes

Затем:

python
if p2.contains(rp):
    inside_other = True

И вложенные контуры удаляются.

Для CAD это неправильно.

Если у стены есть оконный или дверной проём, внутренний контур должен стать отверстием в CAD face.

А здесь он просто исчезает.

В результате:

text
Wall with opening
       ↓
Outer polygon only
       ↓
Solid wall face

То есть алгоритм может закрывать реальные проёмы.

Нужно хранить:

python
outer_loop
inner_loops

И передавать их в:

python
cq.Face.makeFromWires(
    outer_wire,
    inner_wires
)

Но даже это требует проверки: отверстие в растровой маске не всегда является настоящим строительным проёмом.

7. Shared vertices по-прежнему не обеспечивают общие рёбра

В coordinate_shared_vertices():

python
mean = P[ids].mean(0)

Затем:

python
p = project_to_plane(p, s)

Это та же проблема, которую мы уже обнаружили в V22/V23.

После независимой проекции на разные плоскости ранее объединённые вершины снова могут разойтись.

Поэтому V24 не гарантирует отсутствия щелей между CAD-гранями.

Более того, алгоритм объединяет вершины исключительно по расстоянию.

Он не проверяет, действительно ли они принадлежат одной конструктивной вершине.

8. В V24 есть проблема с оценкой геометрической точности

Функция:

python
def consistency_to_cloud(S, cfg):

Считает расстояние от точек поверхности до ближайшего ребра CAD-полигона.

Но для большой плоской поверхности точки внутри грани закономерно могут находиться далеко от её периметра.

Например, точка в центре потолка размером 5 × 5 м будет находиться на расстоянии около 2,5 м от ближайшего края.

Это не ошибка CAD.

Следовательно, такая метрика не измеряет качество приближения CAD-грани облаком.

Она будет наказывать большие правильные поверхности.

Нужна другая метрика

Для каждой исходной точки:

1. Вычислить расстояние до плоскости CAD face.
2. Спроецировать точку в UV.
3. Проверить попадание внутрь полигона.
4. Если точка вне полигона — вычислить расстояние до его границы.

Только такая комбинация позволит отделить ошибку положения плоскости от неполного покрытия.

9. Что в V24 действительно полезно

Несмотря на проблемы, я бы сохранил несколько компонентов.

Первый — discovery missing planes.

Это наиболее интересное новшество.

По аудиту:

text
segment_id=0: 1 047 040 points

Среди них потенциально находятся пропущенные физические поверхности.

Поэтому отдельный поиск новых planar patches оправдан.

Но его нужно использовать как дополнительный этап перед V23, а не как замену structural arrangement.

Второй — разделение P и P_recovered.

Это правильная архитектура:

python
"P": original_segment_points

"P_recovered": recovered_support_points

Плоскость остаётся стабильной.

Восстановленные точки расширяют evidence.

Третий — консервативный merge.

Новые параметры:

python
merge_angle_deg = 3.0
merge_plane_dist = 0.035
merge_spatial_gap = 0.20

Могут помочь сохранить небольшие конструктивные элементы.

Однако менять их следует отдельным экспериментом, потому что они также способны увеличить фрагментацию стен.

Четвёртый — минимальное перемещение при snapping.

python
snap_endpoint_max_move = 0.055

Контроль максимального перемещения полезен.

Но snapping должен работать с конструктивными рёбрами и общей топологией, а не только с измеренными растровыми контурами.

10. Что я бы сделал вместо V24 ROBUST

Я бы взял V23 Structural Arrangement как основную версию и перенёс в неё только модуль поиска пропущенных плоскостей.

Архитектура:

text
FRUGAL SEGMENTED CLOUD
          │
          ▼
ORIGINAL PHYSICAL SURFACES
          │
          ├──────────────┐
          │              │
          ▼              ▼
EXISTING PLANES     SID=0 POINTS
          │              │
          │         PLANAR DISCOVERY
          │              │
          │              ▼
          │         NEW PLANE PATCHES
          │              │
          └──────┬───────┘
                 │
                 ▼
       CONSERVATIVE MERGE
                 │
                 ▼
       STRUCTURAL ADJACENCY
                 │
                 ▼
       PLANE INTERSECTIONS
                 │
                 ▼
       STRUCTURAL ARRANGEMENT
                 │
                 ▼
       CAD CELL EVIDENCE
                 │
                 ▼
       VALIDATED CAD FACES
                 │
                 ▼
          DXF / PLY / STEP

При этом новые плоскости должны быть явно помечены как discovered_sid0.

Их нужно сравнивать с исходными физическими плоскостями, чтобы не создавать дубликаты.

11. Есть ещё один момент, который я бы исправил

В V24:

python
cfg.edge_las = os.path.join(
    base,
    "roughness_edges8.las"
)

Но edge_las нигде фактически не загружается и не используется.

Следовательно, Source B исчез из pipeline полностью.

Это не обязательно ошибка, поскольку у тебя он и раньше был отключён.

Но конфигурация не должна создавать впечатление, что roughness edges влияют на результат.

Итоговая оценка

V24 ROBUST — полезный эксперимент по восстановлению недостающих плоскостей, но геометрически это шаг назад относительно V23 Structural Arrangement.

В нём есть три особенно серьёзных риска:

1. Повторный recovery уничтожает часть ранее восстановленных точек.
2. Discovery допускает перекрывающиеся кандидаты и ошибочное присвоение точек.
3. Основной CAD-контур снова определяется растровой границей наблюдений.

Я бы не пытался исправить V24 ROBUST целиком.

Рациональнее сохранить V23 и разработать отдельный модуль discover_missing_planes(), который добавляет подтверждённые новые физические плоскости, не изменяя существующие.

После этого уже имеющийся structural arrangement сможет использовать новые поверхности для восстановления ступеней, горизонтальных переходов и других деталей.

Так мы сохраним достигнутую чёткость CAD-геометрии и расширим её полноту, не возвращаясь к рваным краям.

%runfile C:/Users/dimma/SAMO/CODE/untitled0.py --wdir
2026-10-09 17:56:44,135 | F2CAD-V24 | INFO | Loaded points=2904376, assigned=1857336, unassigned=1047040
2026-10-09 17:56:47,782 | F2CAD-V24 | INFO | Initial segmented surfaces: 1027
2026-10-09 17:56:47,872 | F2CAD-V24 | INFO | Merge pass 1: 1027 -> 933
2026-10-09 17:56:47,923 | F2CAD-V24 | INFO | Merge pass 2: 933 -> 933
2026-10-09 17:57:14,488 | F2CAD-V24 | INFO | Recovery to existing surfaces: clean=456133 ambiguous=39460
2026-10-09 17:57:18,510 | F2CAD-V24 | WARNING | Discovery produced 31367 quantized groups; processing largest 4000
2026-10-09 17:57:18,690 | F2CAD-V24 | INFO | Missing-plane discovery: added=5 attached=2 used_points=1782
2026-10-09 17:57:45,284 | F2CAD-V24 | INFO | Recovery to existing surfaces: clean=456133 ambiguous=39460
2026-10-09 17:58:00,323 | F2CAD-V24 | INFO | Footprint loops: cells=955, surfaces without loop=2
2026-10-09 17:58:02,478 | F2CAD-V24 | INFO | CAD faces built=955 rejected=0
2026-10-09 17:58:04,628 | F2CAD-V24 | INFO | creating ACAD_COLOR dictionary
2026-10-09 17:58:04,629 | F2CAD-V24 | INFO | creating ACAD_GROUP dictionary
2026-10-09 17:58:04,629 | F2CAD-V24 | INFO | creating ACAD_LAYOUT dictionary
2026-10-09 17:58:04,630 | F2CAD-V24 | INFO | creating ACAD_MATERIAL dictionary
2026-10-09 17:58:04,630 | F2CAD-V24 | INFO | creating ACAD_MLEADERSTYLE dictionary
2026-10-09 17:58:04,631 | F2CAD-V24 | INFO | creating ACAD_MLINESTYLE dictionary
2026-10-09 17:58:04,631 | F2CAD-V24 | INFO | creating ACAD_PLOTSETTINGS dictionary
2026-10-09 17:58:04,632 | F2CAD-V24 | INFO | creating ACAD_PLOTSTYLENAME dictionary
2026-10-09 17:58:04,632 | F2CAD-V24 | INFO | creating ACAD_SCALELIST dictionary
2026-10-09 17:58:04,633 | F2CAD-V24 | INFO | creating ACAD_TABLESTYLE dictionary
2026-10-09 17:58:04,633 | F2CAD-V24 | INFO | creating ACAD_VISUALSTYLE dictionary
2026-10-09 17:58:06,837 | F2CAD-V24 | INFO | did not write header var $INTERFEREOBJVS, value is None. 2026-10-09 17:58:06,839 | F2CAD-V24 | INFO | did not write header var$INTERFEREVPVS, value is None. твой предыдущий скрипт показал хорошие результаты показал но он завис очень долго обрабатывает и не генерирует skeleton dxf и другие выходные данные сделай оптимизацию и улучшениякода и перепиши полный исправленный код