Calculating if a point is within a polygon, or outside of it

Question

Here is the function (along with its support functions):

def findAnglesBetweenTwoVectors(v1s, v2s):
    dot_v1_v2 = np.einsum('ij,ij->i', v1s, v2s)
    dot_v1_v1 = np.einsum('ij,ij->i', v1s, v1s)
    dot_v2_v2 = np.einsum('ij,ij->i', v2s, v2s)

    return np.arccos(dot_v1_v2/(np.sqrt(dot_v1_v1*dot_v2_v2)))

def calculateDirectionOfRayNotParallelToAnyEdgeOfPolygon(vertices, difference = 1e-2, referenceDirectionVector = np.array([1, 0]), maxWhileTime = 0.5):

    relevantVectors = vertices - np.roll(vertices, 1, axis = 0)
    angles = findAnglesBetweenTwoVectors(relevantVectors, np.repeat([referenceDirectionVector], len(relevantVectors), axis=0))
    angles = np.append(angles, np.pi+angles)

    for x in xrange(50):
        direction = random.random()*2*np.pi
        if np.all(np.abs(angles-direction) >= difference):
            return direction

    raise StandardError("No direction found!")

def calculateNumberOfBoundaryIntersections(point, direction, vertices):

    neighbourVertices = np.roll(vertices, 1, axis = 0)
    ss = neighbourVertices - vertices

    r = produceDirectionVectorGivenDirection(direction)

    qs_minus_p = vertices - point
    r_cross_ss = np.cross(r, ss)

    ts = np.cross(qs_minus_p, ss)/r_cross_ss
    if 0 in ts:
        return 1

    us = np.cross((qs_minus_p), r)/r_cross_ss

    return np.sum((ts >= 0) & (us >= 0) & (us < 1))

def isPointWithinPolygon(point, polygonVertexCoords):
    fixedDirection = calculateDirectionOfRayNotParallelToAnyEdgeOfPolygon(polygonVertexCoords)

    numberOfBoundaryIntersections = calculateNumberOfBoundaryIntersections(point, fixedDirection, polygonVertexCoords)

    if numberOfBoundaryIntersections%2 == 0:
        return False
    else:
        return True

The program itself is based on a simple discrete version of the Jordan curve theorem: if a point is inside of a polygon, then a ray emanating from it in a direction that is not parallel to any of the edges of the polygon will cross the polygon boundary an odd number of times. If the point is not inside the polygon, the ray will cross the polygon boundary an even number of times.

The boundary intersection calculation logic is described here and then here.

Earlier, I used a different algorithm which spent a lot of time in findAnglesBetweenTwoVectors. Through a code review I performed here, we determined that it was unlikely we could optimize that function anymore. After I changed the boundary intersection logic, I spend a lot more time in calculateNumberOfBoundaryIntersections, with a small performance gain over the old algorithm.

So, while it is unlikely that performance can be squeezed out of findAnglesBetweenTwoVectors I suppose the other functions are still up for grabs.

Once again, this is a place where my code spends a lot of time (75% of the total time!), so optimizing it could lead to major performance gains for me. Is my NumPy usage good?

EDIT: I just achieved a moderate gain by changing the organization of the code so that it is like this:

def determineUnitVectorNotParallelToAnyEdgeOfPolygon(vertices):
    edges = np.roll(vertices, 1, axis = 0) - vertices

    for x in range(100):
        direction = 2*np.pi*random.random()
        directionVector = np.array([np.cos(direction), np.sin(direction)])

        if 0 not in np.cross(edges, directionVector):
            return directionVector

def calculateNumberOfBoundaryIntersections(point, directionVector, vertices):

    neighbourVertices = np.roll(vertices, 1, axis = 0)
    ss = neighbourVertices - vertices

    r = directionVector

    qs_minus_p = vertices - point
    r_cross_ss = np.cross(r, ss)

    ts = np.cross(qs_minus_p, ss)/r_cross_ss
    if 0 in ts:
        return 1

    us = np.cross((qs_minus_p), r)/r_cross_ss

    return np.sum((ts >= 0) & (us >= 0) & (us < 1))

def isPointWithinPolygon(point, polygonVertexCoords):
    numberOfBoundaryIntersections = calculateNumberOfBoundaryIntersections(point, determineUnitVectorNotParallelToAnyEdgeOfPolygon(polygonVertexCoords), polygonVertexCoords)

    if numberOfBoundaryIntersections%2 == 0:
        return False
    else:
        return True

Answer 1

I had to write an roipoly in python/numpy/matplotlib as I use that a lot for my work (and I liked the one in Matlab). So here is my implementation of roipoly:

def roipoly():
    XY = ginput(-1)

    nrows = sum( abs( r_[ gcf().get_axes()[0].get_xbound() ] ) )
    ncols = sum( abs( r_[ gcf().get_axes()[0].get_ybound() ] ) )
    nrows,ncols = int(nrows), int(ncols)

    path = matplotlib.path.Path( XY )
    mask = array([ path.contains_point((j,i)) for i in range(nrows) for j in range(ncols)  ]).reshape(ncols,nrows)

    return mask  

If I understand the problem, I am guessing you could adapt so it would be something like:

def isPointWithinPolygon(point, polygonVertexCoords):
     path = matplotlib.path.Path( polygonVertexCoords )
     return path.contains_point(point[0], point[1]) 

I am not completely sure point and polygonVertexCoords are going to be in the same format, but maybe this would be a start.

I am not going to claim this is the most efficient method, but is a method. :-)

Answer 2

One may also solve for the generalized barycentric coordinates: at least one of them must be negative for points outside the convex hull of the polytope, see for example this thesis on On Generalized Barycentric Coordinates...