surface_mesh_from_lines.py

# 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))})