113 lines
6.5 KiB
Python
113 lines
6.5 KiB
Python
"""Evidence-based rigid registration of original plan pier geometry.
|
|
|
|
No imagery, text insertion anchors, or user offsets are accepted as controls.
|
|
All pier centers, axes and hold-out geometry must agree; ambiguous fits fail.
|
|
"""
|
|
import itertools,math
|
|
import numpy as np
|
|
|
|
def entities(layout,depth=0):
|
|
for e in layout:
|
|
if e.dxftype()=='INSERT' and depth<6:
|
|
try:yield from entities(e.virtual_entities(),depth+1)
|
|
except Exception:continue
|
|
else:yield e
|
|
|
|
def rectangles(doc,lengths):
|
|
found=[]
|
|
for e in entities(doc.modelspace()):
|
|
if e.dxftype()=='LWPOLYLINE':
|
|
if not e.closed or any(abs(p[4])>1e-10 for p in e.get_points()):continue
|
|
a=np.array([tuple(p) for p in e.vertices_in_wcs()])
|
|
elif e.dxftype()=='POLYLINE':
|
|
if not e.is_closed:continue
|
|
a=np.array([tuple(p) for p in e.points_in_wcs()])
|
|
else:continue
|
|
if len(a)==5 and np.linalg.norm(a[0]-a[-1])<1e-6:a=a[:-1]
|
|
if len(a)!=4 or np.ptp(a[:,2])>1e-6:continue
|
|
a=a[:,:2];edges=np.roll(a,-1,axis=0)-a;ls=np.linalg.norm(edges,axis=1)
|
|
if ls.min()<.01:continue
|
|
# A skew pier is drawn as a parallelogram, not an orthogonal rectangle.
|
|
# Opposite edges must still close and balance precisely.
|
|
if np.linalg.norm(edges[0]+edges[2])>1e-4 or np.linalg.norm(edges[1]+edges[3])>1e-4:continue
|
|
if max(abs(np.dot(edges[i],edges[(i+1)%4]))/(ls[i]*ls[(i+1)%4]) for i in range(4))>.5:continue
|
|
if ls.max()/ls.min()<2.5:continue
|
|
if min(abs(ls.max()-v) for v in lengths)>3:continue
|
|
axis=edges[np.argmax(ls)]/ls.max();center=a.mean(0)
|
|
if any(np.linalg.norm(center-r['center'])<1e-5 for r in found):continue
|
|
found.append(dict(center=center.tolist(),axis=axis.tolist(),length=float(ls.max()),width=float(ls.min()),
|
|
handle=e.dxf.handle,layer=e.dxf.layer,vertices=a.tolist()))
|
|
# Many original drawings store cap rectangles as four separate LINEs.
|
|
# Require a closed parallelogram cycle; a text label or lone centerline is
|
|
# insufficient to establish a physical center.
|
|
lines=[]
|
|
for e in entities(doc.modelspace()):
|
|
if e.dxftype()!='LINE' or not any(t in e.dxf.layer.upper() for t in ('CONC','BRDG','교량')):continue
|
|
a=np.array(tuple(e.dxf.start));b=np.array(tuple(e.dxf.end))
|
|
if abs(a[2]-b[2])>1e-6:continue
|
|
length=np.linalg.norm(a[:2]-b[:2])
|
|
if .1<length<max(lengths)+4:lines.append((a[:2],b[:2],e))
|
|
index={}
|
|
key=lambda p:tuple(np.round(p,4))
|
|
for i,(a,b,e) in enumerate(lines):
|
|
index.setdefault(key(a),[]).append(i);index.setdefault(key(b),[]).append(i)
|
|
for i,(a,b,e) in enumerate(lines):
|
|
length=np.linalg.norm(b-a)
|
|
if min(abs(length-v) for v in lengths)>3:continue
|
|
for j in index.get(key(a),[]):
|
|
if i==j:continue
|
|
c,d,_=lines[j];c=d if np.linalg.norm(c-a)<1e-4 else c
|
|
width=np.linalg.norm(c-a)
|
|
if width<.1 or length/width<2.5 or abs(np.dot(c-a,b-a))/(width*length)>.5:continue
|
|
expected=b+c-a
|
|
for k in index.get(key(b),[]):
|
|
if k==i:continue
|
|
p,q,_=lines[k];d=q if np.linalg.norm(p-b)<1e-4 else p
|
|
if np.linalg.norm(d-expected)>1e-4:continue
|
|
closed=any(min(np.linalg.norm(lines[m][0]-d),np.linalg.norm(lines[m][1]-d))<1e-4 for m in index.get(key(c),[]))
|
|
if not closed:continue
|
|
center=(a+b+c+d)/4
|
|
if any(np.linalg.norm(center-r['center'])<1e-4 for r in found):continue
|
|
found.append(dict(center=center.tolist(),axis=((b-a)/length).tolist(),length=float(length),width=float(width),
|
|
handle=e.dxf.handle,layer=e.dxf.layer,vertices=[p.tolist() for p in (a,b,d,c)],
|
|
evidence_kind='closed_four_LINE_rectangle'))
|
|
return found
|
|
|
|
def fit(source,target):
|
|
source=np.array(source);target=np.array(target);s=source-source.mean(0);t=target-target.mean(0)
|
|
u,_,vt=np.linalg.svd(s.T@t);rotation=u@vt
|
|
if np.linalg.det(rotation)<0:return None
|
|
translation=target.mean(0)-source.mean(0)@rotation
|
|
residual=np.linalg.norm(source@rotation+translation-target,axis=1)
|
|
return rotation,translation,residual
|
|
|
|
def resolve_plan_registration(doc, controls):
|
|
axes=[np.array(c['axis'])[:,::-1] for c in controls]
|
|
target=np.array([a.mean(0) for a in axes]);directions=np.array([(a[1]-a[0])/np.linalg.norm(a[1]-a[0]) for a in axes])
|
|
candidates=rectangles(doc,[np.linalg.norm(a[1]-a[0]) for a in axes]);src=np.array([r['center'] for r in candidates])
|
|
solutions=[]
|
|
if len(target)>=3 and len(src)>=len(target):
|
|
target_delta=target[-1]-target[0];distance=np.linalg.norm(target_delta)
|
|
for i,j in itertools.permutations(range(len(src)),2):
|
|
delta=src[j]-src[i];length=np.linalg.norm(delta)
|
|
if abs(length-distance)>.05:continue
|
|
angle=math.atan2(target_delta[1],target_delta[0])-math.atan2(delta[1],delta[0])
|
|
c,s=math.cos(angle),math.sin(angle);rotation=np.array([[c,s],[-s,c]])
|
|
translation=target[0]-src[i]@rotation
|
|
distances=np.linalg.norm((src@rotation+translation)[:,None,:]-target[None,:,:],axis=2)
|
|
matched=np.argmin(distances,axis=0)
|
|
if len(set(matched))!=len(target) or np.max(distances[matched,np.arange(len(target))])>.1:continue
|
|
fitted=fit(src[matched],target)
|
|
if fitted is None:continue
|
|
rotation,translation,residual=fitted
|
|
if residual.max()>.03:continue
|
|
dots=[abs(np.dot(np.array(candidates[idx]['axis'])@rotation,d)) for idx,d in zip(matched,directions)]
|
|
if min(dots)<.999:continue
|
|
if any(np.max(np.abs(rotation-np.array(v['rotation'])))<1e-7 and np.linalg.norm(translation-v['translation'])<.01 for v in solutions):continue
|
|
solutions.append(dict(rotation=rotation.tolist(),translation=translation.tolist(),residual_m=residual.tolist(),
|
|
controls=[dict(pier=c['id'],**candidates[idx],target_EN=target[k].tolist()) for k,(idx,c) in enumerate(zip(matched,controls))]))
|
|
evidence=dict(candidate_rectangles=candidates,solutions=solutions,status='registered' if len(solutions)==1 else 'ambiguous' if solutions else 'no_verified_frame',
|
|
method='rigid_registration_of_original_plan_pier_rectangle_centers_to_Excel_IP1_IP3_midpoints',
|
|
max_permitted_residual_m=.03)
|
|
return evidence
|