Search code examples
pythontriangulationdelaunay

How to generate edge index after Delaunay triangulation?


I'm using Python 3.7. There is a set of points. I generate Delaunay triangulation through these points.

import numpy as np
points = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1], [1.5, 0.6], [1.2, 0.5], [1.7, 0.9], [1.1, 0.1]])
from scipy.spatial import Delaunay
tri = Delaunay(points)

How can I remove some edges through the edge length threshold? How to plot the new Delaunay triangulation (after remove the edges)? My idea is generating edge index through the point. Like,

[e1, (0,0),(1,0),e1_length], [e2, (0,0),(1,1),e2_length], ...

Solution

  • We need to make three operations: convert triangles from Delaunay object to the set of edges (with removing duplicates), calculate length for each edge and select edges which meets the criterion.

    Creating set of edges and calculation of lengths:

    def less_first(a, b):
        return [a,b] if a < b else [b,a]
    
    def delaunay2edges(tri):
    
        list_of_edges = []
    
        for triangle in tri.simplices:
            for e1, e2 in [[0,1],[1,2],[2,0]]: # for all edges of triangle
                list_of_edges.append(less_first(triangle[e1],triangle[e2])) # always lesser index first
    
        array_of_edges = np.unique(list_of_edges, axis=0) # remove duplicates
    
        list_of_lengths = []
    
        for p1,p2 in array_of_edges:
            x1, y1 = tri.points[p1]
            x2, y2 = tri.points[p2]
            list_of_lengths.append((x1-x2)**2 + (y1-y2)**2)
    
        array_of_lengths = np.sqrt(np.array(list_of_lengths))
    
        return array_of_edges, array_of_lengths
    
    edges, lengths = delaunay2edges(tri)
    

    Selecting edges by criterion (length > 0.5 for example):

    criterion = np.argwhere(lengths > 0.5).flatten()
    
    selected_edges = edges[criterion]
    
    print('Removed', len(edges) - len(selected_edges), 'edges')
    

    Plotting:

    import matplotlib.pyplot as plt
        
    plt.triplot(tri.points[:,0], tri.points[:,1], tri.simplices, color='red')
    
    x_lines = []
    y_lines = []
    
    for p1,p2 in selected_edges:
        x1,y1 = points[p1]
        x2,y2 = points[p2]
        plt.plot([x1,x2],[y1,y2], color='blue')
    
    plt.scatter(points[:,0],points[:,1])
    
    plt.show()