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
scikit-image
has a working implementation of the point in polygon algorithm in Cython here, and you may want to check the discussion here for details on the algorithm, although the code presented is all C. \$\endgroup\$ – Jaime Jun 17 '14 at 7:07