

 This page uses content from Wikipedia. The original article was at Point_in_polygon. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



Given a point and a polygon, check if the point is inside or outside the polygon using the ray-casting algorithm.

A pseudocode can be simply:

 count ← 0
 foreach side in polygon:
   if ray_intersects_segment(P,side) then
     count ← count + 1
 if is_odd(count) then
   return inside
 else
   return outside

Where the function ray_intersects_segment return true if the horizontal ray starting from the point P intersects the side (segment), false otherwise.

An intuitive explanation of why it works is that every time we cross
a border, we change "country" (inside-outside, or outside-inside), but
the last "country" we land on is surely outside (since the inside of the polygon is finite, while the ray continues towards infinity). So, if we crossed an odd number of borders we were surely inside, otherwise we were outside; we can follow the ray backward to see it better: starting from outside, only an odd number of crossing can give an inside: outside-inside, outside-inside-outside-inside, and so on (the - represents the crossing of a border).

So the main part of the algorithm is how we determine if a ray intersects a segment. The following text explain one of the possible ways.

 
Looking at the image on the right, we can easily be convinced of the fact that rays starting from points in the hatched area (like P1 and P2) surely do not intersect the segment AB. We also can easily see that rays starting from points in the greenish area surely intersect the segment AB (like point P3).

So the problematic points are those inside the white area (the box delimited by the points A and B), like P4.

 
 
