diff --git a/pierzainer_ground/pierzainer_ground/survey_terrain.py b/pierzainer_ground/pierzainer_ground/survey_terrain.py index 66d6595..26eec2e 100644 --- a/pierzainer_ground/pierzainer_ground/survey_terrain.py +++ b/pierzainer_ground/pierzainer_ground/survey_terrain.py @@ -17,6 +17,45 @@ from pierzainer_ground.dtm import DTM ROOT=Path(__file__).resolve().parent.parent +def reject_height_outliers(entities, radius=15.0, max_dz=10.0, min_levels=3, same=1e-3): + """Drop height entities that disagree with the surrounding surveyed surface. + + Drawings often keep plan symbols, hatching or design outlines at a stray + constant elevation (e.g. EL 5.624 or -112.719 among 20-40 m contours). + Each entity's height is compared with the median of the *distinct* nearby + height levels (each level votes once, so a dense stray group of short + ticks cannot outvote real contours); the worst height level is removed + first and the test is repeated, so a large stray group cannot mask itself. + Sparse areas fall back to a 3x search radius. + Heights within `same` (1 mm) form one level, absorbing CAD float noise. + Returns (keep flags, [(index, z, deviation)]). + """ + from scipy.spatial import cKDTree + n=len(entities) + keep=np.ones(n,bool) + if n<3:return keep.tolist(),[] + xy=np.vstack([e[0] for e in entities]);owner=np.concatenate([np.full(len(e[0]),i) for i,e in enumerate(entities)]) + zs=np.array([float(e[1]) for e in entities]);sample_z=zs[owner] + tree=cKDTree(xy);rejected=[] + near=lambda e,r:np.unique(np.concatenate(tree.query_ball_point(e[0],r)) if len(e[0]) else np.array([],int)).astype(int) + neighbours=[near(e,radius) for e in entities];wide=[near(e,3*radius) for e in entities] + for _ in range(64): + alive=keep[owner];dev=np.zeros(n) + for i in np.flatnonzero(keep): + nb=neighbours[i];nb=nb[alive[nb]&(np.abs(sample_z[nb]-zs[i])>same)] + if len(np.unique(np.round(sample_z[nb]/same)))same)] + levels=np.unique(np.round(sample_z[nb]/same))*same + # Peaks/pits legitimately exceed the median; require the height to + # lie clearly outside the whole local range as well. + outside=max(levels[0]-zs[i],zs[i]-levels[-1]) if len(levels) else 0.0 + if len(levels)>=min_levels and outside>max_dz/2:dev[i]=zs[i]-float(np.median(levels)) + worst=int(np.argmax(np.abs(dev))) + if abs(dev[worst])<=max_dz:break + group=np.flatnonzero(keep&(np.abs(zs-zs[worst])<=same)&(np.abs(dev)>max_dz)) + keep[group]=False;rejected.extend((int(i),float(zs[i]),float(dev[i])) for i in group) + return keep.tolist(),rejected + + def read_survey_block(source, project, doc=None): """Read actual height geometry in an explicitly selected survey block.""" if doc is None: @@ -54,9 +93,11 @@ def read_survey_block(source, project, doc=None): if e.is_closed: xy.append(xy[0]) contours.append((xy,float(a[0,2]),e.dxf.handle)); types['constant_height_polylines']+=1 # POINT entities have priority over densified contour vertices. - for xy,z,handle in contours: - samples=_densify(xy,project['contour_step_m']) - points.extend((x,y,z) for x,y in samples) + sampled=[(np.asarray(_densify(xy,project['contour_step_m']),dtype=float).reshape(-1,2),z) for xy,z,handle in contours] + sampled=[(np.array([[x,y]]),z) for x,y,z in points]+sampled + keep,rejected=reject_height_outliers(sampled,project.get('outlier_radius_m',15.0),project.get('outlier_dz_m',10.0)) if project.get('reject_outliers',True) else ([True]*len(sampled),[]) + if rejected:types['rejected_height_outliers']=len(rejected) + points=[(x,y,z) for (xy,z),ok in zip(sampled,keep) if ok for x,y in xy] a=np.array(points,dtype=float) if len(a)<3: raise ValueError('선택한 블록에서 실제 표고점·등고선을 3점 이상 읽지 못했습니다.') diff --git a/pierzainer_ground/tests/test_survey_polyline.py b/pierzainer_ground/tests/test_survey_polyline.py index eff83ca..cd9e740 100644 --- a/pierzainer_ground/tests/test_survey_polyline.py +++ b/pierzainer_ground/tests/test_survey_polyline.py @@ -25,4 +25,23 @@ class SurveyPolylineTests(unittest.TestCase): doc=ezdxf.new();doc.modelspace().add_polyline3d([(1,1,20),(2,2,30),(3,3,40)]) with self.assertRaises(ValueError):read_survey_block(Path('unused.dxf'),self.options(),doc=doc) + def contours(self,msp,levels): + for k,z in enumerate(levels): + msp.add_lwpolyline([(0,10*k),(200,10*k)],dxfattribs={'elevation':z}) + + def test_stray_constant_height_symbols_are_rejected(self): + # Ungok case: dense short ticks at EL 5.624 and a few at -112.719 among 20-40 m contours. + doc=ezdxf.new();msp=doc.modelspace();self.contours(msp,[28,29,30,31,32,33]) + for i in range(60):msp.add_lwpolyline([(40+2*i,20),(40+2*i,22)],dxfattribs={'elevation':5.624}) + for i in range(4):msp.add_lwpolyline([(100+i,25),(100+i,27)],dxfattribs={'elevation':-112.719}) + points,counts,_=read_survey_block(Path('unused.dxf'),self.options(),doc=doc) + self.assertEqual(counts['rejected_height_outliers'],64) + self.assertGreaterEqual(points[:,2].min(),28) + + def test_real_peak_above_contours_is_kept(self): + doc=ezdxf.new();msp=doc.modelspace();self.contours(msp,[70,72,74,76,78,80]) + msp.add_lwpolyline([(100,25),(100,27)],dxfattribs={'elevation':84}) + points,counts,_=read_survey_block(Path('unused.dxf'),self.options(),doc=doc) + self.assertNotIn('rejected_height_outliers',counts);self.assertEqual(points[:,2].max(),84) + if __name__=='__main__':unittest.main()