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 was surely inside, otherwise we was 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 ε)
AutoHotkey
<lang ahk> PiP(P,N) ; P is point, N is number of sides {
count := 0 Loop %N% If ray_intersects_segment(P,A_Index) count++ if mod(count,2) return false ; P isn't in the polygon else return true ; P is in the polygon
} </lang>
C
Required includes and definitions:
<lang c>#include <stdio.h>
- include <stdlib.h>
- include <stdbool.h>
- include <math.h>
typedef struct {
double x, y;
} point_t;
typedef struct {
point_t *vertices; int edges[];
} polygon_t;</lang>
Polygons for testing:
<lang c>point_t square_v[] = {
{0,0}, {10,0}, {10,10}, {0,10},
{2.5,2.5}, {7.5,0.1}, {7.5,7.5}, {2.5,7.5}
};
point_t esa_v[] = {
{3,0}, {7,0}, {10,5}, {7,10}, {3,10}, {0,5}
};
polygon_t esa = {
esa_v,
{ 0,1, 1,2, 2,3, 3,4, 4,5, 5,0, -1 }
};
polygon_t square = {
square_v,
{ 0,1, 1,2, 2,3, 3,0, -1 }
};
polygon_t squarehole = {
square_v,
{ 0,1, 1,2, 2,3, 3,0,
4,5, 5,6, 6,7, 7,4, -1 }
};
polygon_t strange = {
square_v,
{ 0,4, 4,3, 3,7, 7,6, 6,2, 2,1, 1,5, 5,0, -1 }
};</lang>
Check for intersection code:
<lang c>#define MAX(A,B) (((A)>(B))?(A):(B))
- define MIN(A,B) (((A)>(B))?(B):(A))
- define minP(A,B,C) ( (((A)->C) > ((B)->C)) ? (B) : (A) )
- define coeff_ang(PA,PB) ( ((PB)->y - (PA)->y) / ((PB)->x - (PA)->x) )
- define EPS 0.00001
inline bool hseg_intersect_seg(point_t *s, point_t *a, point_t *b) {
double eps = 0.0;
if ( s->y == MAX(a->y, b->y) ||
s->y == MIN(a->y, b->y) ) eps = EPS;
if ( (s->y + eps) > MAX(a->y, b->y) ||
(s->y + eps) < MIN(a->y, b->y) ||
s->x > MAX(a->x, b->x) ) return false; if ( s->x <= MIN(a->x, b->x) ) return true; double ca = (a->x != b->x) ? coeff_ang(a,b) : HUGE_VAL; point_t *my = minP(a,b,y); double cp = (s->x - my->x) ? coeff_ang(my, s) : HUGE_VAL; if ( cp >= ca ) return true; return false;
}</lang>
The ray-casting algorithm:
<lang c>bool point_is_inside(polygon_t *poly, point_t *pt) {
int cross = 0, i;
for(i=0; poly->edges[i] != -1 ; i+=2) {
if ( hseg_intersect_seg(pt,
&poly->vertices[ poly->edges[i] ],
&poly->vertices[ poly->edges[i+1] ]) )
cross++;
}
return !(cross%2 == 0);
}</lang>
Testing:
<lang c>#define MAKE_TEST(S) do { \
printf("point (%.5f,%.5f) is ", test_points[i].x, test_points[i].y); \
if ( point_is_inside(&S, &test_points[i]) ) \
printf("INSIDE " #S "\n"); \
else \
printf("OUTSIDE " #S "\n"); \
} while(0);
int main() {
point_t test_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}};
int i;
for(i=0; i < sizeof(test_points)/sizeof(point_t); i++) {
MAKE_TEST(square);
MAKE_TEST(squarehole);
MAKE_TEST(strange);
MAKE_TEST(esa);
printf("\n");
}
return EXIT_SUCCESS;
}</lang>
The test's output reveals the meaning of coherent results: a point on the leftmost vertical side of the square (coordinate 0,5) is considered outside; while a point on the rightmost vertical side of the square (coordinate 10,5) is considered inside, but on the top-right vertex (coordinate 10,10), it is considered outside again.
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.
<lang lisp>(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)))))))))</lang>
Testing code
<lang lisp>(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)))))</lang>
Fortran
The following code uses the Points_Module defined here.
This module defines "polygons".
<lang fortran>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</lang>
The ray casting algorithm module:
<lang fortran>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</lang>
Testing
<lang fortran>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", "exagon" /)
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</lang>
Haskell
<lang haskell>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) l@ is true if the ray {(x, py) | x ≥ px} intersects l. -} intersects (px, _) (Vert xint) = px <= xint intersects (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 == xint onLine (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 -> Bool between x a b | a > b = b <= x && x <= a
| otherwise = a <= x && x <= b
inPolygon :: Point -> Polygon -> Bool inPolygon 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 &&
(not (py == ay) || by < py) &&
(not (py == by) || ay < py)
((ax, ay), (bx, by)) = side
line = carrier side</lang>
J
<lang j>NB.*crossPnP v point in closed polygon, crossing number NB. bool=. points crossPnP polygon crossPnP=: 4 : 0"2
'X Y'=. |:x
'x0 y0 x1 y1'=. |:2 ,/\^:(2={:@$@]) y
p1=. ((y0<:/Y)*. y1>/Y) +. (y0>/Y)*. y1<:/Y
p2=. (x0-/X) < (x0-x1) * (y0-/Y) % (y0 - y1)
2|+/ p1*.p2
)</lang>
Sample data: <lang j>SQUAREV=: 0 0 , 10 0 , 10 10 ,: 0 10 SQUAREV=: 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) , .{ ESAV SQUARE=: (0 1,1 2,2 3,:3 0) , .{ SQUAREV SQUAREHOLE=: (0 1,1 2,2 3,3 0,4 5,5 6,6 7,:7 4) , .{ SQUAREV STRANGE=: (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</lang>
Testing: <lang j> (<POINTS) crossPnP every ESA;SQUARE;SQUAREHOLE;STRANGE 1 1 1 0 0 1 1 1 0 1 0 1 1 1 0 0 1 1 1 0 1 0 0 1 1 0 0 1 1 1 0 1 0 1 0 0 0 0 0 0 1 0 1 0</lang>
OCaml
<lang ocaml>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 begin
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
if cp >= ca then (true) else (false)
end
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;
(not ( (!cross mod 2) = 0))
let make_test p label s =
Printf.printf "point (%.5f,%.5f) is " p.x p.y; if point_is_inside s p then print_string "INSIDE " else print_string "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;
- </lang>
Perl
<lang 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 = $segment->[0]; my $B = $segment->[1];
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;
}</lang>
Testing:
<lang perl># the following are utilities to use the same Fortran data... sub point {
return [shift, shift];
} sub create_polygon {
my ($pts, $sides) = @_;
my @poly = ();
for(my $i = 0; $i < (scalar(@$sides)-1); $i += 2) {
push @poly, [ $pts->[$sides->[$i]-1], $pts->[$sides->[$i+1]-1] ];
} return @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 @squared = create_polygon(\@pts, [ 1,2, 2,3, 3,4, 4,1 ] ); my @squaredhole = create_polygon(\@pts, [ 1,2, 2,3, 3,4, 4,1, 5,6, 6,7, 7,8, 8,5 ] ); my @strange = create_polygon(\@pts, [ 1,5, 5,4, 4,8, 8,7, 7,3, 3,2, 2,5 ] ); my @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 ( "squared", "squaredhole", "strange", "exagon" ) {
print "$pol\n";
my @rp = eval('@' . $pol);
foreach my $tp ( @p ) {
print "\t($tp->[0],$tp->[1]) " .
( point_in_polygon($tp, \@rp) ? "INSIDE" : "OUTSIDE" ) . "\n"; }
}</lang>
Python
from pprint import pprint as pp import sys
Pt = namedtuple('Pt', 'x, y') # Point Edge = namedtuple('Edge', 'a, b') # Polygon edge from a to b Poly = namedtuple('Poly', 'name, edges') # Polygon
_eps = 0.00001 _huge = sys.float_info.max _tiny = sys.float_info.min
def rayintersectseg(p, edge):
takes a point p=Pt() and an edge of two endpoints a,b=Pt() of a line segment returns boolean
a,b = edge
if a.y > b.y:
a,b = b,a
if p.y == a.y or p.y == b.y:
p = Pt(p.x, p.y + _eps)
intersect = False
if (p.y > b.y or p.y < a.y) or (
p.x > max(a.x, b.x)):
return False
if p.x < min(a.x, b.x):
intersect = True
else:
if abs(a.x - b.x) > _tiny:
m_red = (b.y - a.y) / float(b.x - a.x)
else:
m_red = _huge
if abs(a.x - p.x) > _tiny:
m_blue = (p.y - a.y) / float(p.x - a.x)
else:
m_blue = _huge
intersect = m_blue >= m_red
return intersect
def _odd(x): return x%2 == 1
def ispointinside(p, poly):
ln = len(poly)
return _odd(sum(rayintersectseg(p, edge)
for edge in poly.edges ))
def polypp(poly):
print "\n Polygon(name='%s', edges=(" % poly.name
print ' ', ',\n '.join(str(e) for e in poly.edges) + '\n ))'
if __name__ == '__main__':
polys = [
Poly(name='square', edges=(
Edge(a=Pt(x=0, y=0), b=Pt(x=10, y=0)),
Edge(a=Pt(x=10, y=0), b=Pt(x=10, y=10)),
Edge(a=Pt(x=10, y=10), b=Pt(x=0, y=10)),
Edge(a=Pt(x=0, y=10), b=Pt(x=0, y=0))
)),
Poly(name='square_hole', edges=(
Edge(a=Pt(x=0, y=0), b=Pt(x=10, y=0)),
Edge(a=Pt(x=10, y=0), b=Pt(x=10, y=10)),
Edge(a=Pt(x=10, y=10), b=Pt(x=0, y=10)),
Edge(a=Pt(x=0, y=10), b=Pt(x=0, y=0)),
Edge(a=Pt(x=2.5, y=2.5), b=Pt(x=7.5, y=2.5)),
Edge(a=Pt(x=7.5, y=2.5), b=Pt(x=7.5, y=7.5)),
Edge(a=Pt(x=7.5, y=7.5), b=Pt(x=2.5, y=7.5)),
Edge(a=Pt(x=2.5, y=7.5), b=Pt(x=2.5, y=2.5))
)),
Poly(name='strange', edges=(
Edge(a=Pt(x=0, y=0), b=Pt(x=2.5, y=2.5)),
Edge(a=Pt(x=2.5, y=2.5), b=Pt(x=0, y=10)),
Edge(a=Pt(x=0, y=10), b=Pt(x=2.5, y=7.5)),
Edge(a=Pt(x=2.5, y=7.5), b=Pt(x=7.5, y=7.5)),
Edge(a=Pt(x=7.5, y=7.5), b=Pt(x=10, y=10)),
Edge(a=Pt(x=10, y=10), b=Pt(x=10, y=0)),
Edge(a=Pt(x=10, y=0), b=Pt(x=2.5, y=2.5))
)),
Poly(name='exagon', edges=(
Edge(a=Pt(x=3, y=0), b=Pt(x=7, y=0)),
Edge(a=Pt(x=7, y=0), b=Pt(x=10, y=5)),
Edge(a=Pt(x=10, y=5), b=Pt(x=7, y=10)),
Edge(a=Pt(x=7, y=10), b=Pt(x=3, y=10)),
Edge(a=Pt(x=3, y=10), b=Pt(x=0, y=5)),
Edge(a=Pt(x=0, y=5), b=Pt(x=3, y=0))
)),
]
testpoints = (Pt(x=5, y=5), Pt(x=5, y=8),
Pt(x=-10, y=5), Pt(x=0, y=5),
Pt(x=10, y=5), Pt(x=8, y=5),
Pt(x=10, y=10))
print "\n TESTING WHETHER POINTS ARE WITHIN POLYGONS"
for poly in polys:
polypp(poly)
print ' ', '\t'.join("%s: %s" % (p, ispointinside(p, poly))
for p in testpoints[:3])
print ' ', '\t'.join("%s: %s" % (p, ispointinside(p, poly))
for p in testpoints[3:6])
print ' ', '\t'.join("%s: %s" % (p, ispointinside(p, poly))
for p in testpoints[6:])</lang>Sample output
TESTING WHETHER POINTS ARE WITHIN POLYGONS
Polygon(name='square', edges=(
Edge(a=Pt(x=0, y=0), b=Pt(x=10, y=0)),
Edge(a=Pt(x=10, y=0), b=Pt(x=10, y=10)),
Edge(a=Pt(x=10, y=10), b=Pt(x=0, y=10)),
Edge(a=Pt(x=0, y=10), b=Pt(x=0, y=0))
))
Pt(x=5, y=5): True Pt(x=5, y=8): True Pt(x=-10, y=5): False
Pt(x=0, y=5): False Pt(x=10, y=5): True Pt(x=8, y=5): True
Pt(x=10, y=10): False
Polygon(name='square_hole', edges=(
Edge(a=Pt(x=0, y=0), b=Pt(x=10, y=0)),
Edge(a=Pt(x=10, y=0), b=Pt(x=10, y=10)),
Edge(a=Pt(x=10, y=10), b=Pt(x=0, y=10)),
Edge(a=Pt(x=0, y=10), b=Pt(x=0, y=0)),
Edge(a=Pt(x=2.5, y=2.5), b=Pt(x=7.5, y=2.5)),
Edge(a=Pt(x=7.5, y=2.5), b=Pt(x=7.5, y=7.5)),
Edge(a=Pt(x=7.5, y=7.5), b=Pt(x=2.5, y=7.5)),
Edge(a=Pt(x=2.5, y=7.5), b=Pt(x=2.5, y=2.5))
))
Pt(x=5, y=5): False Pt(x=5, y=8): True Pt(x=-10, y=5): False
Pt(x=0, y=5): False Pt(x=10, y=5): True Pt(x=8, y=5): True
Pt(x=10, y=10): False
Polygon(name='strange', edges=(
Edge(a=Pt(x=0, y=0), b=Pt(x=2.5, y=2.5)),
Edge(a=Pt(x=2.5, y=2.5), b=Pt(x=0, y=10)),
Edge(a=Pt(x=0, y=10), b=Pt(x=2.5, y=7.5)),
Edge(a=Pt(x=2.5, y=7.5), b=Pt(x=7.5, y=7.5)),
Edge(a=Pt(x=7.5, y=7.5), b=Pt(x=10, y=10)),
Edge(a=Pt(x=10, y=10), b=Pt(x=10, y=0)),
Edge(a=Pt(x=10, y=0), b=Pt(x=2.5, y=2.5))
))
Pt(x=5, y=5): True Pt(x=5, y=8): False Pt(x=-10, y=5): False
Pt(x=0, y=5): False Pt(x=10, y=5): True Pt(x=8, y=5): True
Pt(x=10, y=10): False
Polygon(name='exagon', edges=(
Edge(a=Pt(x=3, y=0), b=Pt(x=7, y=0)),
Edge(a=Pt(x=7, y=0), b=Pt(x=10, y=5)),
Edge(a=Pt(x=10, y=5), b=Pt(x=7, y=10)),
Edge(a=Pt(x=7, y=10), b=Pt(x=3, y=10)),
Edge(a=Pt(x=3, y=10), b=Pt(x=0, y=5)),
Edge(a=Pt(x=0, y=5), b=Pt(x=3, y=0))
))
Pt(x=5, y=5): True Pt(x=5, y=8): True Pt(x=-10, y=5): False
Pt(x=0, y=5): False Pt(x=10, y=5): True Pt(x=8, y=5): True
Pt(x=10, y=10): False
Helper routine to convert Fortran Polygons and points to Python
<lang python>def _convert_fortran_shapes():
point = Pt
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))
p = (point(5,5), point(5, 8), point(-10, 5), point(0,5), point(10,5),
point(8,5), point(10,10) )
def create_polygon(pts,vertexindex):
return [tuple(Edge(pts[vertexindex[i]-1], pts[vertexindex[i+1]-1])
for i in range(0, len(vertexindex), 2) )]
polys=[]
polys += create_polygon(pts, ( 1,2, 2,3, 3,4, 4,1 ) )
polys += create_polygon(pts, ( 1,2, 2,3, 3,4, 4,1, 5,6, 6,7, 7,8, 8,5 ) )
polys += create_polygon(pts, ( 1,5, 5,4, 4,8, 8,7, 7,3, 3,2, 2,5 ) )
polys += create_polygon(pts, ( 11,12, 12,10, 10,13, 13,14, 14,9, 9,11 ) )
names = ( "square", "square_hole", "strange", "exagon" )
polys = [Poly(name, edges)
for name, edges in zip(names, polys)]
print 'polys = ['
for p in polys:
print " Poly(name='%s', edges=(" % p.name
print ' ', ',\n '.join(str(e) for e in p.edges) + '\n )),'
print ' ]'
_convert_fortran_shapes()</lang>R
<lang R>point_in_polygon <- function(polygon, p) {
count <- 0
for(side in polygon) {
if ( ray_intersect_segment(p, side) ) {
count <- count + 1
}
}
if ( count %% 2 == 1 )
"INSIDE"
else
"OUTSIDE"
}
ray_intersect_segment <- function(p, side) {
eps <- 0.0001
a <- side$A
b <- side$B
if ( a$y > b$y ) {
a <- side$B
b <- side$A
}
if ( (p$y == a$y) || (p$y == b$y) ) {
p$y <- p$y + eps
}
if ( (p$y < a$y) || (p$y > b$y) )
return(FALSE)
else if ( p$x > max(a$x, b$x) )
return(FALSE)
else {
if ( p$x < min(a$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
return( m_blue >= m_red )
}
}
}</lang>
<lang R>######## utility functions #########
point <- function(x,y) list(x=x, y=y)
- pts = list(p1, p2, ... )... coords
- segs = list(c(1,2), c(2,1) ...) indices
createPolygon <- function(pts, segs) {
pol <- list()
for(pseg in segs) {
pol <- c(pol, list(list(A=pts[[pseg[1]]], B=pts[[pseg[2]]])))
}
pol
}</lang>
<lang R>#### testing ####
pts <- list(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))
polygons <-
list(
square = createPolygon(pts, list(c(1,2), c(2,3), c(3,4), c(4,1))),
squarehole = createPolygon(pts, list(c(1,2), c(2,3), c(3,4), c(4,1), c(5,6), c(6,7), c(7,8), c(8,5))),
exagon = createPolygon(pts, list(c(11,12), c(12,10), c(10,13), c(13,14), c(14,9), c(9,11)))
)
testpoints <-
list(
point(5,5), point(5, 8), point(-10, 5), point(0,5), point(10,5),
point(8,5), point(9.9,9.9)
)
for(p in testpoints) {
for(polysi in 1:length(polygons)) {
cat(sprintf("point (%lf, %lf) is %s polygon (%s)\n",
p$x, p$y, point_in_polygon(polygonspolysi, p), names(polygons[polysi])))
}
}</lang>
Smalltalk
The class Segment holds the code to test if a ray starting from a point (and going towards "right") intersects the segment (method doesIntersectRayFrom); while the class Polygon hosts the code to test if a point is inside the polygon (method pointInside).
<lang smalltalk>Object subclass: Segment [
|pts|
Segment class >> new: points [ |a|
a := super new.
^ a init: points
]
init: points [ pts := points copy. ^self ]
endPoints [ ^pts ]
"utility methods"
first [ ^ pts at: 1]
second [ ^ pts at: 2]
leftmostEndPoint [
^ (self first x > self second x) ifTrue: [ self second ] ifFalse: [ self first ]
]
rightmostEndPoint [
^ (self first x > self second x) ifTrue: [ self first ] ifFalse: [ self second ]
]
topmostEndPoint [
^ (self first y > self second y) ifTrue: [ self first ] ifFalse: [ self second ]
]
bottommostEndPoint [
^ (self first y > self second y) ifTrue: [ self second ] ifFalse: [ self first ]
]
slope [
(pts at: 1) x ~= (pts at: 2) x
ifTrue: [ ^ ((pts at: 1) y - (pts at: 2) y) / ((pts at: 1) x - (pts at: 2) x) ]
ifFalse: [ ^ FloatD infinity ]
]
doesIntersectRayFrom: point [ |p A B|
(point y = (pts at: 1) y) | (point y = (pts at: 2) y)
ifTrue: [ p := Point x: (point x) y: (point y) + 0.00001 ]
ifFalse: [ p := point copy ].
A := self bottommostEndPoint.
B := self topmostEndPoint.
(p y < A y) | (p y > B y) | (p x > (self rightmostEndPoint x))
ifTrue: [ ^false ]
ifFalse: [ (p x < (self leftmostEndPoint x))
ifTrue: [ ^true ]
ifFalse: [ |s|
s := Segment new: { A . point }.
(s slope) >= (self slope) ifTrue: [ ^ true ]
]
].
^false
]
].
Object subclass: Polygon [
|polysegs|
Polygon class >> new [ |a| a := super new. ^ a init. ]
Polygon class >> fromSegments: segments [ |a|
a := super new.
^ a initWithSegments: segments
]
Polygon class >> fromPoints: pts and: indexes [ |a|
a := self new.
indexes do: [ :i |
a addSegment: ( Segment new: { pts at: (i at: 1) . pts at: (i at: 2) } )
].
^ a
]
initWithSegments: segments [
polysegs := segments copy. ^self
]
init [ polysegs := OrderedCollection new. ^ self ]
addSegment: segment [ polysegs add: segment ]
pointInside: point [ |cnt|
cnt := 0.
polysegs do: [ :s | (s doesIntersectRayFrom: point)
ifTrue: [ cnt := cnt + 1 ] ].
^ ( cnt \\ 2 = 0 ) not
]
].</lang>
Testing
<lang smalltalk>|points names polys|
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
}.
names := { 'square' . 'square hole' . 'strange' . 'exagon' }.
polys := OrderedCollection new.
polys add:
(
Polygon fromPoints: points
and: { {1 . 2}. {2 . 3}. {3 . 4}. {4 . 1} }
) ;
add:
(
Polygon fromPoints: points
and: { {1 . 2}. {2 . 3}. {3 . 4}. {4 . 1}. {5 . 6}. {6 . 7}. {7 . 8}. {8 . 5} }
) ;
add:
(
Polygon fromPoints: points
and: { {1 . 5}. {5 . 4}. {4 . 8}. {8 . 7}. {7 . 3}. {3 . 2}. {2 . 5} }
) ;
add:
(
Polygon fromPoints: points
and: { {11 . 12}. {12 . 10}. {10 . 13}. {13 . 14}. {14 . 9}. {9 . 11} }
).
{ 5@5 . 5@8 . -10@5 . 0@5 . 10@5 . 8@5 . 10@10 } do: [ :p |
1 to: 4 do: [ :i |
('point %1 inside %2? %3' %
{ p . names at: i. (polys at: i) pointInside: p }) displayNl
].
' ' displayNl.
]</lang>
Tcl
<lang Tcl>package require Tcl 8.5
proc point_in_polygon {point polygon} {
set count 0
foreach side [sides $polygon] {
if {[ray_intersects_line $point $side]} {
incr count
}
}
expr {$count % 2} ;#-- 1 = odd = true, 0 = even = false
} proc sides polygon {
lassign $polygon x0 y0
foreach {x y} [lrange [lappend polygon $x0 $y0] 2 end] {
lappend res [list $x0 $y0 $x $y]
set x0 $x
set y0 $y
}
return $res
} proc ray_intersects_line {point line} {
lassign $point Px Py
lassign $line Ax Ay Bx By
# Reverse line direction if necessary
if {$By < $Ay} {
lassign $line Bx By Ax Ay
}
# Add epsilon to
if {$Py == $Ay || $Py == $By} {
set Py [expr {$Py + abs($Py)/1e6}]
}
# Bounding box checks
if {$Py < $Ay || $Py > $By || $Px > max($Ax,$Bx)} {
return 0
} elseif {$Px < min($Ax,$Bx)} {
return 1
}
# Compare dot products to compare (cosines of) angles
set mRed [expr {$Ax != $Bx ? ($By-$Ay)/($Bx-$Ax) : Inf}]
set mBlu [expr {$Ax != $Px ? ($Py-$Ay)/($Px-$Ax) : Inf}]
return [expr {$mBlu >= $mRed}]
}
foreach {point poly} {
{0 0} {-1 -1 -1 1 1 1 1 -1}
{2 2} {-1 -1 -1 1 1 1 1 -1}
{0 0} {-2 -2 -2 2 2 2 2 -2 2 -1 1 1 -1 1 -1 -1 1 -1 2 -1}
{1.5 1.5} {-2 -2 -2 2 2 2 2 -2 2 -1 1 1 -1 1 -1 -1 1 -1 2 -1}
{5 5} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{5 8} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{2 2} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{0 0} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{10 10} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{2.5 2.5} {0 0 2.5 2.5 0 10 2.5 7.5 7.5 7.5 10 10 10 0 7.5 0.1}
{-5 5} {3 0 7 0 10 5 7 10 3 10 0 5}
} {
puts "$point in $poly = [point_in_polygon $point $poly]"
}</lang>
Ursala
This function takes a point (x,y) and a polygon <(x1,y1)...(xn,yn)> 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. <lang Ursala>#import flo
in =
@lrzyCipPX ~|afatPRZaq ~&EZ+fleq~~lrPrbr2G&& ~&B+fleq~~lrPrbl2G!| -&
~&Y+ ~~lrPrbl2G fleq, ^E(fleq@lrrPX,@rl fleq\0.)^/~&lr ^(~&r,times)^/minus@llPrll2X vid+ minus~~rbbI&-</lang>
This test program tries it on a couple of examples. <lang Ursala>#cast %bL
examples =
in* <
((0.5,0.6),<(0.,0.),(1.,2.),(1.,0.)>), ((0.5,0.6),<(0.,0.),(1.,1.),(1.,0.)>)></lang>
output:
<true,false>