Let us take into account a segment AB (the point A having y coordinate always smaller than B's y coordinate, i.e. point A is always below point B) and a point P. Let us use the cumbersome notation PAX to denote the angle between segment AP and AX, where X is always a point on the horizontal line passing by A with x coordinate bigger than the maximum between the x coordinate of A and the x coordinate of B. As explained graphically by the figures on the right, if PAX is greater than the angle BAX, then the ray starting from P intersects the segment AB. (In the images, the ray starting from PA does not intersect the segment, while the ray starting from PB in the second picture, intersects the segment).

Points on the boundary or "on" a vertex are someway special and through this approach we do not obtain coherent results. They could be treated apart, but it is not necessary to do so.

An algorithm for the previous speech could be (if P is a point, Px 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)
 if Py = Ay or Py = By then
   Py ← Py + ε
 end if
 if Py < Ay or Py > By then 
   return false
 else if Px >= max(Ax, Bx) then 
   return false
 else
   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

(To avoid the "ray on vertex" problem, the point is moved upward of a small quantity   ε.)




1000 PUBLIC NUMERIC x,y
1010 LET x=1
1020 LET y=2
1030 !
1040 DEF isLeft2(L(,),p()) = -SGN(  (L(1,x)-L(2,x))*(p(y)-L(2,y)) - (p(x)-L(2,x))*(L(1,y)-L(2,y)))
1050 !
1060 FUNCTION inpolygon(p1(,),p2())
1070    LET k=UBOUND(p1,1)+1
1080    DIM send (1 TO 2,2)
1090    LET wn=0
1100    FOR n=1 TO UBOUND(p1,1)
1110       LET index=MOD(n, k)
1120       LET nextindex=MOD(n+1, k)
1130       IF nextindex=0 THEN LET nextindex=1
1140       LET send(1,x)=p1(index,x)
1150       LET send(2,x)=p1(nextindex,x)
1160       LET send(1,y)=p1(index,y)
1170       LET send(2,y)=p1(nextindex,y)
1180       IF p1(index,y)<=p2(y) THEN
1190          IF p1(nextindex,y)>p2(y) THEN
1200             IF isleft2(send,p2)>=0 THEN !'=
1210                LET wn=wn+1
1220             END IF
1230          END IF
1240       ELSE
1250          IF p1(nextindex,y)<=p2(y) THEN
1260             IF isleft2(send,p2)<=0 THEN !'=
1270                LET wn=wn-1
1280             END IF
1290          END IF
1300       END IF
1310    NEXT n
1320    LET inpolygon = wn
1330 END FUNCTION
1340 !
1350 DIM type(1 TO 2)
1360 !
1370 DIM square(4,2)
1380 MAT READ square
1390 DATA 0,0,10,0,10,10,0,10
1400 !
1410 DIM hole(4,2)
1420 MAT READ hole
1430 DATA 2.5,2.5,7.5,2.5,7.5,7.5,2.5,7.5
1440 !
1450 DIM strange(8,2)
1460 MAT READ strange
1470 DATA 0,0,2.5,2.5,0,10,2.5,7.5,7.5,7.5,10,10,10,0,2.5,2.5
1480 !
1490 DIM exagon(6,2)  
1500 MAT READ exagon
1510 DATA 3,0,7,0,10,5,7,10,3,10,0,5
1520 !
1530 ! printouts
1540 FOR z=1 TO 4
1550    SELECT CASE z
1560    CASE 1
1570       PRINT "squared"
1580       PRINT "(5,5)  ";TAB(12);
1590       MAT READ type
1600       DATA 5,5
1610       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1620       MAT READ type
1630       DATA 5,8
1640       PRINT "(5,8)  ";TAB(12);
1650       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1660       PRINT "(-10,5)  ";TAB(12);
1670       MAT READ type
1680       DATA -10,5
1690       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1700       Print "(0,5)  ";Tab(12);
1710       MAT READ type
1720       DATA 0,5
1730       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1740       Print "(10,5)  ";Tab(12);
1750       MAT READ type
1760       DATA 10,5
1770       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1780       PRINT "(8,5)  ";TAB(12);
1790       MAT READ type
1800       DATA 8,5
1810       IF inpolygon(square,Type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1820       PRINT "(10,10)  ";TAB(12);
1830       MAT READ type
1840       DATA 10,10
1850       IF inpolygon(square,Type) <> 0 THEN PRINT "in" ELSE PRINT "out"
1860       PRINT
1870    CASE 2
1880       PRINT "squared hole"
1890       PRINT "(5,5)  ";TAB(12);
1900       MAT READ type
1910       DATA 5,5
1920       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
1930       Print "(5,8)  ";Tab(12);
1940       MAT READ type
1950       DATA 5,8
1960       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
1970       PRINT "(-10,5)  ";TAB(12);
1980       MAT READ type
1990       DATA -10,5
2000       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2010       PRINT "(0,5)  ";TAB(12);
2020       MAT READ type
2030       DATA 0,5
2040       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2050       PRINT "(10,5)  ";TAB(12);
2060       MAT READ type
2070       DATA 10,5
2080       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2090       PRINT "(8,5)  ";TAB(12);
2100       MAT READ type
2110       DATA 8,5
2120       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2130       PRINT "(10,10)  ";TAB(12);
2140       MAT READ type
2150       DATA 10,10
2160       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2170       PRINT
2180    CASE 3
2190       PRINT "strange"
2200       PRINT "(5,5)  ";TAB(12);
2210       MAT READ type
2220       DATA 5,5
2230       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2240       PRINT "(5,8)  ";TAB(12);
2250       MAT READ type
2260       DATA 5,8
2270       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2280       PRINT "(-10,5)  ";TAB(12);
2290       MAT READ type
2300       DATA -10,5
2310       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2320       PRINT "(0,5)  ";TAB(12);
2330       MAT READ type
2340       DATA 0,5
2350       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2360       PRINT "(10,5)  ";TAB(12);
2370       MAT READ type
2380       DATA 10,5
2390       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2400       PRINT "(8,5)  ";TAB(12);
2410       MAT READ type
2420       DATA 8,5
2430       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2440       PRINT "(10,10)  ";TAB(12);
2450       MAT READ type
2460       DATA 10,10
2470       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2480       PRINT
2490    CASE 4
2500       PRINT "exagon"
2510       PRINT "(5,5)  ";TAB(12);
2520       MAT READ type
2530       DATA 5,5
2540       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2550       PRINT "(5,8)  ";TAB(12);
2560       MAT READ type
2570       DATA 5,8
2580       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2590       PRINT "(-10,5)  ";TAB(12);
2600       MAT READ type
2610       DATA -10,5
2620       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2630       PRINT "(0,5)  ";TAB(12);
2640       MAT READ type
2650       DATA 0,5
2660       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2670       PRINT "(10,5)  ";TAB(12);
2680       MAT READ type
2690       DATA 10,5
2700       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2710       PRINT "(8,5)  ";TAB(12);
2720       MAT READ type
2730       DATA 8,5
2740       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2750       PRINT "(10,10)  ";TAB(12);
2760       MAT READ type
2770       DATA 10,10
2780       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"
2790       PRINT
2800    END SELECT
2810 NEXT z
2820 END
 