Test Whether A Point Is Inside A Complex Polygon With Spline Curves

Point-In-Spline-Polygon Algorithm — Testing Whether A Point Is Inside A Complex Polygon With Spline Curves
©2006 Darel Rex Finley. This complete article, unmodified, may be freely distributed for educational purposes.

•		If you just need a code sample, skip to the bottom of this page.
•		Before reading this page, you should already be familiar with the ordinary point-in-polygon method found here.
•		For code to determine the horizontal and vertical extent (span) of a spline curve, go here.

A spline curve (also known as a three-point bezier curve) can be simply defined with three points as follows:

spline curve definition

Sx,Sy is the start point, Ex,Ey is the endpoint, and a,b is the spline point. Now imagine that this is a piece of “string art” from the 1970s, and all the points along the first black line [Sx,Sy to a,b] are connected with blue strings to all the points along the second black line [a,b to Ex,Ey]. The lines are connected at their proportional points; i.e. the point that is 10% of the way from Sx,Sy to a,b is connected to the point that is 10% of the way from a,b to Ex,Ey. Now imagine that we place a blue pushpin (represented here as a blue disc) on each string, again at the proportional point. So, for example, on our 10% string we put a pushpin at the point that is 10% along the blue string from the first black line to the second black line.

The result is that the blue pushpins form a nice, graceful curve that happens to be tangent to the two black lines where it touches them at Sx,Sy and Ex,Ey. We want to be able to include such curves in our point-in-polygon function.

One way to accomplish that would be to generate an arbitrary number of pushpin points (say 10 or 20 or 50 of them), then use them as corners of our polygon. That method works, but has a big tradeoff between the smoothness of the curve and the speed with which the polygon can be tested. What we really want is a way of incorporating the mathematical definition of the curve into our point-in-polygon algorithm so it can just test the point against the entire spline at once. Can that be done? Sure it can! Let’s dredge up our old algebra skills and have at it:

If F is the fractional value (along the lines and strings), which ranges from 0 to 1 (e.g. 0.1 would be 10%), then the y-coordinate of the point along the first black line is

    Sy + F (b-Sy)

and the y-coordinate of the point along the second black line is

    b + F (Ey-b)

We will call our pushpin point X,Y — so the y-coordinate of our pushpin point along the blue string will be

    Y = firstPointY + F (secondPointY - firstPointY)

which is

    Y = Sy + F (b-Sy) + F ( b + F (Ey-b) - (Sy + F (b-Sy)) )

Now separate out F as much as possible, preparing for application of the quadratic formula:

    Y = Sy + F b - F Sy + F ( b + F Ey - F b - Sy - F b + F Sy )

    Y = Sy + F b - F Sy + F b + F² Ey - F² b - F Sy - F² b + F² Sy

    0 = (Sy - Y) + F (2 (b-Sy)) + F² (Sy + Ey - 2 b)

The quadratic formula tells us that if

    0 = c + bx + ax²

then

    x = ( -b [+ or -] sqrt(b² - 4ac) ) / (2a)

Applying this to what we have so far, we get

    F = ( -2 (b-Sy) [+ or -] sqrt( (2 (b-Sy))² - 4 (Sy + Ey - 2 b) (Sy - Y) ) ) / (2 (Sy + Ey - 2 b))

We know everything in this formula except F. (We are using the y-coordinate of our test point for Y.) So we can use the formula to get F, and once we have F, we can plug it back into our original formula (for X this time, not Y) to get an X value that can be compared to the x-coordinate of our test point:

    X = Sx + F (a-Sx) + F ( a + F (Ex-a) - (Sx + F (a-Sx)) )

    if X < testX then flip oddNodes

But to write this code we have to be aware of a couple of things. First, when we apply the quadratic formula, if the part of the equation that lies inside the sqrt is negative, then F is undefined, and so our spline does not cross the y-threshold of the test point. Secondly, if F is defined, then it has two values (thanks to the plus-or-minus in the quadratic formula). Each of these two values must be independently tested against our test point — but first each of the two F values must be checked to see if it lies in the range 0 to 1 — if not, it must be ignored, since it is off the end of the spline. (And one more thing: If the two F values are identical, they both still count, so don’t bother comparing them to each other.)

Avoiding Special Cases

What happens if F is exactly 0, or exactly 1? This opens up a whole can of headaches that we’d rather not deal with, for the sake of simpler code and better execution speed. Probably the easiest way to avoid the problem is just to add a very small value (say, 0.000001) to the test point’s y-coordinate before testing the point. That will pretty much guarantee that F will never be exactly 0 or 1.

What about division by zero? Can the quadratic’s “2a” part be zero in this usage, other than in a defective case, like multiple corners on the same spot? Actually, yes it easily can. The “a” part of the quadratic is, in this case, Sy+Ey-2b — that means that if the spline corner is (vertically) exactly inbetween the startpoint and endpoint, which can easily happen (see the house-shaped polygon below), then there will be division by zero! To prevent this, we’ll just do a similar quick-fix stunt: when (and only when) the “a” part is zero, we’ll add 0.0001 to the spline corner’s y-coordinate. Problem solved!

Code Sample

This C code sample is designed around the assumption that the entire polygon resides inside the square that spans from -1 to +1 on each axis. Values outside that range are used as special tags:
    2 = the next coordinate pair is a spline point
    3 = a new loop is starting, but it’s part of the same shape
    -2 = terminator (the shape is done)
See the sample polygons at the bottom of this page for examples of how these tags are used.

