Skip to content

Instantly share code, notes, and snippets.

@jkfurtney
Last active May 4, 2022 20:14
Show Gist options
  • Save jkfurtney/0c95d745d577f3717501858110f9f4bb to your computer and use it in GitHub Desktop.
Save jkfurtney/0c95d745d577f3717501858110f9f4bb to your computer and use it in GitHub Desktop.
# This file was written by Jason Furtney and Wei Fu
# This file is an example of how to generate a surface mesh from lines
# in a DXF file. A vtk file is created of the mesh.
# you may have to install some Python packages to get this to work
import numpy as np
import ezdxf
import time
import math
import rdp
from meshpy.triangle import MeshInfo, build
from scipy.interpolate import LinearNDInterpolator
from scipy.spatial import cKDTree
from scipy.spatial import Delaunay, ConvexHull
from scipy.spatial.distance import pdist
from collections import defaultdict
from itertools import combinations
import meshio
# dxf manipulation functions
def filter_raw_data(b0, b1, polylines, polylines2):
def keep(p):
if b0 is None:
return True
x, y, _ = p
if x > b0[0] and x < b1[0] and y>b0[1] and y<b1[1]:
return True
return False
def add_data(polylines):
old_i = 0
for i, polyline in enumerate(polylines):
lpoints = list(polyline.points())
old_j = -1
for j, p in enumerate(lpoints):
previous_in = False
if j>0:
previous_in = keep(lpoints[j-1])
current_in = keep(p)
next_in = False
if j < len(lpoints)-1:
next_in = keep(lpoints[j+1])
contiguous = j == old_j + 1
same_polyline = i == old_i
old_j = j
old_i = i
use_point = current_in or previous_in or next_in
if use_point:
add_data.npoints += 1
add_data.points.append(p)
add_data.point_data.append(p.z)
add_segment = same_polyline and contiguous and next_in
if add_segment:
add_data.lines.append((add_data.npoints-1, add_data.npoints))
# add this point to next point
add_data.points = []
add_data.point_data = []
add_data.lines = []
add_data.npoints = 0
add_data(polylines)
add_data(polylines2)
npoints = add_data.npoints
# if b0 is not None:
# zmin = np.array(add_data.points)[:,2].min()
# add_data.points.append((b0[0], b0[1], zmin))
# add_data.points.append((b0[0], b1[1], zmin))
# add_data.points.append((b1[0], b1[1], zmin))
# add_data.points.append((b1[0], b0[1], zmin))
# add_data.lines.append((npoints+0, npoints+1))
# add_data.lines.append((npoints+1, npoints+2))
# add_data.lines.append((npoints+2, npoints+3))
# add_data.lines.append((npoints+3, npoints+0))
# for i in range(4):
# add_data.point_data.append(zmin)
return add_data.points, add_data.lines, add_data.point_data
def read_dxf_polylines(filename):
print("reading")
doc = ezdxf.readfile(filename)
print("done")
print("extracting")
msp = doc.modelspace()
#lines = msp.query('LINE[layer=="ROADS"]')
polylines = msp.query('POLYLINE[layer=="ROADS"]')
polylines2 = msp.query('POLYLINE[layer=="BREAKLINES"]')
return polylines, polylines2
def filter_dxf_polylines(b0, b1, polylines, polylines2):
points, lines, point_data = filter_raw_data(b0, b1, polylines, polylines2)
print("done")
points=np.array(points)
cells = [("line", np.array(lines))]
return points, cells, point_data
def simplify(lines, points):
glines = []
plast = -1
zhash = {}
for segment in lines:
p0, p1 = segment
if p0 == plast:
pass
else:
glines.append([])
zhash[(points[p0][0], points[p0][1])] = points[p0][2]
zhash[(points[p1][0], points[p1][1])] = points[p1][2]
glines[-1].append((points[p0][0], points[p0][1]))
glines[-1].append((points[p1][0], points[p1][1]))
plast = p1
slines=[]
for line in glines:
simplified = rdp.rdp(line, 1)
for pp in simplified:
pp.append(zhash[tuple(pp)]) # add the z component back
slines.append(simplified)
return slines
def convert(data):
point_index = 0
points = []
segments = []
for line in data:
for i, point in enumerate(line):
points.append(point)
if i>0:
segments.append((point_index-1, point_index))
point_index += 1
return np.array(points), segments
def mesh_line_segments(lines, points):
print(1)
point_index = 0
reference_point = np.copy(points[0])
points -= reference_point
data = simplify(lines, points)
points, segments = convert(data)
#anp.save("points.npy", points)
#np.save("lines.npy", segments)
zhash = {}
for p in points:
zhash[(p[0], p[1])] = p[2]
print(2)
print("refining line segments")
threshold = 10.0
# add points to line segments to make them at most 10 ft apart
new_points = []
for i, j in segments:
p0, p1 = points[i], points[j]
d = np.linalg.norm(p0-p1)
if d > threshold:
n = int(d/threshold)+1
newx = np.linspace(p0[0], p1[0], n+1, endpoint=False)[1:]
newy = np.linspace(p0[1], p1[1], n+1, endpoint=False)[1:]
newz = np.linspace(p0[2], p1[2], n+1, endpoint=False)[1:]
new_points += np.array((newx, newy, newz)).T.tolist()
points = np.vstack((points, new_points))
print("building height interpolation mesh")
mesh = Delaunay(points[:,:2])
print(4)
tree = cKDTree(points[:,:2])
X, Y = np.array(mesh.points).T
# look up heights of any new points that were created, just in case
_, indices = tree.query(mesh.points)
Z0=[]
for i, _ in zip(indices, mesh.points):
Z0.append(points[:,2][i])
# create a height interpolation object to use later
Z_interp = LinearNDInterpolator(mesh.points, Z0)
# now we make the refined mesh in 2D
print("building refined mesh")
mesh_info = MeshInfo()
mesh_info.set_points(mesh.points)
hull = ConvexHull(mesh.points)
facets = hull.simplices
mesh_info.set_facets(facets)
new_mesh = build(mesh_info,
allow_boundary_steiner=False,
max_volume=300)
X1, Y1 = np.array(new_mesh.points).T
new_z = Z_interp(X1, Y1)
new_points = list(new_mesh.points)
final_points = np.array((*np.array(new_mesh.points).T, new_z.data)).T
# now the new refined mesh has the correct height at the new points
return final_points+reference_point, new_mesh.elements
if __name__=="__main__":
dxf_path = "my_dxf_file.dxf"
poly0, poly1 = read_dxf_polylines(dxf_path)
# filter anything outside the rectangle defined by these points
a0 = 0, 0
a1 = 20_000, 20_000
points, cells, point_data = filter_dxf_polylines(a0, a1, poly0, poly1)
meshio.write_points_cells('data0.vtu', # this is just the poly lines
points=points,
cells=cells,
point_data={"data": np.array(point_data)})
points, cells = mesh_line_segments(lines, points)
meshio.write_points_cells("data1.vtu", # this is the surface mesh
points=points,
cells=[("triangle", cells)],
point_data={"data": np.zeros(len(points))})
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment