Ray-casting algorithm

Ray-casting algorithm
You are encouraged to solve this task according to the task description, using any language you may know.
 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   ε.)

package Polygons is    type Point is record      X, Y : Float;   end record;   type Point_List is array (Positive range <>) of Point;   subtype Segment is Point_List (1 .. 2);   type Polygon is array (Positive range <>) of Segment;    function Create_Polygon (List : Point_List) return Polygon;    function Is_Inside (Who : Point; Where : Polygon) return Boolean; end Polygons;

package body Polygons is   EPSILON : constant := 0.00001;    function Ray_Intersects_Segment     (Who   : Point;      Where : Segment)      return  Boolean   is      The_Point        : Point   := Who;      Above            : Point;      Below            : Point;      M_Red            : Float;      Red_Is_Infinity  : Boolean := False;      M_Blue           : Float;      Blue_Is_Infinity : Boolean := False;   begin      if Where (1).Y < Where (2).Y then         Above := Where (2);         Below := Where (1);      else         Above := Where (1);         Below := Where (2);      end if;      if The_Point.Y = Above.Y or The_Point.Y = Below.Y then         The_Point.Y := The_Point.Y + EPSILON;      end if;      if The_Point.Y < Below.Y or The_Point.Y > Above.Y then         return False;      elsif The_Point.X > Above.X and The_Point.X > Below.X then         return False;      elsif The_Point.X < Above.X and The_Point.X < Below.X then         return True;      else         if Above.X /= Below.X then            M_Red := (Above.Y - Below.Y) / (Above.X - Below.X);         else            Red_Is_Infinity := True;         end if;         if Below.X /= The_Point.X then            M_Blue := (The_Point.Y - Below.Y) / (The_Point.X - Below.X);         else            Blue_Is_Infinity := True;         end if;         if Blue_Is_Infinity then            return True;         elsif Red_Is_Infinity then            return False;         elsif M_Blue >= M_Red then            return True;         else            return False;         end if;      end if;   end Ray_Intersects_Segment;    function Create_Polygon (List : Point_List) return Polygon is      Result : Polygon (List'Range);      Side   : Segment;   begin      for I in List'Range loop         Side (1) := List (I);         if I = List'Last then            Side (2) := List (List'First);         else            Side (2) := List (I + 1);         end if;         Result (I) := Side;      end loop;      return Result;   end Create_Polygon;    function Is_Inside (Who : Point; Where : Polygon) return Boolean is      Count : Natural := 0;   begin      for Side in Where'Range loop         if Ray_Intersects_Segment (Who, Where (Side)) then            Count := Count + 1;         end if;      end loop;      if Count mod 2 = 0 then         return False;      else         return True;      end if;   end Is_Inside; end Polygons;

Example use:

with Ada.Text_IO;with Polygons;procedure Main is   package Float_IO is new Ada.Text_IO.Float_IO (Float);   Test_Points : Polygons.Point_List :=     ((  5.0,  5.0),      (  5.0,  8.0),      (-10.0,  5.0),      (  0.0,  5.0),      ( 10.0,  5.0),      (  8.0,  5.0),      ( 10.0, 10.0));   Square      : Polygons.Polygon    :=     ((( 0.0,  0.0), (10.0,  0.0)),      ((10.0,  0.0), (10.0, 10.0)),      ((10.0, 10.0), ( 0.0, 10.0)),      (( 0.0, 10.0), ( 0.0,  0.0)));   Square_Hole : Polygons.Polygon    :=     ((( 0.0,  0.0), (10.0,  0.0)),      ((10.0,  0.0), (10.0, 10.0)),      ((10.0, 10.0), ( 0.0, 10.0)),      (( 0.0, 10.0), ( 0.0,  0.0)),      (( 2.5,  2.5), ( 7.5,  2.5)),      (( 7.5,  2.5), ( 7.5,  7.5)),      (( 7.5,  7.5), ( 2.5,  7.5)),      (( 2.5,  7.5), ( 2.5,  2.5)));   Strange     : Polygons.Polygon    :=     ((( 0.0,  0.0), ( 2.5,  2.5)),      (( 2.5,  2.5), ( 0.0, 10.0)),      (( 0.0, 10.0), ( 2.5,  7.5)),      (( 2.5,  7.5), ( 7.5,  7.5)),      (( 7.5,  7.5), (10.0, 10.0)),      ((10.0, 10.0), (10.0,  0.0)),      ((10.0,  0.0), ( 2.5,  2.5)));   Exagon      : Polygons.Polygon    :=     ((( 3.0,  0.0), ( 7.0,  0.0)),      (( 7.0,  0.0), (10.0,  5.0)),      ((10.0,  5.0), ( 7.0, 10.0)),      (( 7.0, 10.0), ( 3.0, 10.0)),      (( 3.0, 10.0), ( 0.0,  5.0)),      (( 0.0,  5.0), ( 3.0,  0.0)));begin   Ada.Text_IO.Put_Line ("Testing Square:");   for Point in Test_Points'Range loop      Ada.Text_IO.Put ("Point(");      Float_IO.Put (Test_Points (Point).X, 0, 0, 0);      Ada.Text_IO.Put (",");      Float_IO.Put (Test_Points (Point).Y, 0, 0, 0);      Ada.Text_IO.Put        ("): " &         Boolean'Image (Polygons.Is_Inside (Test_Points (Point), Square)));      Ada.Text_IO.New_Line;   end loop;   Ada.Text_IO.New_Line;   Ada.Text_IO.Put_Line ("Testing Square_Hole:");   for Point in Test_Points'Range loop      Ada.Text_IO.Put ("Point(");      Float_IO.Put (Test_Points (Point).X, 0, 0, 0);      Ada.Text_IO.Put (",");      Float_IO.Put (Test_Points (Point).Y, 0, 0, 0);      Ada.Text_IO.Put        ("): " &         Boolean'Image            (Polygons.Is_Inside (Test_Points (Point), Square_Hole)));      Ada.Text_IO.New_Line;   end loop;   Ada.Text_IO.New_Line;   Ada.Text_IO.Put_Line ("Testing Strange:");   for Point in Test_Points'Range loop      Ada.Text_IO.Put ("Point(");      Float_IO.Put (Test_Points (Point).X, 0, 0, 0);      Ada.Text_IO.Put (",");      Float_IO.Put (Test_Points (Point).Y, 0, 0, 0);      Ada.Text_IO.Put        ("): " &         Boolean'Image (Polygons.Is_Inside (Test_Points (Point), Strange)));      Ada.Text_IO.New_Line;   end loop;   Ada.Text_IO.New_Line;   Ada.Text_IO.Put_Line ("Testing Exagon:");   for Point in Test_Points'Range loop      Ada.Text_IO.Put ("Point(");      Float_IO.Put (Test_Points (Point).X, 0, 0, 0);      Ada.Text_IO.Put (",");      Float_IO.Put (Test_Points (Point).Y, 0, 0, 0);      Ada.Text_IO.Put        ("): " &         Boolean'Image (Polygons.Is_Inside (Test_Points (Point), Exagon)));      Ada.Text_IO.New_Line;   end loop;end Main;

Output:

Testing Square:
Point(5.0,5.0): TRUE
Point(5.0,8.0): TRUE
Point(-10.0,5.0): FALSE
Point(0.0,5.0): FALSE
Point(10.0,5.0): TRUE
Point(8.0,5.0): TRUE
Point(10.0,10.0): FALSE

Testing Square_Hole:
Point(5.0,5.0): FALSE
Point(5.0,8.0): TRUE
Point(-10.0,5.0): FALSE
Point(0.0,5.0): FALSE
Point(10.0,5.0): TRUE
Point(8.0,5.0): TRUE
Point(10.0,10.0): FALSE

Testing Strange:
Point(5.0,5.0): TRUE
Point(5.0,8.0): FALSE
Point(-10.0,5.0): FALSE
Point(0.0,5.0): FALSE
Point(10.0,5.0): TRUE
Point(8.0,5.0): TRUE
Point(10.0,10.0): FALSE

Testing Exagon:
Point(5.0,5.0): TRUE
Point(5.0,8.0): TRUE
Point(-10.0,5.0): FALSE
Point(0.0,5.0): FALSE
Point(10.0,5.0): TRUE
Point(8.0,5.0): TRUE
Point(10.0,10.0): FALSE

ANSI Standard BASIC

Translation of: FreeBASIC
1000 PUBLIC NUMERIC x,y1010 LET x=11020 LET y=21030 !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)+11080    DIM send (1 TO 2,2)1090    LET wn=01100    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=11140       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) THEN1190          IF p1(nextindex,y)>p2(y) THEN1200             IF isleft2(send,p2)>=0 THEN !'=1210                LET wn=wn+11220             END IF1230          END IF1240       ELSE1250          IF p1(nextindex,y)<=p2(y) THEN1260             IF isleft2(send,p2)<=0 THEN !'=1270                LET wn=wn-11280             END IF1290          END IF1300       END IF1310    NEXT n1320    LET inpolygon = wn1330 END FUNCTION1340 !1350 DIM type(1 TO 2)1360 !1370 DIM square(4,2)1380 MAT READ square1390 DATA 0,0,10,0,10,10,0,101400 !1410 DIM hole(4,2)1420 MAT READ hole1430 DATA 2.5,2.5,7.5,2.5,7.5,7.5,2.5,7.51440 !1450 DIM strange(8,2)1460 MAT READ strange1470 DATA 0,0,2.5,2.5,0,10,2.5,7.5,7.5,7.5,10,10,10,0,2.5,2.51480 !1490 DIM exagon(6,2)  1500 MAT READ exagon1510 DATA 3,0,7,0,10,5,7,10,3,10,0,51520 !1530 ! printouts1540 FOR z=1 TO 41550    SELECT CASE z1560    CASE 11570       PRINT "squared"1580       PRINT "(5,5)  ";TAB(12);1590       MAT READ type1600       DATA 5,51610       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"1620       MAT READ type1630       DATA 5,81640       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 type1680       DATA -10,51690       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"1700       Print "(0,5)  ";Tab(12);1710       MAT READ type1720       DATA 0,51730       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"1740       Print "(10,5)  ";Tab(12);1750       MAT READ type1760       DATA 10,51770       IF inpolygon(square,type) <> 0 THEN PRINT "in" ELSE PRINT "out"1780       PRINT "(8,5)  ";TAB(12);1790       MAT READ type1800       DATA 8,51810       IF inpolygon(square,Type) <> 0 THEN PRINT "in" ELSE PRINT "out"1820       PRINT "(10,10)  ";TAB(12);1830       MAT READ type1840       DATA 10,101850       IF inpolygon(square,Type) <> 0 THEN PRINT "in" ELSE PRINT "out"1860       PRINT1870    CASE 21880       PRINT "squared hole"1890       PRINT "(5,5)  ";TAB(12);1900       MAT READ type1910       DATA 5,51920       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 type1950       DATA 5,81960       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 type1990       DATA -10,52000       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 type2030       DATA 0,52040       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 type2070       DATA 10,52080       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 type2110       DATA 8,52120       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 type2150       DATA 10,102160       IF NOT inpolygon(hole,Type)<>0 AND inpolygon(square,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2170       PRINT2180    CASE 32190       PRINT "strange"2200       PRINT "(5,5)  ";TAB(12);2210       MAT READ type2220       DATA 5,52230       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2240       PRINT "(5,8)  ";TAB(12);2250       MAT READ type2260       DATA 5,82270       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2280       PRINT "(-10,5)  ";TAB(12);2290       MAT READ type2300       DATA -10,52310       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2320       PRINT "(0,5)  ";TAB(12);2330       MAT READ type2340       DATA 0,52350       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2360       PRINT "(10,5)  ";TAB(12);2370       MAT READ type2380       DATA 10,52390       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2400       PRINT "(8,5)  ";TAB(12);2410       MAT READ type2420       DATA 8,52430       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2440       PRINT "(10,10)  ";TAB(12);2450       MAT READ type2460       DATA 10,102470       IF inpolygon(strange,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2480       PRINT2490    CASE 42500       PRINT "exagon"2510       PRINT "(5,5)  ";TAB(12);2520       MAT READ type2530       DATA 5,52540       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2550       PRINT "(5,8)  ";TAB(12);2560       MAT READ type2570       DATA 5,82580       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2590       PRINT "(-10,5)  ";TAB(12);2600       MAT READ type2610       DATA -10,52620       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2630       PRINT "(0,5)  ";TAB(12);2640       MAT READ type2650       DATA 0,52660       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2670       PRINT "(10,5)  ";TAB(12);2680       MAT READ type2690       DATA 10,52700       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2710       PRINT "(8,5)  ";TAB(12);2720       MAT READ type2730       DATA 8,52740       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2750       PRINT "(10,10)  ";TAB(12);2760       MAT READ type2770       DATA 10,102780       IF inpolygon(exagon,Type)<>0 THEN PRINT "in" ELSE PRINT "out"2790       PRINT2800    END SELECT2810 NEXT z2820 END

AutoHotkey

Works with: AutoHotkey L
Points :=[{x:  5.0, y: 5.0}		, {x:  5.0, y: 8.0}		, {x:-10.0, y: 5.0}		, {x:  0.0, y: 5.0}		, {x: 10.0, y: 5.0}		, {x:  8.0, y: 5.0}		, {x: 10.0, y:10.0}]Square :=[{x: 0.0, y: 0.0}, {x:10.0, y: 0.0}		, {x:10.0, y: 0.0}, {x:10.0, y:10.0}		, {x:10.0, y:10.0}, {x: 0.0, y:10.0}		, {x: 0.0, y:10.0}, {x: 0.0, y: 0.0}]Sq_Hole:=[{x: 0.0, y: 0.0}, {x:10.0, y: 0.0}		, {x:10.0, y: 0.0}, {x:10.0, y:10.0}		, {x:10.0, y:10.0}, {x: 0.0, y:10.0}		, {x: 0.0, y:10.0}, {x: 0.0, y: 0.0}		, {x: 2.5, y: 2.5}, {x: 7.5, y: 2.5}		, {x: 7.5, y: 2.5}, {x: 7.5, y: 7.5}		, {x: 7.5, y: 7.5}, {x: 2.5, y: 7.5}		, {x: 2.5, y: 7.5}, {x: 2.5, y: 2.5}]Strange:=[{x: 0.0, y: 0.0}, {x: 2.5, y: 2.5}		, {x: 2.5, y: 2.5}, {x: 0.0, y:10.0}		, {x: 0.0, y:10.0}, {x: 2.5, y: 7.5}		, {x: 2.5, y: 7.5}, {x: 7.5, y: 7.5}		, {x: 7.5, y: 7.5}, {x:10.0, y:10.0}		, {x:10.0, y:10.0}, {x:10.0, y: 0.0}		, {x:10.0, y: 0.0}, {x: 2.5, y: 2.5}]Exagon :=[{x: 3.0, y: 0.0}, {x: 7.0, y: 0.0}		, {x: 7.0, y: 0.0}, {x:10.0, y: 5.0}		, {x:10.0, y: 5.0}, {x: 7.0, y:10.0}		, {x: 7.0, y:10.0}, {x: 3.0, y:10.0}		, {x: 3.0, y:10.0}, {x: 0.0, y: 5.0}		, {x: 0.0, y: 5.0}, {x: 3.0, y: 0.0}]Polygons := {"Square":Square, "Sq_Hole":Sq_Hole, "Strange":Strange, "Exagon":Exagon}For j, Poly in Polygons	For i, Point in Points		If (point_in_polygon(Point,Poly))			s.= j " does contain point " i "n"		Else			s.= j " doesn't contain point " i "n"Msgbox %s% point_in_polygon(Point,Poly) {	n:=Poly.MaxIndex()	count:=0	loop, %n% {		if (ray_intersects_segment(Point,Poly[A_Index],Poly[mod(A_Index,n)+1])) {			count++		}	}	if (mod(count,2)) { ; true = inside, false = outside		return true		; P is in the polygon	} else {		return false	; P isn't in the polygon	}} ray_intersects_segment(P,A,B) {	;P = the point from which the ray starts	;A = the end-point of the segment with the smallest y coordinate	;B = the end-point of the segment with the greatest y coordinate	if (A.y > B.y) {		temp:=A		A:=B		B:=temp	}	if (P.y = A.y or P.y = B.y) {		P.y += 0.000001	}	if (P.y < A.y or P.y > B.y) {		return false	} else if (P.x > A.x && P.x > B.x) {		return false	} else {		if (P.x < A.x && P.x < B.x) {			return true		} else {			if (A.x != B.x) {				m_red := (B.y - A.y)/(B.x - A.x)			} else {				m_red := "inf"			}			if (A.x != P.x) {				m_blue := (P.y - A.y)/(P.x - A.x)			} else {				m_blue := "inf"			}			if (m_blue >= m_red) {				return true			} else {				return false			}		}	}}
Output:
---------------------------
Ray-casting_algorithm.ahkl
---------------------------
Exagon does contain point 1
Exagon does contain point 2
Exagon doesn't contain point 3
Exagon doesn't contain point 4
Exagon does contain point 5
Exagon does contain point 6
Exagon doesn't contain point 7
Sq_Hole doesn't contain point 1
Sq_Hole does contain point 2
Sq_Hole doesn't contain point 3
Sq_Hole doesn't contain point 4
Sq_Hole does contain point 5
Sq_Hole does contain point 6
Sq_Hole doesn't contain point 7
Square does contain point 1
Square does contain point 2
Square doesn't contain point 3
Square doesn't contain point 4
Square does contain point 5
Square does contain point 6
Square doesn't contain point 7
Strange does contain point 1
Strange doesn't contain point 2
Strange doesn't contain point 3
Strange doesn't contain point 4
Strange does contain point 5
Strange does contain point 6
Strange doesn't contain point 7

---------------------------
OK
---------------------------

C

#include <stdio.h>#include <stdlib.h>#include <math.h> typedef struct { double x, y; } vec;typedef struct { int n; vec* v; } polygon_t, *polygon; #define BIN_V(op, xx, yy) vec v##op(vec a,vec b){vec c;c.x=xx;c.y=yy;return c;}#define BIN_S(op, r) double v##op(vec a, vec b){ return r; }BIN_V(sub, a.x - b.x, a.y - b.y);BIN_V(add, a.x + b.x, a.y + b.y);BIN_S(dot, a.x * b.x + a.y * b.y);BIN_S(cross, a.x * b.y - a.y * b.x); /* return a + s * b */vec vmadd(vec a, double s, vec b){	vec c;	c.x = a.x + s * b.x;	c.y = a.y + s * b.y;	return c;} /* check if x0->x1 edge crosses y0->y1 edge. dx = x1 - x0, dy = y1 - y0, then   solve  x0 + a * dx == y0 + b * dy with a, b in real   cross both sides with dx, then: (remember, cross product is a scalar)	x0 X dx = y0 X dx + b * (dy X dx)   similarly,	x0 X dy + a * (dx X dy) == y0 X dy   there is an intersection iff 0 <= a <= 1 and 0 <= b <= 1    returns: 1 for intersect, -1 for not, 0 for hard to say (if the intersect   point is too close to y0 or y1)*/int intersect(vec x0, vec x1, vec y0, vec y1, double tol, vec *sect){	vec dx = vsub(x1, x0), dy = vsub(y1, y0);	double d = vcross(dy, dx), a;	if (!d) return 0; /* edges are parallel */ 	a = (vcross(x0, dx) - vcross(y0, dx)) / d;	if (sect)		*sect = vmadd(y0, a, dy); 	if (a < -tol || a > 1 + tol) return -1;	if (a < tol || a > 1 - tol) return 0; 	a = (vcross(x0, dy) - vcross(y0, dy)) / d;	if (a < 0 || a > 1) return -1; 	return 1;} /* distance between x and nearest point on y0->y1 segment.  if the point   lies outside the segment, returns infinity */double dist(vec x, vec y0, vec y1, double tol){	vec dy = vsub(y1, y0);	vec x1, s;	int r; 	x1.x = x.x + dy.y; x1.y = x.y - dy.x;	r = intersect(x, x1, y0, y1, tol, &s);	if (r == -1) return HUGE_VAL;	s = vsub(s, x);	return sqrt(vdot(s, s));} #define for_v(i, z, p) for(i = 0, z = p->v; i < p->n; i++, z++)/* returns 1 for inside, -1 for outside, 0 for on edge */int inside(vec v, polygon p, double tol){	/* should assert p->n > 1 */	int i, k, crosses, intersectResult;	vec *pv;	double min_x, max_x, min_y, max_y; 	for (i = 0; i < p->n; i++) {		k = (i + 1) % p->n;		min_x = dist(v, p->v[i], p->v[k], tol);		if (min_x < tol) return 0;	} 	min_x = max_x = p->v[0].x;	min_y = max_y = p->v[1].y; 	/* calculate extent of polygon */	for_v(i, pv, p) {		if (pv->x > max_x) max_x = pv->x;		if (pv->x < min_x) min_x = pv->x;		if (pv->y > max_y) max_y = pv->y;		if (pv->y < min_y) min_y = pv->y;	}	if (v.x < min_x || v.x > max_x || v.y < min_y || v.y > max_y)		return -1; 	max_x -= min_x; max_x *= 2;	max_y -= min_y; max_y *= 2;	max_x += max_y; 	vec e;	while (1) {		crosses = 0;		/* pick a rand point far enough to be outside polygon */		e.x = v.x + (1 + rand() / (RAND_MAX + 1.)) * max_x;		e.y = v.y + (1 + rand() / (RAND_MAX + 1.)) * max_x; 		for (i = 0; i < p->n; i++) {			k = (i + 1) % p->n;			intersectResult = intersect(v, e, p->v[i], p->v[k], tol, 0); 			/* picked a bad point, ray got too close to vertex.			   re-pick */			if (!intersectResult) break; 			if (intersectResult == 1) crosses++;		}		if (i == p->n) break;	}	return (crosses & 1) ? 1 : -1;} int main(){	vec vsq[] = {	{0,0}, {10,0}, {10,10}, {0,10},			{2.5,2.5}, {7.5,0.1}, {7.5,7.5}, {2.5,7.5}}; 	polygon_t sq = { 4, vsq }, /* outer square */		sq_hole = { 8, vsq }; /* outer and inner square, ie hole */ 	vec c = { 10, 5 }; /* on edge */	vec d = { 5, 5 }; 	printf("%d\n", inside(c, &sq, 1e-10));	printf("%d\n", inside(c, &sq_hole, 1e-10)); 	printf("%d\n", inside(d, &sq, 1e-10));	/* in */	printf("%d\n", inside(d, &sq_hole, 1e-10));  /* out (in the hole) */ 	return 0;}

C++

Works with: C++ version 11
Translation of: D
#include <algorithm>#include <cstdlib>#include <iomanip>#include <iostream>#include <limits> using namespace std; const double epsilon = numeric_limits<float>().epsilon();const numeric_limits<double> DOUBLE;const double MIN = DOUBLE.min();const double MAX = DOUBLE.max(); struct Point { const double x, y; }; struct Edge {    const Point a, b;     bool operator()(const Point& p) const    {        if (a.y > b.y) return Edge{ b, a }(p);        if (p.y == a.y || p.y == b.y) return operator()({ p.x, p.y + epsilon });        if (p.y > b.y || p.y < a.y || p.x > max(a.x, b.x)) return false;        if (p.x < min(a.x, b.x)) return true;        auto blue = abs(a.x - p.x) > MIN ? (p.y - a.y) / (p.x - a.x) : MAX;        auto red = abs(a.x - b.x) > MIN ? (b.y - a.y) / (b.x - a.x) : MAX;        return blue >= red;    }}; struct Figure {    const string  name;    const initializer_list<Edge> edges;     bool contains(const Point& p) const    {        auto c = 0;        for (auto e : edges) if (e(p)) c++;        return c % 2 != 0;    }     template<unsigned char W = 3>    void check(const initializer_list<Point>& points, ostream& os) const    {        os << "Is point inside figure " << name <<  '?' << endl;        for (auto p : points)            os << "  (" << setw(W) << p.x << ',' << setw(W) << p.y << "): " << boolalpha << contains(p) << endl;        os << endl;    }}; int main(){    const initializer_list<Point> points =  { { 5.0, 5.0}, {5.0, 8.0}, {-10.0, 5.0}, {0.0, 5.0}, {10.0, 5.0}, {8.0, 5.0}, {10.0, 10.0} };    const Figure square = { "Square",        {  {{0.0, 0.0}, {10.0, 0.0}}, {{10.0, 0.0}, {10.0, 10.0}}, {{10.0, 10.0}, {0.0, 10.0}}, {{0.0, 10.0}, {0.0, 0.0}} }    };     const Figure square_hole = { "Square hole",        {  {{0.0, 0.0}, {10.0, 0.0}}, {{10.0, 0.0}, {10.0, 10.0}}, {{10.0, 10.0}, {0.0, 10.0}}, {{0.0, 10.0}, {0.0, 0.0}},           {{2.5, 2.5}, {7.5, 2.5}}, {{7.5, 2.5}, {7.5, 7.5}}, {{7.5, 7.5}, {2.5, 7.5}}, {{2.5, 7.5}, {2.5, 2.5}}        }    };     const Figure strange = { "Strange",        {  {{0.0, 0.0}, {2.5, 2.5}}, {{2.5, 2.5}, {0.0, 10.0}}, {{0.0, 10.0}, {2.5, 7.5}}, {{2.5, 7.5}, {7.5, 7.5}},           {{7.5, 7.5}, {10.0, 10.0}}, {{10.0, 10.0}, {10.0, 0.0}}, {{10.0, 0}, {2.5, 2.5}}        }    };     const Figure exagon = { "Exagon",        {  {{3.0, 0.0}, {7.0, 0.0}}, {{7.0, 0.0}, {10.0, 5.0}}, {{10.0, 5.0}, {7.0, 10.0}}, {{7.0, 10.0}, {3.0, 10.0}},           {{3.0, 10.0}, {0.0, 5.0}}, {{0.0, 5.0}, {3.0, 0.0}}        }    };     for(auto f : {square, square_hole, strange, exagon})        f.check(points, cout);     return EXIT_SUCCESS;}
Output:

As D.

CoffeeScript

Takes a polygon as a list of points joining segments, and creates segments between them.

  Point = (@x,@y) ->   pointInPoly = (point,poly) ->    segments = for pointA, index in poly                 pointB = poly[(index + 1) % poly.length]                 [pointA,pointB]    intesected = (segment for segment in segments when rayIntesectsSegment(point,segment))    intesected.length % 2 != 0   rayIntesectsSegment = (p,segment) ->    [p1,p2] = segment    [a,b] = if p1.y < p2.y              [p1,p2]            else              [p2,p1]    if p.y == b.y || p.y == a.y      p.y += Number.MIN_VALUE     if p.y > b.y || p.y < a.y      false    else if p.x > a.x && p.x > b.x      false    else if p.x < a.x && p.x < b.x      true    else      mAB = (b.y - a.y) / (b.x - a.x)      mAP = (p.y - a.y) / (p.x - a.x)      mAP > mAB

Common Lisp

Points are represented as cons cells whose car is an x value and whose cdr is a y value. A line segment is a cons cell of two points. A polygon is a list of line segments.

(defun point-in-polygon (point polygon)  (do ((in-p nil)) ((endp polygon) in-p)    (when (ray-intersects-segment point (pop polygon))      (setf in-p (not in-p))))) (defun ray-intersects-segment (point segment &optional (epsilon .001))  (destructuring-bind (px . py) point    (destructuring-bind ((ax . ay) . (bx . by)) segment      (when (< ay by)        (rotatef ay by)        (rotatef ax bx))      (when (or (= py ay) (= py by))        (incf py epsilon))      (cond       ;; point is above, below, or to the right of the rectangle       ;; determined by segment; ray does not intesect the segment.       ((or (> px (max ax bx)) (> py (max ay by)) (< py (min ay by)))        nil)       ;; point is to left of the rectangle; ray intersects segment       ((< px (min ax bx))        t)       ;; point is within the rectangle...       (t (let ((m-red (if (= ax bx) nil                         (/ (- by ay) (- bx ax))))                (m-blue (if (= px ax) nil                          (/ (- py ay) (- px ax)))))            (cond             ((null m-blue) t)             ((null m-red) nil)             (t (>= m-blue m-red)))))))))

Testing code

(defparameter *points*  #((0 . 0) (10 . 0) (10 . 10) (0 . 10)      (2.5 . 2.5) (7.5 . 2.5) (7.5 . 7.5) (2.5 . 7.5)     (0 . 5) (10 . 5) (3 . 0) (7 . 0)    (7 . 10) (3 . 10))) (defun create-polygon (indices &optional (points *points*))  (loop for (a b) on indices by 'cddr        collecting (cons (aref points (1- a))                         (aref points (1- b))))) (defun square ()  (create-polygon '(1 2 2 3 3 4 4 1))) (defun square-hole ()  (create-polygon '(1 2 2 3 3 4 4 1 5 6 6 7 7 8 8 5))) (defun strange ()  (create-polygon '(1 5 5 4 4 8 8 7 7 3 3 2 2 5))) (defun exagon ()  (create-polygon '(11 12 12 10 10 13 13 14 14 9 9 11))) (defparameter *test-points*  #((5 . 5) (5 . 8) (-10 . 5) (0 . 5)    (10 . 5) (8 . 5) (10 . 10))) (defun test-pip ()  (dolist (shape '(square square-hole strange exagon))    (print shape)    (loop with polygon = (funcall shape)          for test-point across *test-points*          do (format t "~&~w ~:[outside~;inside ~]."                     test-point                     (point-in-polygon test-point polygon)))))

D

import std.stdio, std.math, std.algorithm; immutable struct Point { double x, y; }immutable struct Edge { Point a, b; }immutable struct Figure {    string name;    Edge[] edges;} bool contains(in Figure poly, in Point p) pure nothrow @safe @nogc {    static bool raySegI(in Point p, in Edge edge)    pure nothrow @safe @nogc {        enum double epsilon = 0.00001;        with (edge) {            if (a.y > b.y)                //swap(a, b); // if edge is mutable                return raySegI(p, Edge(b, a));            if (p.y == a.y || p.y == b.y)                //p.y += epsilon; // if p is mutable                return raySegI(Point(p.x, p.y + epsilon), edge);            if (p.y > b.y || p.y < a.y || p.x > max(a.x, b.x))                return false;            if (p.x < min(a.x, b.x))                return true;            immutable blue = (abs(a.x - p.x) > double.min_normal) ?                             ((p.y - a.y) / (p.x - a.x)) :                             double.max;            immutable red = (abs(a.x - b.x) > double.min_normal) ?                            ((b.y - a.y) / (b.x - a.x)) :                            double.max;            return blue >= red;        }    }     return poly.edges.count!(e => raySegI(p, e)) % 2;} void main() {    immutable Figure[] polys = [  {"Square", [    {{ 0.0,  0.0}, {10.0,  0.0}},  {{10.0,  0.0}, {10.0, 10.0}},    {{10.0, 10.0}, { 0.0, 10.0}},  {{ 0.0, 10.0}, { 0.0,  0.0}}]},  {"Square hole", [    {{ 0.0,  0.0}, {10.0,  0.0}},  {{10.0,  0.0}, {10.0, 10.0}},    {{10.0, 10.0}, { 0.0, 10.0}},  {{ 0.0, 10.0}, { 0.0,  0.0}},    {{ 2.5,  2.5}, { 7.5,  2.5}},  {{ 7.5,  2.5}, { 7.5,  7.5}},    {{ 7.5,  7.5}, { 2.5,  7.5}},  {{ 2.5,  7.5}, { 2.5,  2.5}}]},  {"Strange", [    {{ 0.0,  0.0}, { 2.5,  2.5}},  {{ 2.5,  2.5}, { 0.0, 10.0}},    {{ 0.0, 10.0}, { 2.5,  7.5}},  {{ 2.5,  7.5}, { 7.5,  7.5}},    {{ 7.5,  7.5}, {10.0, 10.0}},  {{10.0, 10.0}, {10.0,  0.0}},    {{10.0,  0},   { 2.5,  2.5}}]},  {"Exagon", [    {{ 3.0,  0.0}, { 7.0,  0.0}},  {{ 7.0,  0.0}, {10.0,  5.0}},    {{10.0,  5.0}, { 7.0, 10.0}},  {{ 7.0, 10.0}, { 3.0, 10.0}},    {{ 3.0, 10.0}, { 0.0,  5.0}},  {{ 0.0,  5.0}, { 3.0,  0.0}}]}];     immutable Point[] testPoints = [{ 5, 5}, {5, 8}, {-10,  5}, {0, 5},                                    {10, 5}, {8, 5}, { 10, 10}];     foreach (immutable poly; polys) {        writefln(Is point inside figure "%s"?, poly.name);        foreach (immutable p; testPoints)            writefln("  (%3s, %2s): %s", p.x, p.y, contains(poly, p));        writeln;    }}
Output:
Is point inside figure "Square"?
(  5,  5): true
(  5,  8): true
(-10,  5): false
(  0,  5): false
( 10,  5): true
(  8,  5): true
( 10, 10): false

Is point inside figure "Square hole"?
(  5,  5): false
(  5,  8): true
(-10,  5): false
(  0,  5): false
( 10,  5): true
(  8,  5): true
( 10, 10): false

Is point inside figure "Strange"?
(  5,  5): true
(  5,  8): false
(-10,  5): false
(  0,  5): false
( 10,  5): true
(  8,  5): true
( 10, 10): false

Is point inside figure "Exagon"?
(  5,  5): true
(  5,  8): true
(-10,  5): false
(  0,  5): false
( 10,  5): true
(  8,  5): true
( 10, 10): false

Factor

To test whether a ray intersects a line, we test that the starting point is between the endpoints in y value, and that it is to the left of the point on the segment with the same y value. Note that this implementation does not support polygons with horizontal edges.

USING: kernel prettyprint sequences arrays math math.vectors ;IN: raycasting : between ( a b x -- ? ) [ last ] [email protected] [ < ] curry [email protected] xor ; : lincomb ( a b x -- w )  3dup [ last ] [email protected]  [ - ] curry [email protected]  [ drop ] 2dip  neg 2dup + [ / ] curry [email protected]  [ [ v*n ] curry ] [email protected] bi*  v+ ;: leftof ( a b x -- ? ) dup [ lincomb ] dip [ first ] [email protected] > ; : ray ( a b x -- ? ) [ between ] [ leftof ] 3bi and ; : raycast ( poly x -- ? )  [ dup first suffix [ rest-slice ] [ but-last-slice ] bi ] dip  [ ray ] curry 2map  f [ xor ] reduce ;

Usage:

( scratchpad ) CONSTANT: square { { -2 -1 } { 1 -2 } { 2 1 } { -1 2 } }( scratchpad ) square { 0 0 } raycast .t( scratchpad ) square { 5 5 } raycast .f( scratchpad ) square { 2 0 } raycast .f

Fortran

Works with: Fortran version 95 and later

The following code uses the Points_Module defined here.

This module defines "polygons".

module Polygons  use Points_Module  implicit none   type polygon     type(point), dimension(:), allocatable :: points     integer, dimension(:), allocatable :: vertices  end type polygon contains   function create_polygon(pts, vt)    type(polygon) :: create_polygon    type(point), dimension(:), intent(in) :: pts    integer, dimension(:), intent(in) :: vt     integer :: np, nv     np = size(pts,1)    nv = size(vt,1)     allocate(create_polygon%points(np), create_polygon%vertices(nv))    create_polygon%points = pts    create_polygon%vertices = vt   end function create_polygon   subroutine free_polygon(pol)    type(polygon), intent(inout) :: pol     deallocate(pol%points, pol%vertices)   end subroutine free_polygon end module Polygons

The ray casting algorithm module:

module Ray_Casting_Algo  use Polygons  implicit none   real, parameter, private :: eps = 0.00001  private :: ray_intersects_seg contains   function ray_intersects_seg(p0, a0, b0) result(intersect)    type(point), intent(in) :: p0, a0, b0    logical :: intersect     type(point) :: a, b, p    real :: m_red, m_blue     p = p0    ! let variable "a" be the point with smallest y coordinate    if ( a0%y > b0%y ) then       b = a0       a = b0    else       a = a0       b = b0    end if     if ( (p%y == a%y) .or. (p%y == b%y) ) p%y = p%y + eps     intersect = .false.     if ( (p%y > b%y) .or. (p%y < a%y) ) return    if ( p%x > max(a%x, b%x) ) return     if ( p%x < min(a%x, b%x) ) then       intersect = .true.    else       if ( abs(a%x - b%x) > tiny(a%x) ) then          m_red = (b%y - a%y) / (b%x - a%x)       else          m_red = huge(m_red)       end if       if ( abs(a%x - p%x) > tiny(a%x) ) then          m_blue = (p%y - a%y) / (p%x - a%x)       else          m_blue = huge(m_blue)       end if       if ( m_blue >= m_red ) then          intersect = .true.       else          intersect = .false.       end if    end if   end function ray_intersects_seg   function point_is_inside(p, pol) result(inside)    logical :: inside    type(point), intent(in) :: p    type(polygon), intent(in) :: pol     integer :: i, cnt, pa, pb     cnt = 0    do i = lbound(pol%vertices,1), ubound(pol%vertices,1), 2       pa = pol%vertices(i)       pb = pol%vertices(i+1)       if ( ray_intersects_seg(p, pol%points(pa), pol%points(pb)) ) cnt = cnt + 1    end do     inside = .true.    if ( mod(cnt, 2) == 0 ) then       inside = .false.    end if   end function point_is_inside end module Ray_Casting_Algo

Testing

program Pointpoly  use Points_Module  use Ray_Casting_Algo  implicit none   character(len=16), dimension(4) :: names  type(polygon), dimension(4) :: polys  type(point), dimension(14) :: pts  type(point), dimension(7) :: p   integer :: i, j   pts = (/ point(0,0), point(10,0), point(10,10), point(0,10), &           point(2.5,2.5), point(7.5,2.5), point(7.5,7.5), point(2.5,7.5), &           point(0,5), point(10,5), &           point(3,0), point(7,0), point(7,10), point(3,10) /)   polys(1) = create_polygon(pts, (/ 1,2, 2,3, 3,4, 4,1 /) )  polys(2) = create_polygon(pts, (/ 1,2, 2,3, 3,4, 4,1, 5,6, 6,7, 7,8, 8,5 /) )  polys(3) = create_polygon(pts, (/ 1,5, 5,4, 4,8, 8,7, 7,3, 3,2, 2,5 /) )  polys(4) = create_polygon(pts, (/ 11,12, 12,10, 10,13, 13,14, 14,9, 9,11 /) )   names = (/ "square", "square hole", "strange", "hexagon" /)   p = (/ point(5,5), point(5, 8), point(-10, 5), point(0,5), point(10,5), &         point(8,5), point(10,10) /)   do j = 1, size(p)     do i = 1, size(polys)        write(*, "('point (',F8.2,',',F8.2,') is inside ',A,'? ', L)") &             p(j)%x, p(j)%y, names(i), point_is_inside(p(j), polys(i))     end do     print *, ""  end do   do i = 1, size(polys)     call free_polygon(polys(i))  end do end program Pointpoly

FreeBASIC

Inpolygon by Winding number method

Type Point  As Single x,yEnd Type Function inpolygon(p1() As Point,p2 As Point) As Integer  #Macro 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))  #EndMacro  Dim As Integer index,nextindex  Dim k As Integer=UBound(p1)+1  Dim send (1 To 2) As Point  Dim wn As Integer=0  For n As Integer=1 To UBound(p1)    index=n Mod k:nextindex=(n+1) Mod k    If nextindex=0 Then nextindex=1    send(1).x=p1(index).x:send(2).x=p1(nextindex).x    send(1).y=p1(index).y:send(2).y=p1(nextindex).y    If p1(index).y<=p2.y Then      If p1(nextindex).y>p2.y Then        If isleft2(send,p2)>=0 Then '=          wn=wn+1        End If      End If    Else      If p1(nextindex).y<=p2.y Then        If isleft2(send,p2)<=0 Then'=          wn=wn-1        End If      End If    End If  Next n  Return wnEnd Function  Dim As Point square(1 To 4)  ={(0,0),(10,0),(10,10),(0,10)} Dim As Point hole(1 To 4)    ={(2.5,2.5),(7.5,2.5),(7.5,7.5),(2.5,7.5)} Dim As Point strange(1 To 8) ={(0,0),(2.5,2.5),(0,10),(2.5,7.5),_                              (7.5,7.5),(10,10),(10,0),(2.5,2.5)} Dim As Point exagon(1 To 6)  ={(3,0),(7,0),(10,5),(7,10),(3,10),(0,5)} 'printoutsFor z As Integer=1 To 4  Select Case z    Case 1: Print "squared"      Print "(5,5)  " ;Tab(12);      If inpolygon(square(),Type<Point>(5,5)) Then Print "in" Else Print "out"      Print "(5,8)  " ;Tab(12);      If inpolygon(square(),Type<Point>(5,8)) Then Print "in" Else Print "out"      Print "(-10,5)  " ;Tab(12);      If inpolygon(square(),Type<Point>(-10,5)) Then Print "in" Else Print "out"      Print "(0,5)  " ;Tab(12);      If inpolygon(square(),Type<Point>(0,5)) Then Print "in" Else Print "out"      Print "(10,5)  " ;Tab(12);      If inpolygon(square(),Type<Point>(10,5)) Then Print "in" Else Print "out"      Print "(8,5)  " ;Tab(12);      If inpolygon(square(),Type<Point>(8,5)) Then Print "in" Else Print "out"      Print "(10,10)  " ;Tab(12);      If inpolygon(square(),Type<Point>(10,10)) Then Print "in" Else Print "out"      Print    Case 2:Print "squared hole"      Print "(5,5)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(5,5)) And inpolygon(square(),Type<Point>(5,5)) Then Print "in" Else Print "out"      Print "(5,8)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(5,8)) And inpolygon(square(),Type<Point>(5,8))Then Print "in" Else Print "out"      Print "(-10,5)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(-10,5))And inpolygon(square(),Type<Point>(-10,5)) Then Print "in" Else Print "out"      Print "(0,5)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(0,5))And inpolygon(square(),Type<Point>(0,5)) Then Print "in" Else Print "out"      Print "(10,5)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(10,5))And inpolygon(square(),Type<Point>(10,5)) Then Print "in" Else Print "out"      Print "(8,5)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(8,5))And inpolygon(square(),Type<Point>(8,5)) Then Print "in" Else Print "out"      Print "(10,10)  " ;Tab(12);      If Not inpolygon(hole(),Type<Point>(10,10))And inpolygon(square(),Type<Point>(10,10)) Then Print "in" Else Print "out"      Print    Case 3:Print "strange"      Print "(5,5)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(5,5)) Then Print "in" Else Print "out"      Print "(5,8)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(5,8)) Then Print "in" Else Print "out"      Print "(-10,5)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(-10,5)) Then Print "in" Else Print "out"      Print "(0,5)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(0,5)) Then Print "in" Else Print "out"      Print "(10,5)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(10,5)) Then Print "in" Else Print "out"      Print "(8,5)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(8,5)) Then Print "in" Else Print "out"      Print "(10,10)  " ;Tab(12);      If inpolygon(strange(),Type<Point>(10,10)) Then Print "in" Else Print "out"      Print    Case 4:Print "exagon"      Print "(5,5)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(5,5)) Then Print "in" Else Print "out"      Print "(5,8)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(5,8)) Then Print "in" Else Print "out"      Print "(-10,5)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(-10,5)) Then Print "in" Else Print "out"      Print "(0,5)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(0,5)) Then Print "in" Else Print "out"      Print "(10,5)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(10,5)) Then Print "in" Else Print "out"      Print "(8,5)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(8,5)) Then Print "in" Else Print "out"      Print "(10,10)  " ;Tab(12);      If inpolygon(exagon(),Type<Point>(10,10)) Then Print "in" Else Print "out"      Print  End SelectNext zSleep

Output:

squared
(5,5)      in
(5,8)      in
(-10,5)    out
(0,5)      out
(10,5)     in
(8,5)      in
(10,10)    out

squared hole
(5,5)      out
(5,8)      in
(-10,5)    out
(0,5)      out
(10,5)     in
(8,5)      in
(10,10)    out

strange
(5,5)      in
(5,8)      out
(-10,5)    out
(0,5)      out
(10,5)     in
(8,5)      in
(10,10)    out

exagon
(5,5)      in
(5,8)      in
(-10,5)    out
(0,5)      out
(10,5)     in
(8,5)      in
(10,10)    out

Go

The first solution given here follows the model of most other solutions on the page in defining a polygon as a list of segments. Unfortunately this representation does not require that the polygon is closed. Input to the ray-casting algorithm, as noted in the WP article though, is specified to be a closed polygon. The "strange" shape defined here is not a closed polygon and so gives incorrect results against some points. (Graphically it may appear closed but mathematically it needs an additional segment returning to the starting point.)

package main import (    "fmt"    "math") type xy struct {    x, y float64} type seg struct {    p1, p2 xy} type poly struct {    name  string    sides []seg} func inside(pt xy, pg poly) (i bool) {    for _, side := range pg.sides {        if rayIntersectsSegment(pt, side) {            i = !i        }    }    return} func rayIntersectsSegment(p xy, s seg) bool {    var a, b xy    if s.p1.y < s.p2.y {        a, b = s.p1, s.p2    } else {        a, b = s.p2, s.p1    }    for p.y == a.y || p.y == b.y {        p.y = math.Nextafter(p.y, math.Inf(1))    }    if p.y < a.y || p.y > b.y {        return false    }    if a.x > b.x {        if p.x > a.x {            return false        }        if p.x < b.x {            return true        }    } else {        if p.x > b.x {            return false        }        if p.x < a.x {            return true        }    }    return (p.y-a.y)/(p.x-a.x) >= (b.y-a.y)/(b.x-a.x)} var (    p1  = xy{0, 0}    p2  = xy{10, 0}    p3  = xy{10, 10}    p4  = xy{0, 10}    p5  = xy{2.5, 2.5}    p6  = xy{7.5, 2.5}    p7  = xy{7.5, 7.5}    p8  = xy{2.5, 7.5}    p9  = xy{0, 5}    p10 = xy{10, 5}    p11 = xy{3, 0}    p12 = xy{7, 0}    p13 = xy{7, 10}    p14 = xy{3, 10}) var tpg = []poly{    {"square", []seg{{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1}}},    {"square hole", []seg{{p1, p2}, {p2, p3}, {p3, p4}, {p4, p1},        {p5, p6}, {p6, p7}, {p7, p8}, {p8, p5}}},    {"strange", []seg{{p1, p5},        {p5, p4}, {p4, p8}, {p8, p7}, {p7, p3}, {p3, p2}, {p2, p5}}},    {"exagon", []seg{{p11, p12}, {p12, p10}, {p10, p13},        {p13, p14}, {p14, p9}, {p9, p11}}},} var tpt = []xy{    // test points common in other solutions on this page    {5, 5}, {5, 8}, {-10, 5}, {0, 5}, {10, 5}, {8, 5}, {10, 10},    // test points that show the problem with "strange"    {1, 2}, {2, 1},} func main() {    for _, pg := range tpg {        fmt.Printf("%s:\n", pg.name)        for _, pt := range tpt {            fmt.Println(pt, inside(pt, pg))        }    }}
Output:
square:
{5 5} true
{5 8} true
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} true
{2 1} true
square hole:
{5 5} false
{5 8} true
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} true
{2 1} true
strange:
{5 5} true
{5 8} false
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} true
{2 1} false
exagon:
{5 5} true
{5 8} true
{-10 5} false
{0 5} false
{10 5} true
{8 5} true
{10 10} false
{1 2} false
{2 1} false


Closed polygon solution

Here input is given as a list of N vertices defining N segments, where one segment extends from each vertex to the next, and one more extends from the last vertex to the first. In the case of the "strange" shape, this mathematically closes the polygon and allows the program to give correct results.

package main import (    "math"    "fmt") type xy struct {    x, y float64} type closedPoly struct {    name string    vert []xy} func inside(pt xy, pg closedPoly) bool {    if len(pg.vert) < 3 {        return false    }    in := rayIntersectsSegment(pt, pg.vert[len(pg.vert)-1], pg.vert[0])    for i := 1; i < len(pg.vert); i++ {        if rayIntersectsSegment(pt, pg.vert[i-1], pg.vert[i]) {            in = !in        }    }    return in} func rayIntersectsSegment(p, a, b xy) bool {    if a.y > b.y {        a, b = b, a    }    for p.y == a.y || p.y == b.y {        p.y = math.Nextafter(p.y, math.Inf(1))    }    if p.y < a.y || p.y > b.y {        return false    }    if a.x > b.x {        if p.x > a.x {            return false        }        if p.x < b.x {            return true        }    } else {        if p.x > b.x {            return false        }        if p.x < a.x {            return true        }    }    return (p.y-a.y)/(p.x-a.x) >= (b.y-a.y)/(b.x-a.x)} var tpg = []closedPoly{    {"square", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}}},    {"square hole", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}, {0, 0},        {2.5, 2.5}, {7.5, 2.5}, {7.5, 7.5}, {2.5, 7.5}, {2.5, 2.5}}},    {"strange", []xy{{0, 0}, {2.5, 2.5}, {0, 10}, {2.5, 7.5}, {7.5, 7.5},        {10, 10}, {10, 0}, {2.5, 2.5}}},    {"exagon", []xy{{3, 0}, {7, 0}, {10, 5}, {7, 10}, {3, 10}, {0, 5}}},} var tpt = []xy{{1, 2}, {2, 1}} func main() {    for _, pg := range tpg {        fmt.Printf("%s:\n", pg.name)        for _, pt := range tpt {            fmt.Println(pt, inside(pt, pg))        }    }}
Output:
square:
{1 2} true
{2 1} true
square hole:
{1 2} true
{2 1} true
strange:
{1 2} false
{2 1} false
exagon:
{1 2} false
{2 1} false


PNPoly algorithm

This solution replaces the rayIntersectsSegment function above with the expression from the popular PNPoly algorithm described at https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html. The expression is not only simpler but more accurate.

This solution is preferred over the two above.

package main import "fmt" type xy struct {    x, y float64} type closedPoly struct {    name string    vert []xy} func inside(pt xy, pg closedPoly) bool {    if len(pg.vert) < 3 {        return false    }    in := rayIntersectsSegment(pt, pg.vert[len(pg.vert)-1], pg.vert[0])    for i := 1; i < len(pg.vert); i++ {        if rayIntersectsSegment(pt, pg.vert[i-1], pg.vert[i]) {            in = !in        }    }    return in} func rayIntersectsSegment(p, a, b xy) bool {    return (a.y > p.y) != (b.y > p.y) &&        p.x < (b.x-a.x)*(p.y-a.y)/(b.y-a.y)+a.x} var tpg = []closedPoly{    {"square", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}}},    {"square hole", []xy{{0, 0}, {10, 0}, {10, 10}, {0, 10}, {0, 0},        {2.5, 2.5}, {7.5, 2.5}, {7.5, 7.5}, {2.5, 7.5}, {2.5, 2.5}}},    {"strange", []xy{{0, 0}, {2.5, 2.5}, {0, 10}, {2.5, 7.5}, {7.5, 7.5},        {10, 10}, {10, 0}, {2.5, 2.5}}},    {"exagon", []xy{{3, 0}, {7, 0}, {10, 5}, {7, 10}, {3, 10}, {0, 5}}},} var tpt = []xy{{1, 2}, {2, 1}} func main() {    for _, pg := range tpg {        fmt.Printf("%s:\n", pg.name)        for _, pt := range tpt {            fmt.Println(pt, inside(pt, pg))        }    }}

import Data.Ratio type Point = (Rational, Rational)type Polygon = [Point]data Line = Sloped {lineSlope, lineYIntercept :: Rational} |            Vert {lineXIntercept :: Rational} polygonSides :: Polygon -> [(Point, Point)]polygonSides poly@(p1 : ps) = zip poly $ps ++ [p1] intersects :: Point -> Line -> Bool{- @intersects (px, py) [email protected] is true if the ray {(x, py) | x ≥ px}intersects l. -}intersects (px, _) (Vert xint) = px <= xintintersects (px, py) (Sloped m b) | m < 0 = py <= m * px + b | otherwise = py >= m * px + b onLine :: Point -> Line -> Bool{- Is the point on the line? -}onLine (px, _) (Vert xint) = px == xintonLine (px, py) (Sloped m b) = py == m * px + b carrier :: (Point, Point) -> Line{- Finds the line containing the given line segment. -}carrier ((ax, ay), (bx, by)) | ax == bx = Vert ax | otherwise = Sloped slope yint where slope = (ay - by) / (ax - bx) yint = ay - slope * ax between :: Ord a => a -> a -> a -> Boolbetween x a b | a > b = b <= x && x <= a | otherwise = a <= x && x <= b inPolygon :: Point -> Polygon -> BoolinPolygon p@(px, py) = f 0 . polygonSides where f n [] = odd n f n (side : sides) | far = f n sides | onSegment = True | rayIntersects = f (n + 1) sides | otherwise = f n sides where far = not$ between py ay by                onSegment | ay == by  = between px ax bx                          | otherwise = p onLine line                rayIntersects =                    intersects p line &&                    (py /= ay || by < py) &&                    (py /= by || ay < py)                ((ax, ay), (bx, by)) = side                line = carrier side

J

NB.*crossPnP v point in closed polygon, crossing numberNB.  bool=. points crossPnP polygoncrossPnP=: 4 : 0"2  'X Y'=. |:x  'x0 y0 x1 y1'=. |:2 ,/\^:(2={:@[email protected]]) y  p1=. ((y0<:/Y)*. y1>/Y) +. (y0>/Y)*. y1<:/Y  p2=. (x0-/X) < (x0-x1) * (y0-/Y) % (y0 - y1)  2|+/ p1*.p2)

Sample data:

SQUAREV=:          0   0   , 10  0   , 10  10  ,: 0   10SQUAREV=: SQUAREV, 2.5 2.5 , 7.5 0.1 , 7.5 7.5 ,: 2.5 7.5 ESAV=: 3 0 , 7 0 , 10 5 , 7 10 , 3 10 ,: 0 5 ESA=:        (0 1,1 2,2 3,3 4,4 5,:5 0) , .{ ESAVSQUARE=:     (0 1,1 2,2 3,:3 0)         , .{ SQUAREVSQUAREHOLE=: (0 1,1 2,2 3,3 0,4 5,5 6,6 7,:7 4) , .{ SQUAREVSTRANGE=:    (0 4,4 3,3 7,7 6,6 2,2 1,1 5,:5 0) , .{ SQUAREV POINTS=: 5 5,5 8,2 2,0 0,10 10,2.5 2.5,0.01 5,2.2 7.4,0 5,10 5,:_4 10

Testing:

   (<POINTS) crossPnP every ESA;SQUARE;SQUAREHOLE;STRANGE1 1 1 0 0 1 1 1 0 1 01 1 1 0 0 1 1 1 0 1 00 1 1 0 0 1 1 1 0 1 01 0 0 0 0 0 0 1 0 1 0

Java

import static java.lang.Math.*; public class RayCasting {     static boolean intersects(int[] A, int[] B, double[] P) {        if (A[1] > B[1])            return intersects(B, A, P);         if (P[1] == A[1] || P[1] == B[1])            P[1] += 0.0001;         if (P[1] > B[1] || P[1] < A[1] || P[0] >= max(A[0], B[0]))            return false;         if (P[0] < min(A[0], B[0]))            return true;         double red = (P[1] - A[1]) / (double) (P[0] - A[0]);        double blue = (B[1] - A[1]) / (double) (B[0] - A[0]);        return red >= blue;    }     static boolean contains(int[][] shape, double[] pnt) {        boolean inside = false;        int len = shape.length;        for (int i = 0; i < len; i++) {            if (intersects(shape[i], shape[(i + 1) % len], pnt))                inside = !inside;        }        return inside;    }     public static void main(String[] a) {        double[][] testPoints = {{10, 10}, {10, 16}, {-20, 10}, {0, 10},        {20, 10}, {16, 10}, {20, 20}};         for (int[][] shape : shapes) {            for (double[] pnt : testPoints)                System.out.printf("%7s ", contains(shape, pnt));            System.out.println();        }    }     final static int[][] square = {{0, 0}, {20, 0}, {20, 20}, {0, 20}};     final static int[][] squareHole = {{0, 0}, {20, 0}, {20, 20}, {0, 20},    {5, 5}, {15, 5}, {15, 15}, {5, 15}};     final static int[][] strange = {{0, 0}, {5, 5}, {0, 20}, {5, 15}, {15, 15},    {20, 20}, {20, 0}};     final static int[][] hexagon = {{6, 0}, {14, 0}, {20, 10}, {14, 20},    {6, 20}, {0, 10}};     final static int[][][] shapes = {square, squareHole, strange, hexagon};}
   true    true   false   false    true    true   false
false    true   false   false    true    true   false
true   false   false   false    true    true   false
true    true   false   false    true    true   false 

JavaScript

 /** * @return {boolean} true if (lng, lat) is in bounds */function contains(bounds, lat, lng) {    //https://rosettacode.org/wiki/Ray-casting_algorithm    var count = 0;    for (var b = 0; b < bounds.length; b++) {        var vertex1 = bounds[b];        var vertex2 = bounds[(b + 1) % bounds.length];        if (west(vertex1, vertex2, lng, lat))            ++count;    }    return count % 2;     /**     * @return {boolean} true if (x,y) is west of the line segment connecting A and B     */    function west(A, B, x, y) {        if (A.y <= B.y) {            if (y <= A.y || y > B.y ||                x >= A.x && x >= B.x) {                return false;            } else if (x < A.x && x < B.x) {                return true;            } else {                return (y - A.y) / (x - A.x) > (B.y - A.y) / (B.x - A.x);            }        } else {            return west(B, A, x, y);        }    }} var square = {name: 'square', bounds: [{x: 0, y: 0}, {x: 20, y: 0}, {x: 20, y: 20}, {x: 0, y: 20}]};var squareHole = {    name: 'squareHole',    bounds: [{x: 0, y: 0}, {x: 20, y: 0}, {x: 20, y: 20}, {x: 0, y: 20}, {x: 5, y: 5}, {x: 15, y: 5}, {x: 15, y: 15}, {x: 5, y: 15}]};var strange = {    name: 'strange',    bounds: [{x: 0, y: 0}, {x: 5, y: 5}, {x: 0, y: 20}, {x: 5, y: 15}, {x: 15, y: 15}, {x: 20, y: 20}, {x: 20, y: 0}]};var hexagon = {    name: 'hexagon',    bounds: [{x: 6, y: 0}, {x: 14, y: 0}, {x: 20, y: 10}, {x: 14, y: 20}, {x: 6, y: 20}, {x: 0, y: 10}]}; var shapes = [square, squareHole, strange, hexagon];var testPoints = [{lng: 10, lat: 10}, {lng: 10, lat: 16}, {lng: -20, lat: 10},    {lng: 0, lat: 10}, {lng: 20, lat: 10}, {lng: 16, lat: 10}, {lng: 20, lat: 20}]; for (var s = 0; s < shapes.length; s++) {    var shape = shapes[s];    for (var tp = 0; tp < testPoints.length; tp++) {        var testPoint = testPoints[tp];        console.log(JSON.stringify(testPoint) + '\tin ' + shape.name + '\t' + contains(shape.bounds, testPoint.lat, testPoint.lng));    }}

Julia

Works with: Julia version 0.6
Translation of: Python

Module:

module RayCastings export Point struct Point{T}    x::T    y::TendBase.show(io::IO, p::Point) = print(io, "($(p.x),$(p.y))") const Edge = Tuple{Point{T}, Point{T}} where TBase.show(io::IO, e::Edge) = print(io, "$(e[1]) ∘-∘$(e[2])") function rayintersectseg(p::Point{T}, edge::Edge{T}) where T    a, b = edge    if a.y > b.y        a, b = b, a    end    if p.y ∈ (a.y, b.y)        p = Point(p.x, p.y + eps(p.y))    end     rst = false    if (p.y > b.y || p.y < a.y) || (p.x > max(a.x, b.x))        return false    end     if p.x < min(a.x, b.x)        rst = true    else        mred = (b.y - a.y) / (b.x - a.x)        mblu = (p.y - a.y) / (p.x - a.x)        rst = mblu ≥ mred    end     return rstend isinside(poly::Vector{Tuple{Point{T}, Point{T}}}, p::Point{T}) where T =    isodd(count(edge -> rayintersectseg(p, edge), poly)) connect(a::Point{T}, b::Point{T}...) where T =    [(a, b) for (a, b) in zip(vcat(a, b...), vcat(b..., a))] end  # module RayCastings

Main:

let A = Point(0.0, 0.0),    B = Point(0.0, 10.0),    C = Point(10.0, 10.0),    D = Point(10.0, 0.0),    E = Point(2.5, 2.5),    F = Point(2.5, 7.5),    G = Point(7.5, 7.5),    H = Point(7.5, 2.5),    I = Point(3.0, 0.0),    J = Point(7.0, 0.0),    K = Point(10.0, 5.0),    L = Point(7.0, 10.0),    M = Point(3.0, 10.0),    N = Point(0.0, 5.0),    testpts = (Point(5.0, 5.0), Point(5.0, 8.0), Point(-10.0, 5.0), Point(0.0, 5.0),        Point(10.0, 5.0), Point(8.0, 5.0), Point(10.0, 10.0))     square = RayCastings.connect(A, B, C, D)    square_withhole = vcat(square, RayCastings.connect(E, F, G, H))    strange = RayCastings.connect(A, E, B, F, G, C, D, E)    exagon = RayCastings.connect(I, J, K, L, M, N)     println("\n# TESTING WHETHER POINTS ARE WITHIN POLYGONS")    for poly in (square, square_withhole, strange, exagon)        println("\nEdges: \n - ", join(poly, "\n - "))        println("Inside/outside:")        for p in testpts            @printf(" - %-12s is %s\n", p, RayCastings.isinside(poly, p) ? "inside" : "outside")        end    endend
Output:
# TESTING WHETHER POINTS ARE WITHIN POLYGONS

Edges:
- (0.0, 0.0) ∘-∘ (0.0, 10.0)
- (0.0, 10.0) ∘-∘ (10.0, 10.0)
- (10.0, 10.0) ∘-∘ (10.0, 0.0)
- (10.0, 0.0) ∘-∘ (0.0, 0.0)
Inside/outside:
- (5.0, 5.0)   is inside
- (5.0, 8.0)   is inside
- (-10.0, 5.0) is outside
- (0.0, 5.0)   is outside
- (10.0, 5.0)  is inside
- (8.0, 5.0)   is inside
- (10.0, 10.0) is outside

Edges:
- (0.0, 0.0) ∘-∘ (0.0, 10.0)
- (0.0, 10.0) ∘-∘ (10.0, 10.0)
- (10.0, 10.0) ∘-∘ (10.0, 0.0)
- (10.0, 0.0) ∘-∘ (0.0, 0.0)
- (2.5, 2.5) ∘-∘ (2.5, 7.5)
- (2.5, 7.5) ∘-∘ (7.5, 7.5)
- (7.5, 7.5) ∘-∘ (7.5, 2.5)
- (7.5, 2.5) ∘-∘ (2.5, 2.5)
Inside/outside:
- (5.0, 5.0)   is outside
- (5.0, 8.0)   is inside
- (-10.0, 5.0) is outside
- (0.0, 5.0)   is outside
- (10.0, 5.0)  is inside
- (8.0, 5.0)   is inside
- (10.0, 10.0) is outside

Edges:
- (0.0, 0.0) ∘-∘ (2.5, 2.5)
- (2.5, 2.5) ∘-∘ (0.0, 10.0)
- (0.0, 10.0) ∘-∘ (2.5, 7.5)
- (2.5, 7.5) ∘-∘ (7.5, 7.5)
- (7.5, 7.5) ∘-∘ (10.0, 10.0)
- (10.0, 10.0) ∘-∘ (10.0, 0.0)
- (10.0, 0.0) ∘-∘ (2.5, 2.5)
- (2.5, 2.5) ∘-∘ (0.0, 0.0)
Inside/outside:
- (5.0, 5.0)   is inside
- (5.0, 8.0)   is outside
- (-10.0, 5.0) is outside
- (0.0, 5.0)   is outside
- (10.0, 5.0)  is inside
- (8.0, 5.0)   is inside
- (10.0, 10.0) is outside

Edges:
- (3.0, 0.0) ∘-∘ (7.0, 0.0)
- (7.0, 0.0) ∘-∘ (10.0, 5.0)
- (10.0, 5.0) ∘-∘ (7.0, 10.0)
- (7.0, 10.0) ∘-∘ (3.0, 10.0)
- (3.0, 10.0) ∘-∘ (0.0, 5.0)
- (0.0, 5.0) ∘-∘ (3.0, 0.0)
Inside/outside:
- (5.0, 5.0)   is inside
- (5.0, 8.0)   is inside
- (-10.0, 5.0) is outside
- (0.0, 5.0)   is outside
- (10.0, 5.0)  is inside
- (8.0, 5.0)   is inside
- (10.0, 10.0) is outside

Kotlin

Translation of: D
import java.lang.Double.MAX_VALUEimport java.lang.Double.MIN_VALUEimport java.lang.Math.abs data class Point(val x: Double, val y: Double) data class Edge(val s: Point, val e: Point) {    operator fun invoke(p: Point) : Boolean = when {        s.y > e.y -> Edge(e, s).invoke(p)        p.y == s.y || p.y == e.y -> invoke(Point(p.x, p.y + epsilon))        p.y > e.y || p.y < s.y || p.x > Math.max(s.x, e.x) -> false        p.x < Math.min(s.x, e.x) -> true        else -> {            val blue = if (abs(s.x - p.x) > MIN_VALUE) (p.y - s.y) / (p.x - s.x) else MAX_VALUE            val red = if (abs(s.x - e.x) > MIN_VALUE) (e.y - s.y) / (e.x - s.x) else MAX_VALUE            blue >= red        }    }     val epsilon = 0.00001} class Figure(val name: String, val edges: Array<Edge>) {    operator fun contains(p: Point) = edges.count({ it(p) }) % 2 != 0} object Ray_casting {    fun check(figures : Array<Figure>, points : List<Point>) {        println("points: " + points)        figures.forEach { f ->            println("figure: " + f.name)            f.edges.forEach { println("        " + it) }            println("result: " + (points.map { it in f }))        }    }}
fun main(args: Array<String>) {    val figures = arrayOf(Figure("Square", arrayOf(Edge(Point(0.0, 0.0), Point(10.0, 0.0)), Edge(Point(10.0, 0.0), Point(10.0, 10.0)),            Edge(Point(10.0, 10.0), Point(0.0, 10.0)),Edge(Point(0.0, 10.0), Point(0.0, 0.0)))),    Figure("Square hole", arrayOf(Edge(Point(0.0, 0.0), Point(10.0, 0.0)), Edge(Point(10.0, 0.0), Point(10.0, 10.0)),            Edge(Point(10.0, 10.0), Point(0.0, 10.0)), Edge(Point(0.0, 10.0), Point(0.0, 0.0)), Edge(Point(2.5, 2.5), Point(7.5, 2.5)),            Edge(Point(7.5, 2.5), Point(7.5, 7.5)),Edge(Point(7.5, 7.5), Point(2.5, 7.5)), Edge(Point(2.5, 7.5), Point(2.5, 2.5)))),    Figure("Strange", arrayOf(Edge(Point(0.0, 0.0), Point(2.5, 2.5)), Edge(Point(2.5, 2.5), Point(0.0, 10.0)),            Edge(Point(0.0, 10.0), Point(2.5, 7.5)), Edge(Point(2.5, 7.5), Point(7.5, 7.5)), Edge(Point(7.5, 7.5), Point(10.0, 10.0)),            Edge(Point(10.0, 10.0), Point(10.0, 0.0)), Edge(Point(10.0, 0.0), Point(2.5, 2.5)))),    Figure("Exagon", arrayOf(Edge(Point(3.0, 0.0), Point(7.0, 0.0)), Edge(Point(7.0, 0.0), Point(10.0, 5.0)), Edge(Point(10.0, 5.0), Point(7.0, 10.0)),            Edge(Point(7.0, 10.0), Point(3.0, 10.0)), Edge(Point(3.0, 10.0), Point(0.0, 5.0)), Edge(Point(0.0, 5.0), Point(3.0, 0.0)))))     val points = listOf(Point(5.0, 5.0), Point(5.0, 8.0), Point(-10.0, 5.0), Point(0.0, 5.0),            Point(10.0, 5.0), Point(8.0, 5.0), Point(10.0, 10.0))     Ray_casting.check(figures, points)}
Output:
points: [Point(x=5.0, y=5.0), Point(x=5.0, y=8.0), Point(x=-10.0, y=5.0), Point(x=0.0, y=5.0), Point(x=10.0, y=5.0), Point(x=8.0, y=5.0), Point(x=10.0, y=10.0)]
figure: Square
Edge(s=Point(x=0.0, y=0.0), e=Point(x=10.0, y=0.0))
Edge(s=Point(x=10.0, y=0.0), e=Point(x=10.0, y=10.0))
Edge(s=Point(x=10.0, y=10.0), e=Point(x=0.0, y=10.0))
Edge(s=Point(x=0.0, y=10.0), e=Point(x=0.0, y=0.0))
result: [true, true, false, false, true, true, false]
figure: Square hole
Edge(s=Point(x=0.0, y=0.0), e=Point(x=10.0, y=0.0))
Edge(s=Point(x=10.0, y=0.0), e=Point(x=10.0, y=10.0))
Edge(s=Point(x=10.0, y=10.0), e=Point(x=0.0, y=10.0))
Edge(s=Point(x=0.0, y=10.0), e=Point(x=0.0, y=0.0))
Edge(s=Point(x=2.5, y=2.5), e=Point(x=7.5, y=2.5))
Edge(s=Point(x=7.5, y=2.5), e=Point(x=7.5, y=7.5))
Edge(s=Point(x=7.5, y=7.5), e=Point(x=2.5, y=7.5))
Edge(s=Point(x=2.5, y=7.5), e=Point(x=2.5, y=2.5))
result: [false, true, false, false, true, true, false]
figure: Strange
Edge(s=Point(x=0.0, y=0.0), e=Point(x=2.5, y=2.5))
Edge(s=Point(x=2.5, y=2.5), e=Point(x=0.0, y=10.0))
Edge(s=Point(x=0.0, y=10.0), e=Point(x=2.5, y=7.5))
Edge(s=Point(x=2.5, y=7.5), e=Point(x=7.5, y=7.5))
Edge(s=Point(x=7.5, y=7.5), e=Point(x=10.0, y=10.0))
Edge(s=Point(x=10.0, y=10.0), e=Point(x=10.0, y=0.0))
Edge(s=Point(x=10.0, y=0.0), e=Point(x=2.5, y=2.5))
result: [true, false, false, false, true, true, false]
figure: Exagon
Edge(s=Point(x=3.0, y=0.0), e=Point(x=7.0, y=0.0))
Edge(s=Point(x=7.0, y=0.0), e=Point(x=10.0, y=5.0))
Edge(s=Point(x=10.0, y=5.0), e=Point(x=7.0, y=10.0))
Edge(s=Point(x=7.0, y=10.0), e=Point(x=3.0, y=10.0))
Edge(s=Point(x=3.0, y=10.0), e=Point(x=0.0, y=5.0))
Edge(s=Point(x=0.0, y=5.0), e=Point(x=3.0, y=0.0))
result: [true, true, false, false, true, true, false]


Liberty BASIC

Translated from C code at: http://alienryderflex.com/polygon/

Displays interactively on-screen.

NoMainWinGlobal sw, sh, verts sw = 640 :   sh = 480WindowWidth  = sw+8 : WindowHeight = sh+31UpperLeftX = (DisplayWidth -sw)/2UpperLeftY = (DisplayHeight-sh)/2Open"Ray Casting Algorithm" For Graphics_nf_nsb As #g#g "Down; TrapClose [halt]"h$= "#g" Dim xp(15),yp(15)#g "when leftButtonDown [halt];when mouseMove checkPoint"#g "when rightButtonDown [Repeat]" [Repeat] #g "Cls;Fill 32 160 255; Color white;BackColor 32 160 255" #g "Place 5 460;\L-click to exit" #g "Place 485 460;\R-click for new polygon" 'generate polygon from random points numPoints = rand(4,15) verts = numPoints For i = 0 To numPoints-1 xp(i) = rand(20,620) yp(i) = rand(40,420) Next Call drawPoly h$, verts, "white"    #g "Flush"    Wait [halt]Close #gEnd 'Point In Polygon FunctionFunction pnp(x, y, numSides)    j= numSides-1: oddNodes = 0    For i = 0 To numSides-1        If ((yp(i)<y) And (yp(j)>=y)) Or ((yp(j)<y) And (yp(i)>=y)) Then            f1 = y - yp(i):f2 = yp(j) - yp(i): f3 = xp(j) - xp(i)            If (xp(i) + f1 / f2 * f3) < x Then oddNodes = 1 - oddNodes        End If        j = i    Next    pnp = oddNodesEnd Function 'draw the polygonSub drawPoly h$, verts, colour$    #h$, "Color ";colour$    j = verts-1    For i = 0 To verts-1        #h$"Line ";xp(j);" ";yp(j);" ";xp(i);" ";yp(i) j = i NextEnd Sub 'change message and color of polygonSub checkPoint h$, x, y    If pnp(x,y,verts) Then        #h$"Color 32 160 255;BackColor 32 160 255" #h$ "Place 5 0;BoxFilled 150 20;Color white"        #h$"Place 7 15;\Mouse In Polygon" Call drawPoly h$, verts, "red"    Else        #h$"Color 32 160 255;BackColor 32 160 255" #h$ "Place 5 0;BoxFilled 150 20;Color white"        #h$"Place 7 15;\Mouse Not In Polygon" Call drawPoly h$, verts, "white"    End IfEnd Sub Function rand(loNum,hiNum)    rand = Int(Rnd(0)*(hiNum-loNum+1)+loNum)End Function

OCaml

Translation of: C
type point = { x:float; y:float } type polygon = {  vertices: point array;  edges: (int * int) list;} let p x y = { x=x; y=y } let square_v = [|  (p 0. 0.); (p 10. 0.); (p 10. 10.); (p 0. 10.);  (p 2.5 2.5); (p 7.5 0.1); (p 7.5 7.5); (p 2.5 7.5)|] let esa_v = [|  (p 3. 0.); (p 7. 0.); (p 10. 5.); (p 7. 10.); (p 3. 10.); (p 0. 5.)|] let esa = {  vertices = esa_v;  edges = [ (0,1); (1,2); (2,3); (3,4); (4,5); (5,0) ]} let square = {  vertices = square_v;  edges = [ (0,1); (1,2); (2,3); (3,0) ]} let squarehole = {  vertices = square_v;  edges = [ (0,1); (1,2); (2,3); (3,0); (4,5); (5,6); (6,7); (7,4) ]} let strange = {  vertices = square_v;  edges = [ (0,4); (4,3); (3,7); (7,6); (6,2); (2,1); (1,5); (5,0) ]}  let min_y ~a ~b = if a.y > b.y then (b) else (a) let coeff_ang ~pa ~pb = (pb.y -. pa.y) /. (pb.x -. pa.x) let huge_val = infinity let hseg_intersect_seg ~s ~a ~b =  let _eps =    if s.y = (max a.y b.y) ||       s.y = (min a.y b.y) then 0.00001 else 0.0  in  if  (s.y +. _eps) > (max a.y b.y) ||      (s.y +. _eps) < (min a.y b.y) ||       s.x > (max a.x b.x) then (false)  else if s.x <= (min a.x b.x) then (true)  else    let ca = if a.x <> b.x then (coeff_ang a b) else (huge_val) in    let my = min_y ~a ~b in    let cp = if (s.x -. my.x) <> 0.0 then (coeff_ang my s) else (huge_val) in    (cp >= ca);;  let point_is_inside ~poly ~pt =  let cross = ref 0 in  List.iter (fun (a,b) ->    if hseg_intersect_seg pt             poly.vertices.(a)             poly.vertices.(b)    then incr cross  ) poly.edges;  ( (!cross mod 2) <> 0);;  let make_test p label s =  Printf.printf "point (%.5f,%.5f) is " p.x p.y;  print_string (if point_is_inside s p                then "INSIDE "                else "OUTSIDE ");  print_endline label;;;  let () =  let test_points = [    (p 5. 5.); (p 5. 8.); (p 2. 2.); (p 0. 0.);    (p 10. 10.); (p 2.5 2.5); (p 0.01 5.);    (p 2.2 7.4); (p 0. 5.); (p 10. 5.); (p (-4.) 10.) ] in   List.iter (fun p ->    make_test p "square"     square;    make_test p "squarehole" squarehole;    make_test p "strange"    strange;    make_test p "esa"        esa;    print_newline()  ) test_points;;;

Perl

use strict;use List::Util qw(max min); sub point_in_polygon{    my ( $point,$polygon ) = @_;     my $count = 0; foreach my$side ( @$polygon ) {$count++ if ray_intersect_segment($point,$side);    }    return ($count % 2 == 0) ? 0 : 1;} my$eps = 0.0001;my $inf = 1e600; sub ray_intersect_segment{ my ($point, $segment) = @_; my ($A, $B) = @$segment;     my @P = @$point; # copy it ($A, $B) = ($B, $A) if$A->[1] > $B->[1];$P[1] += $eps if ($P[1] == $A->[1]) || ($P[1] == $B->[1]); return 0 if ($P[1] < $A->[1]) || ($P[1] > $B->[1]) || ($P[0] > max($A->[0],$B->[1]) );    return 1 if $P[0] < min($A->[0], $B->[0]); my$m_red = ($A->[0] !=$B->[0]) ? ( $B->[1] -$A->[1] )/($B->[0] -$A->[0]) : $inf; my$m_blue = ($A->[0] !=$P[0]) ? ( $P[1] -$A->[1] )/($P[0] -$A->[0]) : $inf; return ($m_blue >= $m_red) ? 1 : 0;} Testing: # the following are utilities to use the same Fortran data...sub point{ [shift, shift];}sub create_polygon{ my ($pts, $sides) = @_; my @poly; for(my$i = 0; $i <$#$sides;$i += 2) {	push @poly, [ $pts->[$sides->[$i]-1],$pts->[$sides->[$i+1]-1] ];    }    \@poly;} my @pts = ( point(0,0), point(10,0), point(10,10), point(0,10), 	    point(2.5,2.5), point(7.5,2.5), point(7.5,7.5), point(2.5,7.5), 	    point(0,5), point(10,5), 	    point(3,0), point(7,0), point(7,10), point(3,10) ); my %pgs = (    squared => create_polygon(\@pts, [ 1,2, 2,3, 3,4, 4,1 ] ),    squaredhole => create_polygon(\@pts, [ 1,2, 2,3, 3,4, 4,1, 5,6, 6,7, 7,8, 8,5 ] ),    strange => create_polygon(\@pts, [ 1,5, 5,4, 4,8, 8,7, 7,3, 3,2, 2,5 ] ),    exagon => create_polygon(\@pts, [ 11,12, 12,10, 10,13, 13,14, 14,9, 9,11 ]) ,); my @p = ( point(5,5), point(5, 8), point(-10, 5), point(0,5), point(10,5), &	  point(8,5), point(10,10) ); foreach my $pol ( sort keys %pgs ) { no strict 'refs'; print "$pol\n";    my $rp =$pgs{$pol}; foreach my$tp ( @p ) {	print "\t($tp->[0],$tp->[1]) " .            ( point_in_polygon($tp,$rp) ? "INSIDE" : "OUTSIDE" ) . "\n";    }}

Ursala

This function takes a point ${\displaystyle (x,y)}$ and a polygon ${\displaystyle \langle (x_{1},y_{1})\dots (x_{n},y_{n})\rangle }$ to a true value if the point is enclosed by the polygon and a false value if it's outside, using the algorithm described above. For points on the boundary the result is unspecified.

#import flo in = @lrzyCipPX ~|afatPRZaq ~&EZ+fleq~~lrPrbr2G&& ~&B+fleq~~lrPrbl2G!| -&   ~&Y+ ~~lrPrbl2G fleq,   ^E([email protected],@rl fleq\0.)^/~&lr ^(~&r,times)^/[email protected] vid+ minus~~rbbI&-

This test program tries it on a couple of examples.

#cast %bL examples =  in* <   ((0.5,0.6),<(0.,0.),(1.,2.),(1.,0.)>),   ((0.5,0.6),<(0.,0.),(1.,1.),(1.,0.)>)>

output:

<true,false>

Visual Basic .NET

Translation of: Java
Imports System.Math Module RayCasting     Private square As Integer()() = {New Integer() {0, 0}, New Integer() {20, 0}, New Integer() {20, 20}, New Integer() {0, 20}}    Private squareHole As Integer()() = {New Integer() {0, 0}, New Integer() {20, 0}, New Integer() {20, 20}, New Integer() {0, 20}, New Integer() {5, 5}, New Integer() {15, 5}, New Integer() {15, 15}, New Integer() {5, 15}}    Private strange As Integer()() = {New Integer() {0, 0}, New Integer() {5, 5}, New Integer() {0, 20}, New Integer() {5, 15}, New Integer() {15, 15}, New Integer() {20, 20}, New Integer() {20, 0}}    Private hexagon As Integer()() = {New Integer() {6, 0}, New Integer() {14, 0}, New Integer() {20, 10}, New Integer() {14, 20}, New Integer() {6, 20}, New Integer() {0, 10}}    Private shapes As Integer()()() = {square, squareHole, strange, hexagon}     Public Sub Main()        Dim testPoints As Double()() = {New Double() {10, 10}, New Double() {10, 16}, New Double() {-20, 10}, New Double() {0, 10}, New Double() {20, 10}, New Double() {16, 10}, New Double() {20, 20}}         For Each shape As Integer()() In shapes            For Each point As Double() In testPoints                Console.Write(String.Format("{0} ", Contains(shape, point).ToString.PadLeft(7)))            Next            Console.WriteLine()        Next    End Sub     Private Function Contains(shape As Integer()(), point As Double()) As Boolean         Dim inside As Boolean = False        Dim length As Integer = shape.Length         For i As Integer = 0 To length - 1            If Intersects(shape(i), shape((i + 1) Mod length), point) Then                inside = Not inside            End If        Next         Return inside    End Function     Private Function Intersects(a As Integer(), b As Integer(), p As Double()) As Boolean         If a(1) > b(1) Then Return Intersects(b, a, p)        If p(1) = a(1) Or p(1) = b(1) Then p(1) += 0.0001        If p(1) > b(1) Or p(1) < a(1) Or p(0) > Max(a(0), b(0)) Then Return False        If p(0) < Min(a(0), b(0)) Then Return True        Dim red As Double = (p(1) - a(1)) / (p(0) - a(0))        Dim blue As Double = (b(1) - a(1)) / (b(0) - a(0))         Return red >= blue    End FunctionEnd Module
Output:
   True    True   False   False    True    True   False
False    True   False   False    True    True   False
True   False   False   False    True    True   False
True    True   False   False    True    True   False


zkl

Translation of: Perl 6
const E = 0.0001; fcn rayHitsSeg([(Px,Py)],[(Ax,Ay)],[(Bx,By)]){ // --> Bool   if(Py==Ay or Py==By) Py+=E;   if(Py<Ay or Py>By or Px>Ax.max(Bx)) False   else if(Px<Ax.min(Bx))              True   else try{ ( (Py - Ay)/(Px - Ax) )>=( (By - Ay)/(Bx - Ax) ) } //blue>=red        catch(MathError){ Px==Ax } // divide by zero == ∞, only blue?} fcn pointInPoly(point, polygon){ // --> Bool, polygon is ( (a,b),(c,d).. )   polygon.filter('wrap([(ab,cd)]){ a,b:=ab; c,d:=cd;      if(b<=d) rayHitsSeg(point,ab,cd); // left point has smallest y coordinate      else     rayHitsSeg(point,cd,ab);   })   .len().isOdd; // True if crossed an odd number of borders ie inside polygon}
polys:=T( //(name,( ((a,b),(c,d)),((a,b),(c,d))... ), ... )==(nm,(ln,ln..) ..)    T("squared",	T(T(T( 0.0,  0.0), T(10.0,  0.0)),	  T(T(10.0,  0.0), T(10.0, 10.0)),	  T(T(10.0, 10.0), T( 0.0, 10.0)),	  T(T( 0.0, 10.0), T( 0.0,  0.0)))),    T("squaredhole",	T(T(T( 0.0,  0.0), T(10.0,  0.0)),	  T(T(10.0,  0.0), T(10.0, 10.0)),	  T(T(10.0, 10.0), T( 0.0, 10.0)),	  T(T( 0.0, 10.0), T( 0.0,  0.0)),	  T(T( 2.5,  2.5), T( 7.5,  2.5)),	  T(T( 7.5,  2.5), T( 7.5,  7.5)),	  T(T( 7.5,  7.5), T( 2.5,  7.5)),	  T(T( 2.5,  7.5), T( 2.5,  2.5)))),    T("strange",	T(T(T( 0.0,  0.0), T( 2.5,  2.5)),	  T(T( 2.5,  2.5), T( 0.0, 10.0)),	  T(T( 0.0, 10.0), T( 2.5,  7.5)),	  T(T( 2.5,  7.5), T( 7.5,  7.5)),	  T(T( 7.5,  7.5), T(10.0, 10.0)),	  T(T(10.0, 10.0), T(10.0,  0.0)),	  T(T(10.0,  0.0), T( 2.5,  2.5)),	  T(T( 2.5,  2.5), T( 0.0,  0.0)))),  # conjecturally close polygon    T("exagon",	T(T(T( 3.0,  0.0), T( 7.0,  0.0)),	  T(T( 7.0,  0.0), T(10.0,  5.0)),	  T(T(10.0,  5.0), T( 7.0, 10.0)),	  T(T( 7.0, 10.0), T( 3.0, 10.0)),	  T(T( 3.0, 10.0), T( 0.0,  5.0)),	  T(T( 0.0,  5.0), T( 3.0,  0.0)))),); testPoints:=T(	  T(  5.0,  5.0),	  T(  5.0,  8.0),	  T(-10.0,  5.0),	  T(  0.0,  5.0),	  T( 10.0,  5.0),	  T(  8.0,  5.0),	  T( 10.0, 10.0)); foreach name,polywanna in (polys){   name.println();   foreach testPoint in (testPoints){      println("\t(%6.1f,%6.1f)\t".fmt(testPoint.xplode()),         pointInPoly(testPoint,polywanna) and "IN" or "OUT");   }}
Output:
squared
(   5.0,   5.0)	IN
(   5.0,   8.0)	IN
( -10.0,   5.0)	OUT
(   0.0,   5.0)	OUT
(  10.0,   5.0)	IN
(   8.0,   5.0)	IN
(  10.0,  10.0)	OUT
squaredhole
(   5.0,   5.0)	OUT
(   5.0,   8.0)	IN
( -10.0,   5.0)	OUT
(   0.0,   5.0)	OUT
(  10.0,   5.0)	IN
(   8.0,   5.0)	IN
(  10.0,  10.0)	OUT
strange
(   5.0,   5.0)	IN
(   5.0,   8.0)	OUT
( -10.0,   5.0)	OUT
(   0.0,   5.0)	OUT
(  10.0,   5.0)	IN
(   8.0,   5.0)	IN
(  10.0,  10.0)	OUT
exagon
(   5.0,   5.0)	IN
(   5.0,   8.0)	IN
( -10.0,   5.0)	OUT
(   0.0,   5.0)	OUT
(  10.0,   5.0)	IN
(   8.0,   5.0)	IN
(  10.0,  10.0)	OUT