Ray-casting algorithm: Difference between revisions
(+ Tcl) |
m (→{{header|C}}: removed old-work comments; alignment of \ into define) |
||
Line 129: | Line 129: | ||
if ( s->x <= MIN(a->x, b->x) ) return true; |
if ( s->x <= MIN(a->x, b->x) ) return true; |
||
double ca = (a->x != b->x) ? coeff_ang(a,b) : |
double ca = (a->x != b->x) ? coeff_ang(a,b) : |
||
HUGE_VAL |
HUGE_VAL; |
||
point_t *my = minP(a,b,y); |
point_t *my = minP(a,b,y); |
||
double cp = (s->x - my->x) ? coeff_ang(my, s) : |
double cp = (s->x - my->x) ? coeff_ang(my, s) : |
||
HUGE_VAL |
HUGE_VAL; |
||
if ( cp >= ca ) return true; |
if ( cp >= ca ) return true; |
||
return false; |
return false; |
||
Line 153: | Line 153: | ||
Testing: |
Testing: |
||
<lang c>#define MAKE_TEST(S) do { \ |
<lang c>#define MAKE_TEST(S) do { \ |
||
printf("point (%.5f,%.5f) is ", test_points[i].x, test_points[i].y); \ |
printf("point (%.5f,%.5f) is ", test_points[i].x, test_points[i].y); \ |
||
if ( point_is_inside(&S, &test_points[i]) ) \ |
if ( point_is_inside(&S, &test_points[i]) ) \ |
||
printf("INSIDE " #S "\n"); \ |
printf("INSIDE " #S "\n"); \ |
||
else \ |
else \ |
||
printf("OUTSIDE " #S "\n"); \ |
printf("OUTSIDE " #S "\n"); \ |
||
} while(0); |
} while(0); |
||
Revision as of 09:15, 24 May 2009
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 Ray-casting algorithm. 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.
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 ≤ B 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
(As it can be seen, if the ray passes through a vertex, it is not considered intersection; there exists a set of "situations" where this behaviour brings to wrong results; to avoid this, we should "move" the point up or down by 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}
};
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) )
inline bool hseg_intersect_seg(point_t *s, point_t *a, point_t *b) {
if ( (s->y >= MAX(a->y, b->y) || s->y <= 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); 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.
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 {
foreach {x0 y0} $polygon break 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} {
foreach {x y} $point break foreach {xa xb ya yb} $line break set xout [expr {max($xa,$xb) + 1}] lines_cross [list $x $y $xout $y] $line
} proc lines_cross {line1 line2} {
foreach {xa ya xb yb} $line1 break foreach {xc yc xd yd} $line2 break set det [expr {($xb - $xa) * ($yd - $yc) - ($xd - $xc) * ($yb - $ya)}] if {$det == 0} {return 0} ;#-- parallel or collinear set nump [expr {($xd - $xc) * ($yb + $ya) - ($xb + $xa) * ($yd - $yc) + 2. * ($xc * $yd - $xd * $yc)}] if {abs($nump) > abs($det)} {return 0} ;#-- cross outside line1 set numq [expr {($xd + $xc) * ($yb - $ya) - ($xb - $xa) * ($yd + $yc) + 2. * ($xb * $ya - $xa * $yb)}] if {abs($numq) > abs($det)} {return 0} ;#-- cross outside line2 return 1
}
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}
} {
puts "$point in $poly = [point_in_polygon $point $poly]"
} </lang>