Method for determining if a line segment is an external edge of a Delauney triangulation?

Question

I've created a Delauney triangulation of a set of points. Now I want to iterate over the triangulation and remove any line segments/edges on the exterior for which the following are true:

1. The outside edge is the hypotenuse of the triangle
2. Call the hypotenuse the "base" and determine the "height" of the triangle. Then if base/height > x (where x is some criteria value, say 4), delete that edge/triangle (depending on data structure/object definitions).

In essence I'm trying to remove edges from the triangulation to create a more characteristic shape by deleting edges that are long, skinny, and on the outside of the triangulation. I might add in a check later for size of the triangle relative to the average of the rest, too.

My question is: is there an easy way to check if a given edge is on the exterior of the triangulation?

I figured I could build in to the triangulation process itself a design that will make every line segment its own object, then have each line segment "know" its neighboring polygons, and if it has an empty neighbor slot then it is external. This seemed inefficient, though, and the code I currently have just spits out a bunch of triangle objects with its vertices's indexes stored. Calculating the line segments shouldn't be too bad, but then I only want to do the deletion process on external edges as opposed to all of them.

Not sure if this question belonged in math or programming, apologies if it should not be in here or if I provided too little/too much information.

I'm working in Python, but my concern here is about theoretical approach/algorithm as opposed to code so I don't feel it necessary to include code.

EDIT to note: I'll want to iterate over the exterior multiple times in the event that after deleting one exterior triangle the first round through it reveals a new exterior that fits the deletion parameters. Therefore I'd like to avoid completely recalculating the exterior edges each time -- thus the question about a quick-and-easy, if it exists.

Answer 1

A triangulation of your point set gives you for each point a list of "adjacent" points, the neighbours to which a point is connected.

Once you have calculated the triangulation, you can determine the outer edges by applying an algorithm similar to the classic gift-wrapping algorithm for convex hulls. You start with, for example, the "rightmost" point (which is part of the convex hull and so at the outer edge). Then step from point to point using only outer edges of the triangulation. The only difference to the "convex hull" algorithm is, when having points p_(n-1) and p_(n), you choose p_(n+1) as the point among all adjacent points minimizing the angle between p_(n)->p_(n-1) and p_(n) ->p_(n+1).

Hope this helps.

Answer 2

It's fortunate you can start with a list of triangle objects, because that makes it easy:

• flatMap the list of triangles into a list of edges.
• Your outer edges appear only once in the list.

You can see how this works by looking at a picture. Each inner edge is shared by exactly two triangles. Each outer edge belongs to exactly one triangle.

When you remove an edge, the other two edges in that triangle always become outer edges, if they weren't already. Just add them to the set of outer edges before you remove the triangle from its list. This can be sped up by keeping a map from edges back to triangles.

Answer 3

For future readers, here is a python function (using numpy) that creates a list of the indices of points on the boundary.

I am starting with a general network of connected points (not necessarily triangulated), so my data is in the form of (xy, NL, KL) with NL being the neighbor list and KL being the connectivity list (see below).

def extract_boundary_from_NL(xy,NL,KL):
    '''
    Extract the boundary of a 2D network (xy,NL,KL).

    Parameters
    ----------
    xy : #pts x 2 float array
        point set in 2D
    NL : #pts x max # nearest neighbors int array
        Neighbor list. The ith row contains the indices of xy that are the bonded pts to the ith pt.
        Nonexistent bonds are replaced by zero.
    KL : #pts x max # nearest neighbors int array
        Connectivity list. The jth column of the ith row ==1 if pt i is bonded to pt NL[i,j].
        The jth column of the ith row ==0 if pt i is not bonded to point NL[i,j].

    Returns
    ----------
    boundary : #points on boundary x 1 int array
        indices of points living on boundary of the network

    '''
    # Initialize the list of boundary indices to be larger than necessary
    bb = np.zeros(len(xy), dtype=int)

    # Start with the rightmost point, which is guaranteed to be 
    # at the convex hull and thus also at the outer edge.
    # Then take the first step to be along the minimum angle bond
    rightIND = np.where(xy[:,0]== max(xy[:,0]))[0]
    # If there are more than one rightmost point, choose one
    if rightIND.size >1:
        rightIND = rightIND[0]
    # Grab the true neighbors of this starting point
    neighbors = NL[rightIND,np.argwhere(KL[rightIND]).ravel()]
    # Compute the angles of the neighbor bonds 
    angles = np.mod(np.arctan2(xy[neighbors,1]-xy[rightIND,1],xy[neighbors,0]-xy[rightIND,0]).ravel(), 2*np.pi)
    # Take the second particle to be the one with the lowest bond angle (will be >= pi/2)
    nextIND = neighbors[angles==min(angles)][0]
    bb[0] = rightIND

    dmyi = 1
    # as long as we haven't completed the full outer edge/boundary, add nextIND
    while nextIND != rightIND:
        bb[dmyi] = nextIND
        n_tmp = NL[nextIND,np.argwhere(KL[nextIND]).ravel()]
        # Exclude previous boundary particle from the neighbors array
        # since its angle with itself is zero!
        neighbors = np.delete(n_tmp, np.where(n_tmp == bb[dmyi-1])[0])
        angles = np.mod( np.arctan2(xy[neighbors,1]-xy[nextIND,1],xy[neighbors,0]-xy[nextIND,0]).ravel() \
                - np.arctan2( xy[bb[dmyi-1],1]-xy[nextIND,1], xy[bb[dmyi-1],0]-xy[nextIND,0] ).ravel(), 2*np.pi)
        nextIND = neighbors[angles==min(angles)][0]
        dmyi += 1

    # Truncate the list of boundary indices
    boundary = bb[0:dmyi]

    return boundary

