This page is released under the GNU Free Documentation License, as any other text related to programming in this site; Images are released under the CreativeCommons Attribution - Noncommercial - Share Alike 3.0.
I've started this work as stream-of-consciousness to write a new task for RosettaCode. In that page I've used thoughts here expressed, altogether with images generated for this article. Due to collaborative efforts typical of a wiki, the page on RosettaCode might contain more updates than this article.

Points in polygons

One of the way to know if a point is inside a polygon is rather simple and worth studying. Let's suppose we want to know if the point P is inside a polygon; then we count how may times a horizontal ray starting from P (going towards increasing x coordinates) intersect the sides of the figure. If this number is even then P is outside, otherwise it is inside the polygon (the number 0 is considered even).

We are not interested in a proof (an intuitive one is here); rather, it can be recognized that if we want to implement such an algorithm, the main part of our implementation will be in the way we check if a ray intersect a segment. Having this subroutine, the code can be simply summarized:

  count ← 0
  foreach side in polygon:
    if ray_to_incr_x_from(P) intersects side then
      count ← count + 1
  if count is odd then
    return inside
  else
    return outside  

So the whole matter is implementing an intersects operator / function. Being not concerned with efficiency, let's see a simple way how we can see if a ray intersect or not a segment.

The segment we are considering can be bounded into a box delimited by the end points (A and B) of the segment. If the point is above, under or on the right of the box, then it surely does not intersect the segment (see points P1 and P2). So the hatched area contains all points not intersecting the segment for sure. On the contrary, points inside the greenish area (like P3), intersect for sure the segment.

A little bit harder are the triangular areas "up" and "down" the segment. One way to check if the point is above is to find the equation of the line passing through A and B, then assign to the "x" of the equation the x coordinate of the point we are checking (e.g. P4; this means to find the intersection point between the vertical line passing by P4 and the segment from A to B). If the y coordinate of the interesection point between the vertical line passing through P4 and the segment is less than the y coordinate of the point, then the point is in the triangle from where the ray will intersect the segment. Otherwise, it is in the other triangle. Of course the speech is the opposite if the slope of the segment is negative (point A is on the left of the point B); in this case, the triangle hosting points that will "cause" the intersection is the one in the bottom.

There are some quirks to take into account. First, we must consider the special case of vertical segments (when the x coordinate of A is the same of the x coordinate of B); and the special case of horizontal segments. They are two different problem. The case of the vertical line can be easily resolved. The case of a horizontal line creates more problems. The "green" area would be "crushed" in a line, and there will exist infinite point of intersection between the ray and the segment (since a segment contains an infinite number of points); to say it in other words, it does not exist a single intersection point. Do these special points in the green area intersect, and how many time? In general, the best way is to "shift" a little bit the point up or down, and this means that there's not a green area. So we can say that in such a situation, there's no intersection.

Points on the border delimiting the polygons are another trouble. For sure, we can say there exists (at least) one intersection point, and it is the same point we are checking. But back to our problem of determining if a point is inside or outside a polygon: a point on the border of the polygon is inside or outside it? We can say both, and be not wrong whatever we've choosen. Points on the borders or on the vertices are a problem; the algorithms here exposed want not to be coherent in the answer when such points are considered.

Another way to know if the point in those triangles will cause intersection with the segment or not, is to compare the slopes of the segment and of a segment from point A (the one placed below) and the point we are considering. This method is sketched in the following pictures.

The ray (from the point towards right horizontally) intersects the red segment if the angle of the blue segment (with the x axis, of course!) is greater than the angle of the red segment (rightmost image); if it is less than, then the ray does not intersect (leftmost image).

To find the angle, we compute the angular coefficient and then we should compute the inverse of the trigonometric function tangent (often arctan or atan in computer languages), which returns an angle in the range [-π/2, π/2]. To make things simpler we can "map" the range [-π/2, 0) to (π/2, π]. Simply, if the given angle A is negative, we compute A := 180 - A.

This way, we must compute 4 atan, that is a little bit time consuming and after all not necessary: we can use directly the angular coefficient, which is computed as

Where the subscripts refer to two distinct points on the line we want to know the angular coefficient of. We note that if the line is vertical, the x coordinates are the same and m is infinity. The atan of a very big number (infinity from a computer point of view) tends to ±π/2 (both indetify the same line).

The m of the horizontal line is of course 0.

Another important fact is that the angular coefficient (slope) of a line is negative when the angle between the ray ("semi-line" or half-line) and the x axis is greater than 90 degrees. For a little bit more than 90 degrees, it is huge and negative, and it becomes smaller ("towards" 0) as the angle approaches 180 degrees (i.e. the horizontal line).

So let us take another look at the rightmost image. The m of the blue segment with one end-point PB, is "closer" to 0 than the one of the red segment. Id est, m(blue) > m(red) so that we know that the ray will intersect the red segment.

