How to find delaunay triangulation facets containing the given points

Question

I have plotted n random points (the black points) and used delaunay triangulation, now I want to interpolate m random evaluation points (the red points) so I need to calculate which triangle the evaluation point is inside.

What is the approach for calculating the vertices of the triangle for each point?

Answer 1

For a given triangle, ABC, a point is inside the triangle if it is on the same side of line AB as point C is, on the same side of line BC as point A is, and on the same side of line AC as point B is. You can pre-optimize this check for each triangle and check them all until you find the triangle(s) it is in. See this page for more details.

To save computation, you can compute the minimum and maximum X and Y coordinates of the points for each triangle. If the X and Y coordinates of a point are not within the minimum and maximum values, you can immediately skip checking that triangle. The point cannot be inside it if it isn't inside the rectangle that bounds the triangle.

Answer 2

I'll assume that triangles do not intersect except of common edges.

You don't want to check every triangle (or subset of them) independently. The main reason is computation errors - due to them you may get answer "inside" for more than one triangle (or zero) which may break logic of your program.

More robust way is:

1. Find closest edge to the point
2. Select one of triangles touching this edge
3. Make one check for that triangle (the point lies on the same side as the third triangle vertex)
4. If "inside" - return this triangle
5. If "outside" - return another triangle on this edge (or "nothing" if there is no other triangle)

Even if you will return wrong triangle because of computation error, there still be exactly one triangle and point will lie close enough to it to accept such mistakes.

For #1 you can use something like quad-tree as Michael Wild suggests.

Answer 3

This simple example triangulates 10 random points, a further 3 random points are generated and if they fall inside a triangle the vertices are given:

import numpy as np
from pyhull.delaunay import DelaunayTri

def sign(a,b,c):
    return (a[0]-c[0])*(b[1]-c[1])-(b[0]-c[0])*(a[1]-c[1])

def findfacet(p,simplice):
    c,b,a = simplice.coords
    b1 = sign(p,a,b) < 0.0
    b2 = sign(p,b,c) < 0.0
    b3 = sign(p,c,a) < 0.0
    return b1 == b2 == b3

data = np.random.randn(10, 2)
dtri = DelaunayTri(data)

interpolate = np.random.randn(3, 2)

for point in interpolate:
    for triangle in dtri.simplices:
        if findfacet(point,triangle):
            print "Point",point,"inside",triangle.coords
        break

Using matplotlib to visualize (code omitted):

The dotted cyan lines now connect the points to interpolate with the vertices of triangle it lays within. The black lines are the convex hull, and the solid cyan lines are the delaunay triangulation.

Answer 4

A Delaunay triangulation is in itself a search data structure. Your Delaunay triangulation implementation probably has location functions. How have you computed the Delaunay triangulation of your points?

CGAL has an implementation of 2D and 3D triangulations. The resulting data structure is able to localize any point using a walk from a given point. See for example that chapter of the manual. CGAL is a C++ library, but it has python bindings.