Speed: This code sample is fine for testing a point or two, but if you want to test a lot of points (such as the massive, raster testing I used to make the sample polygon images below), then this function would be horribly inefficient, because it performs the same calculations on the polygon’s data every time it’s called. In that case, you would probably want to use the efficient fill technique described here.

When this function encounters two consecutive spline corners, it interpolates a hard corner equidistant between them. That’s pretty standard for most programs that deal with splines.

//  public domain function by Darel Rex Finley,
2006



bool pointInSplinePoly(double *poly, double X, double Y) {



  #define  SPLINE     2.

  #define  NEW_LOOP   3.

  #define  END       -2.



  double  Sx, Sy, Ex, Ey, a, b, sRoot, F, plusOrMinus, topPart, bottomPart, xPart ;

  int     i=0, j, k, start=0 ;

  bool    oddNodes=NO ;



  Y+=.000001;   //  prevent the need for special tests when F is exactly 0 or
1



  while (poly[i]!=END) {



    j=i+2; if (poly[i]==SPLINE) j++;

    if (poly[j]==END    || poly[j]==NEW_LOOP) j=start;



    if (poly[i]!=SPLINE && poly[j]!=SPLINE  ) {   //  STRAIGHT
LINE

      if (poly[i+1]<Y && poly[j+1]>=Y || poly[j+1]<Y && poly[i+1]>=Y) {

        if (poly[i]+(Y-poly[i+1])/(poly[j+1]-poly[i+1])*(poly[j]-poly[i])<X)
oddNodes=!oddNodes; }}



    else if (poly[j]==SPLINE)
{                   //  SPLINE CURVE

      a=poly[j+1]; b=poly[j+2]; k=j+3; if (poly[k]==END || poly[k]==NEW_LOOP) k=start;

      if (poly[i]!=SPLINE) {

        Sx=poly[i]; Sy=poly[i+1]; }

      else {   //  interpolate hard corner

        Sx=(poly[i+1]+poly[j+1])/2.; Sy=(poly[i+2]+poly[j+2])/2.; }

      if (poly[k]!=SPLINE) {

        Ex=poly[k]; Ey=poly[k+1]; }

      else {   //  interpolate hard corner

        Ex=(poly[j+1]+poly[k+1])/2.; Ey=(poly[j+2]+poly[k+2])/2.; }

      bottomPart=2.*(Sy+Ey-b-b);

      if (bottomPart==0.) {   //  prevent
division-by-zero

        b+=.0001; bottomPart=-.0004; }

      sRoot=2.*(b-Sy); sRoot*=sRoot; sRoot-=2.*bottomPart*(Sy-Y);

      if (sRoot>=0.) {

        sRoot=sqrt(sRoot); topPart=2.*(Sy-b);

        for (plusOrMinus=-1.; plusOrMinus<1.1; plusOrMinus+=2.) {

          F=(topPart+plusOrMinus*sRoot)/bottomPart;

          if (F>=0. && F<=1.) {

            xPart=Sx+F*(a-Sx); if (xPart+F*(a+F*(Ex-a)-xPart)<X)
oddNodes=!oddNodes; }}}}



    if (poly[i]==SPLINE) i++;

    i+=2;

    if (poly[i]==NEW_LOOP) {

      i++; start=i; }}



  return oddNodes; }

Polygon Samples

These polygon samples are designed for use with the above function.

This five-cornered polygon makes a simple house shape: `1,0, 1,-1, -1,-1, -1,0, 0,1, END`
This is the same house shape, but with the top corner turned into a spline point: `1,0, 1,-1, -1,-1, -1,0, SPLINE,0,1, END`
Same house shape, but with the two roof’s-eve corners turned into spline points: `SPLINE,1,0, 1,-1, -1,-1, SPLINE,-1,0, 0,1, END` Note that this case would not work without the division-by-zero protection, because the eve points are vertically at the midpoint between the top and bottom of the house.
Same house shape, with the two bottom corner points turned into spline points: `1,0, SPLINE,1,-1, SPLINE,-1,-1, -1,0, 0,1, END`
Same house shape, with all five corners turned into spline points: `SPLINE,1,0, SPLINE,1,-1, SPLINE,-1,-1, SPLINE,-1,0, SPLINE,0,1, END`
Here’s an example of using the “NEW_LOOP” tag to put two polygons in one. In this case, the smaller polygon is a hole in the larger one, since it is inside it: `SPLINE,1,1, SPLINE,1,-1, SPLINE,-1,-1, SPLINE,-1,1, NEW_LOOP, SPLINE,.5,.5, SPLINE,.5,-.5, SPLINE,-.5,-.5, SPLINE,-.5,.5, END`
This uses the “NEW_LOOP” tag to make two loops, but they don’t overlap each other: `-.25,0, .125,0, .125,-.875, .25,-.875, .25,-1, -.25,-1, -.25,-.875, -.125,-.875, -.125,-.125, -.25,-.125, NEW_LOOP, SPLINE,-.125,.375, SPLINE,.125,.375, SPLINE,.125,.25, SPLINE,-.125,.25, END`
You can’t do a real circle with splines, but this eight-point polygon makes a sweet fake: `SPLINE,.4142,1, SPLINE,1,.4142, SPLINE,1,-.4142, SPLINE,.4142,-1, SPLINE,-.4142,-1, SPLINE,-1,-.4142, SPLINE,-1,.4142, SPLINE,-.4142,1, END` FYI, .4142 is sqrt(2)-1.

Send me an e-mail!

Does the brace style in the above code sample freak you out? Click here to see it explained in a new window.