How to change color of Delaunay points in Python?

Viewed 261

I am fairly new to Python, and I have been tasked to make a script that creates a plot of a Delaunay triangulation. I have successfully done so, but I need the colors of the points to change based on how many neighbors the plots have. I have a function that defines the number of neighbors and puts all points into a dictionary with the number of neighbors as the key. However, I don't know how to apply that dictionary to my plot - is there a specific function where you can take the plot created with triplot and change the colors based on an outside dictionary? Here is my code so far:

from curved_analysis import read_xyz, read_nfo, num_neighbors
from scipy.spatial import Delaunay
import matplotlib.tri as mtri
import numpy as np
import matplotlib.pyplot as plt


coords = np.array(read_xyz("traj0.xyz"))
for k in range(coords.shape[0]):
    points = coords[k]
    tri = Delaunay(points[:, :2], qhull_options=('Qz'))




neigh = num_neighbors(tri)
for key in neigh:
    if key <=5:
        plt.triplot(neigh[key], color = 'green')
    if key == 6:
        plt.triplot(neigh[key], color = 'red')
    if key >= 7:
        plt.triplot(neigh[key], color = 'yellow')



plt.triplot(points[:,0], points[:,1], tri.simplices)
plt.plot(points[:,0], points[:,1], 'o')
plt.show()

Thank you for your help!

2 Answers

You want the coloring to be done on the vertices (i.e., the points) of the triangulation, so you should do it with the plt.plot function (replacing the line before last).

Here is one way to do it (continuing your code). Replace the line before last with the following lines (the for key loop above it should also be removed).

for k, count in neigh.items():
    if count <= 5:
        plt.plot(tri.points[k, 0], tri.points[k, 1], "og")
    elif count == 6:
        plt.plot(tri.points[k, 0], tri.points[k, 1], "or")
    else:
        assert count >= 7
        plt.plot(tri.points[k, 0], tri.points[k, 1], "oy")

On a triangulation from a random sample of 200 points I get the following figure:

enter image description here

Note: for computing the neigh dictionary I used the following function (you didn't give your implementation). The key of the dictionary is the index of the vertex in the triangulation points. Note also that there are more efficient ways to do it. For example, storing the count in a numpy array instead of a dictionary (e.g., using np.diff(indptr)) and using vectorized operations for each case instead of plotting one by one in a loop.

def num_neighbors(tri):
    indptr, indices = tri.vertex_neighbor_vertices
    assert len(indptr) == len(tri.points) + 1
    vertex_order_map = {}
    for k in range(len(indptr) - 1):
        num_neighbors_k = indptr[k+1] - indptr[k]
        vertex_order_map[k] = num_neighbors_k
    return vertex_order_map 

EDIT: Below is the more efficient vectorized code. The results are the same (up to the order of coloring of the points).

indptr, indices = tri.vertex_neighbor_vertices
neighbor_counter_array = np.diff(indptr)
cond_le_5 = (neighbor_counter_array <= 5)
cond_eq_6 = (neighbor_counter_array == 6)
cond_ge_7 = (neighbor_counter_array >= 7)
plt.plot(tri.points[cond_le_5, 0], tri.points[cond_le_5, 1], "og")
plt.plot(tri.points[cond_eq_6, 0], tri.points[cond_eq_6, 1], "or")
plt.plot(tri.points[cond_ge_7, 0], tri.points[cond_ge_7, 1], "oy")

Numpy allows to filter points via a condition. For example points[num_neigh <= 5] will collect all the points with 5 or less neighbors, provided num_neigh is a numpy array that contains the number of neighbors for each point.

num_neigh can be computed as the difference between successive entries in tri.vertex_neighbor_vertices[0], so num_neigh = np.diff(tri.vertex_neighbor_vertices[0]).

The following code uses plt.scatter to draw the points with some larger dot. The code supposes the first dimension in coords represents some kind of layer. Note that in Python it is recommended to write loops directly with the elements instead of via indices. The question is not so clear whether each layer should be a separate plot, or whether they should be combined in subplots.

import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
import numpy as np

coords = np.random.uniform(0, 10, size=(1, 40, 2))
for points in coords:
    tri = Delaunay(points[:, :2], qhull_options=('Qz'))
    num_neigh = np.diff(tri.vertex_neighbor_vertices[0])

    plt.triplot(points[:, 0], points[:, 1], tri.simplices, linestyle='--')
    for filter, color in zip([num_neigh <= 5, num_neigh == 6, num_neigh >= 7], ['lime', 'crimson', 'gold']):
        plt.scatter(points[filter, 0], points[filter, 1], s=80, color=color, zorder=2)
    plt.gca().set_aspect('equal', 'box') # equal axes needed because Delaunay depends on the scales
    plt.show()

example plot

Related