Ray-casting algorithm: Difference between revisions

Content added Content deleted

Inline

Revision as of 16:58, 24 May 2009

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 P₁ and P₂) 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 P₃).

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

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 P_A does not intersect the segment, while the ray starting from P_B 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 ε)

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.

Fortran

Works with: Fortran version 95 and later

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>

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>