Point-in-polygon: Jordan Curve Theorem

From Sidvind
Jump to: navigation, search

Calculating whenever a point is inside a polygon can sometimes be a hard and costly calculation. This article describes a quite cheap solution to calculate whenever a point is inside ANY closed polygon. In an open polygon it's hard to determine what's in and out so naturally it won't work.

A closed polygon with 3 points marked

The Jordan Curve Theorem states that a point is inside a polygon if the number of crossings from an arbitrary direction is odd. An image explains more than a thousand words so lets take a look at the picture. As you can see point 1 and 3 is inside the polygon but point 2 isn't. Follow the rays from each point and count each time you cross a line-segment. In this article I only deal with 2D polygons but it can easily used in a 3D-environment.

Find crossings[edit]

Casting a ray[edit]

One of the first things to do is to cast a ray from the point in an arbitrary direction. I use a ray along the Y axis (pointing upwards as in the picture) for simplicity. Along X-axis is good too, but use one of those or it gets a lot harder. Remember that I use the ray mentioned above throughout this article.

Finding the equation of a line-segment[edit]

As a first step, if the ray is along y-axis, check if the point x-coordinate is between the two points connecting the line. If not it don't cross it either. You can also check if the y-coordinate is above both points. The next step is to find the equation of the line. Hopefully you remember this from grade school. The equation of a straight line is $ y=kx+m\, $ (Swedish notation). The slope is $ k=\frac{\Delta\ y}{\Delta\ x} $ and offset is $ m=y-kx\, $. Do the math we have the equation. Now, insert the x-coordinate of the point into the equation. If the result is larger than the y-coordinate the ray does not cross the line-segment.

Repeat this for each line-segment.

C/C++ implementation[edit]

Code: Sample (naive) implementation

  1. /* The points creating the polygon. */
  2. float x[8];
  3. float y[8];
  4. float x1,x2;
  5.  
  6. /* The coordinates of the point */
  7. float px, py;
  8.  
  9. /* How many times the ray crosses a line-segment */
  10. int crossings = 0;
  11.  
  12. /* Coordinates of the points */
  13. x[0] = 100;     y[0] = 100;
  14. x[1] = 200;     y[1] = 200;
  15. x[2] = 300;     y[2] = 200;
  16. x[3] = 300;     y[3] = 170;
  17. x[4] = 240;     y[4] = 170;
  18. x[5] = 240;     y[5] = 90;
  19. x[6] = 330;     y[6] = 140;
  20. x[7] = 270;     y[7] = 30;
  21.  
  22. /* Iterate through each line */
  23. for ( int i = 0; i < 8; i++ ){
  24.        
  25.         /* This is done to ensure that we get the same result when
  26.            the line goes from left to right and right to left */
  27.         if ( x[i] < x[ (i+1)%8 ] ){
  28.                 x1 = x[i];
  29.                 x2 = x[(i+1)%8];
  30.         } else {
  31.                 x1 = x[(i+1)%8];
  32.                 x2 = x[i];
  33.         }
  34.        
  35.         /* First check if the ray is possible to cross the line */
  36.         if ( px > x1 && px <= x2 && ( py < y[i] || py <= y[(i+1)%8] ) ) {
  37.                 static const float eps = 0.000001;
  38.  
  39.                 /* Calculate the equation of the line */
  40.                 float dx = x[(i+1)%8] - x[i];
  41.                 float dy = y[(i+1)%8] - y[i];
  42.                 float k;
  43.  
  44.                 if ( fabs(dx) < eps ){
  45.                         k = INFINITY;   // math.h
  46.                 } else {
  47.                         k = dy/dx;
  48.                 }
  49.  
  50.                 float m = y[i] - k * x[i];
  51.                
  52.                 /* Find if the ray crosses the line */
  53.                 float y2 = k * px + m;
  54.                 if ( py <= y2 ){
  55.                         crossings++;
  56.                 }
  57.         }
  58. }
  59.  
  60. printf("The point is crossing %d lines", crossings);
  61. if ( crossings % 2 == 1 ){
  62.         printf(" thus it is inside the polygon");
  63. }
  64. printf("\n");