mesh_to_poly

mesh_to_poly.py

import subprocess
import numpy as np
from scipy.sparse import lil_matrix
from matplotlib.path import Path

# ------------------------------
def read_triangle_mesh(filename):
    """
    Read in a mesh in Triangle's format.

    Parameters
    ==========
    filename: the stem of the filename for the mesh, so if the mesh is stored
              in the files <meshname>.X.node, <meshname>.X.ele, ... then
              the proper argument will be <meshname>.X.

    Returns:
    =======
    x, y: coordinates of the mesh points
    ele: (num_triangles, 3)-array of nodes in each triangle
    edge: (num_edges, 2)-array of nodes on each edge
    neigh: (num_triangles, 3)-array of neighboring triangles of each triangle
    bnd: boundary marker of each node
    edge_bnd: boundary marker of each edge
    """

    # Read in the nodes
    fid = open(filename + ".node", "r")
    nn = int(fid.readline().split()[0])
    x = np.zeros(nn, dtype = np.float64)
    y = np.zeros(nn, dtype = np.float64)
    bnd = np.zeros(nn, dtype = np.int32)
    for i in range(nn):
        line = fid.readline().split()
        x[i] = float(line[1])
        y[i] = float(line[2])
        bnd[i] = int(line[3])
    fid.close()

    # Read in the triangles
    fid = open(filename + ".ele", "r")
    nt = int(fid.readline().split()[0])
    ele = np.zeros((nt, 3), dtype = np.int32)
    for i in range(nt):
        ele[i, :] = map(lambda k: k-1,
                        map(int, fid.readline().split()[1:]))
    fid.close()

    # Read in the edges
    fid = open(filename + ".edge", "r")
    ne = int(fid.readline().split()[0])
    edge = np.zeros((ne, 2), dtype = np.int32)
    edge_bnd = np.zeros(ne, dtype = np.int32)
    for i in range(ne):
        line = map(int, fid.readline().split())
        edge[i,:] = map(lambda k:k-1, line[1:3])
        edge_bnd[i] = line[3]
    fid.close()

    # Read in the triangle neighbors
    fid = open(filename + ".ele", "r")
    neigh = np.zeros((nt, 3), dtype = np.int32)
    fid.readline()
    for i in range(nt):
        neigh[i,:] = map(lambda k: max(-1, k-1),
                         map(int, fid.readline().split()[1:]))
    fid.close()

    return x, y, ele, edge, neigh, bnd, edge_bnd


# ---------------------------------------------------------------
def mesh_to_boundary_graph(x, y, edge, bnd, edge_bnd):
    """
    Given a mesh, create a PSLG of its boundary
    """
    num_nodes = len(x)
    num_bnd_points = sum(bnd != 0)
    index = np.zeros(num_bnd_points, dtype = np.int32)
    index_inverse = np.zeros(num_nodes, dtype = np.int32)
    count = 0
    for i in range(num_nodes):
        if bnd[i] != 0:
            index[count] = i
            count += 1

    for i in range(num_bnd_points):
        index_inverse[index[i]] = i

    num_edges = np.shape(edge)[0]
    A = lil_matrix((num_bnd_points, num_bnd_points), dtype = np.int32)
    for n in range(num_edges):
        if edge_bnd[n] != 0:
            i, j = index_inverse[edge[n, :]]
            A[i, j] = edge_bnd[n]
            A[j, i] = edge_bnd[n]

    return A, index


# -------------------------
def connected_components(A):
    num_bnd_nodes = np.shape(A)[0]
    components = []
    visited = np.zeros(num_bnd_nodes, dtype = bool)
    while not all(visited):
        i = np.argmin(visited)
        visited[i] = True
        stack = [i]
        p = []
        while stack:
            i = stack.pop()
            p.append(i)
            neighbors = A[i,:].nonzero()[1]

            for j in neighbors:
                if not visited[j]:
                    visited[j] = True
                    stack.append(j)

        components.append(p)

    return components


# ------------------------
class GoofyPSLG(Exception):
    pass

