Point in polygon test is one of the most fundamental building block in ray tracing. Given a point on a surface and a vertices array representing a polygon, is the point in the polygon? It also has wide application in GIS systems, forming one of the sustrates of spatial joins.

Implementing the algorithm efficiently is crucial as ray-tracing is perpetually computation hungry. Eric Heines gave one efficient implementation in Graphics Gem online source.

int CrossingsMultiplyTest(double pgon[][2], int numverts, double point[2])
{
    int j, yflag0, yflag1, inside_flag;
    double ty, tx, *vtx0, *vtx1;

    tx = point[X];
    ty = point[Y];

    vtx0 = pgon[numverts - 1];
    yflag0 = (vtx0[Y] > ty);
    vtx1 = pgon[0];

    inside_flag = 0;
    for (j = numverts + 1; --j; ) {

        yflag1 = (vtx1[Y] > ty);
        if (yflag0 != yflag1) {
            if (((vtx1[Y] - ty) * (vtx0[X] - vtx1[X]) >
                (vtx1[X] - tx) * (vtx0[Y] - vtx1[Y])) == yflag1) {
                inside_flag = !inside_flag;
            }
        }

        yflag0 = yflag1;
        vtx0 = vtx1;
        vtx1 += 2;
    }

    return(inside_flag);
}

At first glance, part of the algorithm does look familiar, the inner-most if clause seems to be doing some handedness test by computing the determinent. But there are also many parts that’s alien, what is yflags?

Back to Basics

Let’s derive the ray-crossings test. Step 1 is to check if the test point is with the y range of the segment. Set test point \((t_x, t_y)\), segment as \((a_x, a_y)\) -> \((b_x, b_y)\). This is equivalent to checking:

\[a_y<t_y<b_y \space \lor \space b_y<t_y<a_y\]

From the test, we know that \(yflag0 \iff t_y>a_y\), \(yflag1 \iff t_y>by\). By writing a truth table between the relationship of yflags and our target:

                ty > ay    ty > by   ay<ty<by or by<ty<ay   ty>ay != ty>by
ty < ay < by       F          F               F                   F
ty < by < ay       F          F               F                   F
ay < by < ty       T          T               F                   F
by < ay < ty       T          T               F                   F
ay < ty < by       T          F               T                   T
by < ty < ay       F          T               T                   T

We thus proved that the checking yflag0 != yflag1 is equivalent to checking that ty is in the y range of the segment.

Note that in the discussion above, I have confounded boudnary conditions by omitting = sign in the truth table. Boudnary condition for this point in polygon algorithm is not gaurenteed to start with, so we tend to not discuss case when a point falls right on the boundary of the y range. In ray-tracing cases there is an infinitestimally small probability for the point to fall exactly on the boundary.

Step 2 is ray-segment intersection.

The ray is extending to the right direction of x axis, parallel to x axis. Set ray and segment equations:

\[A_1x+B_1y=C_1 \\ A_2x+B_2y=C_2\]

set \(rise = b_y-a_y \\ run = b_x-a_x\)

then \(A_1 = 0 \\ B_1 = 1 \\ C_1 = t_y \\ A_2 = rise \\ B_2 = -run \\ C_2 = rise\cdot a_x-run\cdot a_y \\\)

x coordinate of the intersection is:

\[x' = \frac{ \det\begin{bmatrix} B_1 & C_1 \\ B_2 & C_2 \\ \end{bmatrix}}{ \det\begin{bmatrix} A_1 & B_1 \\ A_2 & B_2 \\ \end{bmatrix} } = \frac{ \det\begin{bmatrix} B_1 & C_1 \\ B_2 & C_2 \\ \end{bmatrix}}{ \det\begin{bmatrix} A_1 & B_1 \\ A_2 & B_2 \\ \end{bmatrix} }=\frac{run\cdot a_y + (rise\cdot a_x - run \cdot a_y)}{rise}\]

Ray and segment have 1 intersection iff \(tx < x' \iff tx < \frac{run\cdot a_y + (rise\cdot a_x - run \cdot a_y)}{rise}\)

Since division is expensive to compute, depending on the sign of rise, this can be further reduced to:

\[tx < x' \iff (rise > 0 \land tx \cdot rise < run\cdot a_y + (rise\cdot a_x - run \cdot a_y)) \lor \\ (rise < 0 \land rise > run\cdot a_y + (rise\cdot a_x - run \cdot a_y))\]

However, the or predicate can introduce a branching in the code. In GPU kernel, branching introduces divergence and could reduce performance.

Notice that the sign of \(rise\) and the direction of the comparator is correlated.

Given \(rise>0\) is true, we expect \(tx \cdot rise < run\cdot a_y + (rise\cdot a_x - run \cdot a_y)\) to be true. And the second half of the predicate basically flips the direction of the comparator at the same time. Still ignoring the equal case, we can quasi-transcribe the second half as “given \(rise>0\) is false, we expect \(tx \cdot rise < run\cdot a_y + (rise\cdot a_x - run \cdot a_y)\) to be false”. That is, both inequaility should be true or false at the same time. This leads to the simplified form of the ray-crossings test:

\[tx < x' \iff tx \cdot rise < run\cdot a_y + (rise\cdot a_x - run \cdot a_y) == rise > 0\]

\(rise>0\) is the yflag1 we used previously, replacing it results in the finalized form.

cuSpatial Implementation

cuSpatial impelements this fast point-in-polygon primitive and is available as both C++ and Python API. A more advanced version utilized quadtree to index the points and filters out the qualified points before performing point in polygon test and is able to perform 1.3M+ points in 27 ROI in the matter of ~1ms (source). To know more about cuspatial, see documentation to get started.