To put a triangulation into the form (NL, KL) to use the above, you could do:

BL = TRI2BL(TRI)
NL,KL = BL2NLandKL(BL,NN=10)



def BL2NLandKL(BL, NP='auto', NN=6):
    '''Convert bond list (#bonds x 2) to neighbor list (NL) and connectivity list (KL) for lattice of bonded points.
    Returns KL as ones where there is a bond and zero where there is not.
    (Even if you just want NL from BL, you have to compute KL anyway.)
    Note that this makes no attempt to preserve any previous version of NL, which in the philosophy of these simulations should remain constant during a simulation.
    If NL is known, use BL2KL instead, which creates KL according to the existing NL. 

    Parameters
    ----------
    BL : array of dimension #bonds x 2
        Each row is a bond and contains indices of connected points
    NP : int
        number of points (defines the length of NL and KL)
    NN : int
        maximum number of neighbors (defines the width of NL and KL)

    Returns
    ----------
    NL : array of dimension #pts x max(#neighbors)
        The ith row contains indices for the neighbors for the ith point.
    KL :  array of dimension #pts x (max number of neighbors)
        Spring constant list, where 1 corresponds to a true connection while 0 signifies that there is not a connection.
    '''
    if NP=='auto':
        if BL.size>0:
            NL = np.zeros((max(BL.ravel())+1,NN),dtype=np.intc)
            KL = np.zeros((max(BL.ravel())+1,NN),dtype=np.intc)
        else:
            raise RuntimeError('ERROR: there is no BL to use to define NL and KL, so cannot run BL2NLandKL()')
    else:
        NL = np.zeros((NP,NN),dtype=np.intc)
        KL = np.zeros((NP,NN),dtype=np.intc)

    if BL.size > 0:
        for row in BL:
            col = np.where(KL[row[0],:]==0)[0][0]
            NL[row[0],col] = row[1]
            KL[row[0],col] = 1
            col = np.where(KL[row[1],:]==0)[0][0]
            NL[row[1],col] = row[0]
            KL[row[1],col] = 1

    return NL, KL

def TRI2BL(TRI):
    '''
    Convert triangulation index array (Ntris x 3) to Bond List (Nbonds x 2) array.

    Parameters
    ----------
    TRI : Ntris x 3 int array
        Triangulation of a point set. Each row gives indices of vertices of single triangle.

    Returns
    ----------
    BL : Nbonds x 2 int array 
        Bond list

    '''
    # each edge is shared by 2 triangles unless at the boundary.
    # each row contains 3 edges.
    # An upper bound on the number bonds is 3*len(TRI)
    BL = np.zeros((3*len(TRI),2),dtype = int)

    dmyi = 0
    for row in TRI:
        BL[dmyi] = [row[0], row[1]]
        BL[dmyi+1] = [row[1], row[2]]
        BL[dmyi+2] = [row[0], row[2]]
        dmyi += 3

    # Sort each row to be ascending
    BL_sort = np.sort(BL, axis=1)
    BLtrim = unique_rows(BL_sort)

    return BLtrim