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

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.

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. • So in this approach each point would "know" its neighbors (sorry for possible lack of appropriate terminology), and then select the point from that list such that it forms the smallest angle as the next point? Sounds good so far, haven't thought of a reason why this wouldn't work.But since I'd like to run through the loop of deleting external triangles multiple times, is there any faster way to check just a single edge? I guess I could compute the hull once at the beginning, and then every time something is deleted compute it again for the chunk that was deleted. Jan 16, 2015 at 19:20

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.

• Nice, your approach is simpler than mine (and I like simpler solutions). It will produce different results if there are "holes" in the polygon (it is actually not clear from the question if the OP has to deal with that case, or what he expects for it). My approach does not only find the outer edges, it gives them an order, maybe that's more than the OP needs, I don't know. Jan 25, 2015 at 8:06

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)

# 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]))
# If there are more than one rightmost point, choose one
if rightIND.size >1:
rightIND = rightIND
# 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)]
bb = 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]))
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)]
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)
NL[row,col] = row
KL[row,col] = 1
col = np.where(KL[row,:]==0)
NL[row,col] = row
KL[row,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, row]
BL[dmyi+1] = [row, row]
BL[dmyi+2] = [row, row]
dmyi += 3

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

return BLtrim
``````
• You've got a lot of code there and not much explanation. Programmers.SE is a site that focuses on the conceptual side of software development rather than the raw code. Just posting raw code to solve a problem seems to me to be something more appropriate to a gist or repository on a code sharing site.
– user40980
Feb 8, 2016 at 14:11