def outermost_path(components, index, x, y):
    ps = map(lambda c: Path(zip(x[index[c]], y[index[c]]),
                            closed = True),
             components)

    num_components = len(ps)
    for i in range(num_components):
        for j in range(i + 1, num_components):
            if ps[i].contains_path(ps[j]):
                return i

    raise GoofyPSLG


# --------------------------------------
def holes(components, comp, index, x, y):
    xh = []
    yh = []

    for n in range(len(components)):
        if n != comp:
            p = Path(zip(x[index[components[n]]],
                         y[index[components[n]]]),
                     closed = True)

            X = 0.0
            Y = 0.0

            i = 0
            j = len(p)/2 - 1

            while not p.contains_point((X, Y)):
                j += 1
                X, Y = 0.5 * (p.vertices[i,:] + p.vertices[j,:])

            xh.append(X)
            yh.append(Y)

    return np.array(xh), np.array(yh)


# ----------------------------------------
def simplify_pslg(x, y, index, components):
    tol = np.cos(2 * np.pi / 100)
    def collinear(X, Y, Z):
        P = X - Y
        Q = Y - Z
        return np.dot(P, Q) / (np.linalg.norm(P) * np.linalg.norm(Q)) > tol

    mask = np.ones(len(index), dtype = bool)

    for component in components:
        length = len(component)

        for l in range(length):
            i = component[(l - 1) % length]
            j = component[l]
            k = component[(l + 1) % length]

            mask[j] = not collinear(np.array((x[index[i]], y[index[i]])),
                                    np.array((x[index[j]], y[index[j]])),
                                    np.array((x[index[k]], y[index[k]])))

    new_index = []

    new_components = []
    k = 0
    for component in components:
        new_component = []
        for l in range(len(component)):
            if mask[component[l]]:
                new_index.append(index[component[l]])
                new_component.append(k)
                k += 1

        new_components.append(new_component)

    return np.array(new_index), new_components


# --------------------------------------------------------------------------------
def write_to_poly(filename, x, y, index, components, outer_component, xh, yh, bnd):
    num_nodes = len(index)
    num_holes = len(xh)

    with open(filename, "w") as poly_file:
        # Write out the descriptor for the number of nodes & dimension
        poly_file.write("{0} 2 0 1\n".format(num_nodes))

        # Write out the node positions
        k = 0
        for n in range(len(components)):
            for i in index[components[n]]:
                k += 1
                poly_file.write("{0} {1} {2} {3}\n"
                                .format(k, x[i], y[i], bnd[i]))

        # Write out the descriptor for the number of edges
        poly_file.write("{0} 1\n".format(num_nodes))

        # Write out the PSLG edges
        k = 1
        for n in range(len(components)):
            component_length = len(components[n])
            for i in range(component_length):
                j = (i + 1) % component_length
                indx = index[components[n][i]]
                poly_file.write("{0} {1} {2} {3}\n"
                                .format(i + k, i + k, j + k, bnd[indx]))
            k += component_length

        # Write out the holes
        poly_file.write("{0}\n".format(num_holes))
        for n in range(num_holes):
            poly_file.write("{0} {1} {2}\n".format(n + 1, xh[n], yh[n]))

    return


# -----------------
def download_mesh():
    url = "http://students.washington.edu/shapero/meshes/ross/ross.1"
    for ext in [".node", ".ele", ".edge", ".neigh"]:
        subprocess.call(["wget", "-nc", url + ext])


# -----------------------
if __name__ == "__main__":

    # Download the mesh files
    download_mesh()

    # Read in the Triangle mesh
    x, y, ele, edge, neigh, bnd, edge_bnd = read_triangle_mesh("ross.1")

    # Construct the boundary connectivity graph
    A, index = mesh_to_boundary_graph(x, y, edge, bnd, edge_bnd)

    # Find the connected components on the boundary, e.g. each PSLG segment
    components = connected_components(A)
    index, components = simplify_pslg(x, y, index, components)

    # Find the outermost boundary segment
    comp = outermost_path(components, index, x, y)

    # Find a point inside each hole
    xh, yh = holes(components, comp, index, x, y)

    # Write out the .poly file
    write_to_poly("ross.poly", x, y, index, components, comp, xh, yh, bnd)