4
\$\begingroup\$

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
\$\endgroup\$
  • \$\begingroup\$ Can we use matplotlib/numpy, or do you want pure python? \$\endgroup\$ – brechmos Jun 17 '14 at 1:10
  • \$\begingroup\$ @brechmos No rules at all! Let's use whatever we can to get as much performance as we can! \$\endgroup\$ – user89 Jun 17 '14 at 1:13
  • 1
    \$\begingroup\$ You are unnecessarily complicating things. 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
  • \$\begingroup\$ @Jaime Hm. How come that algorithm works by always choosing the same fixed direction? \$\endgroup\$ – user89 Jun 17 '14 at 14:37
  • \$\begingroup\$ @Jaime Ah! It seems that that methodology doesn't work for convex polygons... \$\endgroup\$ – user89 Jun 17 '14 at 14:52
1
\$\begingroup\$

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. :-)

\$\endgroup\$
0
\$\begingroup\$

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...

\$\endgroup\$
  • 1
    \$\begingroup\$ would you mind including at least pseudocode? I think no one really wants to skim a 160 pages dissertation for an algorithm. Currently this answer is also no more than a link. While I believe you might be correct, I'd also like to see the contents of the link summarized, in case the link goes down... \$\endgroup\$ – Vogel612 Jul 1 '14 at 6:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.