Triangular mesh queries in Python

Theese kinds of queries became quite easy and effiecient with improved face based data structure which is used by CGAL. Here I have implemented code to valk around one specific vertex:

# The demonstration of improved face based data structure

from numpy import array

triangles = array([[ 5,  7, 10],
                  [ 7,  5,  6],
                  [ 4,  0,  3],
                  [ 0,  4,  6],
                  [ 4,  7,  6],
                  [ 4,  9, 10],
                  [ 7,  4, 10],
                  [ 0,  2,  1],
                  [ 2,  0,  6],
                  [ 2,  5,  1],
                  [ 5,  2,  6],
                  [ 8,  4,  3],
                  [ 4, 11,  9],
                  [ 8, 11,  4],
                  [ 9, 11,  3],
                  [11,  8,  3]], dtype=int)

points = array([[ 0.95448092,  0.45655774],
                [ 0.86370317,  0.02141752],
                [ 0.53821089,  0.16915935],
                [ 0.97218064,  0.72769053],
                [ 0.55030382,  0.70878147],
                [ 0.34692982,  0.08765148],
                [ 0.46289581,  0.29827649],
                [ 0.21159925,  0.39472549],
                [ 0.61679844,  0.79488884],
                [ 0.4272861 ,  0.93375762],
                [ 0.12451604,  0.54267654],
                [ 0.45974728,  0.91139648]])

import pylab as plt

fig = plt.figure()
pylab.triplot(points[:,0],points[:,1],triangles)

for i,tri in enumerate(triangles):
    v1,v2,v3 = points[tri]
    vavg = (v1 + v2 + v3)/3
    plt.text(vavg[0],vavg[1],i)

#plt.show()

## constructing improved face based data structure

def edge_search(v1,v2,skip):
    """
    Which triangle has edge with verticies i and j and aren't triangle <skip>?
    """
    neigh = -1
    for i,tri in enumerate(triangles):
        if (v1 in tri) and (v2 in tri):
            if i is skip:
                continue
            else:
                neigh = i
                break

    return(neigh)


def triangle_search(i):
    """
    For given vertex with index i return any triangle from neigberhood
    """
    for i,tri in enumerate(triangles):
        if i in tri:
            return(i)

neighberhood = []
for i,tri in enumerate(triangles):

    v1, v2, v3 = tri

    t3 = edge_search(v1,v2,i)
    t1 = edge_search(v2,v3,i)
    t2 = edge_search(v3,v1,i)

    neighberhood.append([t1,t2,t3])

neighberhood = array(neighberhood,dtype=int)

faces = []
for vi,_ in enumerate(points):
    faces.append(triangle_search(vi))

## Now walking over first ring can be implemented
def triangle_ring(vertex):

    tri_start = faces[vertex]
    tri = tri_start

    ## with asumption that vertex is not on the boundary
    for i in range(10):

        yield tri

        boolindx = triangles[tri]==vertex

        # permutating to next and previous vertex
        w = boolindx[[0,1,2]]
        cw = boolindx[[2,0,1]]
        ccw = boolindx[[1,2,0]]

        ct = neighberhood[tri][cw][0]

        if ct==tri_start:
            break
        else:
            tri=ct

for i in triangle_ring(6):
    print(i)

## Using it for drawing lines on plot

vertex = 6
ring_points = []

for i in triangle_ring(vertex):
    vi = triangles[i]
    cw = (vi==vertex)[[2,0,1]] 
    print("v={}".format(vi[cw][0]))
    ring_points.append(vi[cw][0])

data = array([points[i] for i in ring_points])
plt.plot(data[:,0],data[:,1],"ro")
#plt.savefig("topology.png")
plt.show()

input("Press Enter to continue...")
plt.close("all")