If we take another look at leftmost image, we note that here m(blue) < m(red)... and in fact the ray from the point PA does not intersect the red segment.

This almost conclude our simplistic trip into the problem. We can now write an algorithm to check if a ray walking from a given point towards right intersect a segment (defined by two given end points). (If A is a point, Ax is its x coordinate)

  ray_intersects_segment:
     P : the point from which the ray starts
     A : the end-point of the segment with the smallest y coordinate
         (A must be "below" B)
     B : the end-point of the segment with the greatest y coordinate
         (B must be "above" A)
  % move the point upward of an ε to avoid the ray passing on the
  % vertex
  if Py = Ay or Py = By then
    Py ← Py + ε
  end if
  % check if P is in the hatched areas
  if Py > Ay or Py < By then 
    return false
  else if Px > max(Ax, Bx) then 
    return false
  else
    % if it is not in the hatched areas, check if it is in the
    % green area
    if Px < min(Ax, Bx) then
      return true
    else
      if Ax ≠ Bx then
        m_red ← (By - Ay)/(Bx - Ax)
      else
        m_red ← ∞
      end if
      if Ax ≠ Px then
        m_blue ← (Py - Ay)/(Px - Ax)
      else
        m_blue ← ∞
      end if
      if m_blue > m_red then
        return true
      else
        return false
      end if
    end if
  end if

Note: the angular coefficient of the vertical line is ±∞; if we consider always +∞ we could have problem: imagine the vertical red segment and a point to its left, so that in theory the blue segment should have a greater angular coefficient, according to our speech, i.e. we expect that m(blue) > ∞. This is of course false (it would have been true with -∞ instead of ∞) and apparently the algorithm fails. But this is not so, since the case of a point on the left of the segment is caught by the "greenish area" case, and similarly it's true for a point on the right. The only cases where we should worry about are when the point is on the segment. But I've already written we don't care for incoherent results in these special cases.

When you read incoherent you don't have to think random, nor that the fact can create serious problems. As said before, it depends on how we consider a point on the border of the polygon: is it inside or outside? If we answer once "inside" and once "outside" (for two different points), it does not happen nothing special.

At the end, we've just added the check to avoid a divion by zero that can never occur, save in few special cases. It is left to the reader to understand and enumerate them.

Why the "even-odd rule" works?

The reason why the things work this way can be intuitively understood, without too many efforts. And in practice it is also a proof. I will explain it the following way (like a tale since otherwise it would be extremely simple...!).

On the world you are on, Flatworld, two coutries exist: one is the land of Nobody, and the other is the land of Someone. You normally travel by feet during the night, don't pay attention to the landscape in the dark... You like this way: to be in some place, without knowing exactly where!

Now you are sitting on a chair-like rock, sun-bathing and chewing a licorice root, with your eyes closed. Suddenly you hear horse footsteps approaching you and you half-open your eyes, just to check if it's a bandit. You see a harmless stranger.

«Good morning sir» he says.

«Good morning» you reply. He looks like you imagine rich people look. «Can I help you?» you ask.

«Indeed, you can help me. I must know in which land we are. It's a very important issue, I can pay to know it».

You have enough money, but life is so strange, and it's better not to miss an opportunity. So your mind start thinking deeply. You currently don't know where you are... you have a smart idea, but you need an horse and time.

«Sir» you say, «I can help you, but you must lend me your horse, and I will say where we are in two days». First the stranger looks at you as you were a bandit. Then he talk:

«You clearly don't know your position exactly like me. But you have an idea to know it. I will lend my horse, and pay you, even if you explain me how you want to proceed to know our position»

So you explain: «You know that Flatworld has an outer border that can't be crossed: there's a wall, a very tall wall, which separates Flatland from the Void Flat Space. Maybe you ignore, but I can assure you it is true, that the land near the outer border is the land of Nobody. So my idea is the following: I will ride always in the same direction until I will reach the outer border. I will count the borders I will cross. If we are in the land of Nobody, the first border I will cross will let me enter the land of Someone, the second border I will cross will let me enter the land of Nobody and so on. So that, when I will be in the land of Nobody near the outer border which can't be crossed, if I've crossed an even number of borders, then the land we are on must be the land of Nobody. Otherwise, it is the land of Someone». While you explained it, you drew some lines in the sand.

«I see it. Go and ride until the extreme borders of Flatworld... but wait! Is this the proof that this planet is really flat?»

«No sir. It just suggests that we could be prisoners in a bigger spherical world!» You hit the horse that jumps forward poiting to the extreme borders of Flatworld running with the flat wind.

«Wait!» the stranger shouts. But you are far and can't hear him. «If it is so, it's not important to know where I am!» he says in a sad weak whisper when he realized you are gone.

There's no need to explain more, except maybe the fact that a border separates always a land from a different land... internal borders do not exist, or they can be distinguished from "country borders".


Home page