1 ! Copyright (c) 2012 Anonymous
2 ! See http://factorcode.org/license.txt for BSD license.
3 USING: kernel prettyprint sequences arrays math math.vectors ;
4 IN: rosetta-code.raycasting
7 ! http://rosettacode.org/wiki/Ray-casting_algorithm
9 ! Given a point and a polygon, check if the point is inside or
10 ! outside the polygon using the ray-casting algorithm.
12 ! A pseudocode can be simply:
15 ! foreach side in polygon:
16 ! if ray_intersects_segment(P,side) then
18 ! if is_odd(count) then
23 ! Where the function ray_intersects_segment return true if the
24 ! horizontal ray starting from the point P intersects the side
25 ! (segment), false otherwise.
27 ! An intuitive explanation of why it works is that every time we
28 ! cross a border, we change "country" (inside-outside, or
29 ! outside-inside), but the last "country" we land on is surely
30 ! outside (since the inside of the polygon is finite, while the
31 ! ray continues towards infinity). So, if we crossed an odd number
32 ! of borders we was surely inside, otherwise we was outside; we
33 ! can follow the ray backward to see it better: starting from
34 ! outside, only an odd number of crossing can give an inside:
35 ! outside-inside, outside-inside-outside-inside, and so on (the -
36 ! represents the crossing of a border).
38 ! So the main part of the algorithm is how we determine if a ray
39 ! intersects a segment. The following text explain one of the
42 ! Looking at the image on the right, we can easily be convinced
43 ! of the fact that rays starting from points in the hatched area
44 ! (like P1 and P2) surely do not intersect the segment AB. We also
45 ! can easily see that rays starting from points in the greenish
46 ! area surely intersect the segment AB (like point P3).
48 ! So the problematic points are those inside the white area (the
49 ! box delimited by the points A and B), like P4.
51 ! Let us take into account a segment AB (the point A having y
52 ! coordinate always smaller than B's y coordinate, i.e. point A is
53 ! always below point B) and a point P. Let us use the cumbersome
54 ! notation PAX to denote the angle between segment AP and AX,
55 ! where X is always a point on the horizontal line passing by A
56 ! with x coordinate bigger than the maximum between the x
57 ! coordinate of A and the x coordinate of B. As explained
58 ! graphically by the figures on the right, if PAX is greater than
59 ! the angle BAX, then the ray starting from P intersects the
60 ! segment AB. (In the images, the ray starting from PA does not
61 ! intersect the segment, while the ray starting from PB in the
62 ! second picture, intersects the segment).
64 ! Points on the boundary or "on" a vertex are someway special
65 ! and through this approach we do not obtain coherent results.
66 ! They could be treated apart, but it is not necessary to do so.
68 ! An algorithm for the previous speech could be (if P is a
69 ! point, Px is its x coordinate):
71 ! ray_intersects_segment:
72 ! P : the point from which the ray starts
73 ! A : the end-point of the segment with the smallest y coordinate
74 ! (A must be "below" B)
75 ! B : the end-point of the segment with the greatest y coordinate
76 ! (B must be "above" A)
77 ! if Py = Ay or Py = By then
80 ! if Py < Ay or Py > By then
82 ! else if Px > max(Ax, Bx) then
85 ! if Px < min(Ax, Bx) then
89 ! m_red ← (By - Ay)/(Bx - Ax)
94 ! m_blue ← (Py - Ay)/(Px - Ax)
98 ! if m_blue ≥ m_red then
106 ! (To avoid the "ray on vertex" problem, the point is moved
107 ! upward of a small quantity ε)
109 : between ( a b x -- ? ) [ last ] tri@ [ < ] curry bi@ xor ;
111 : lincomb ( a b x -- w )
115 neg 2dup + [ / ] curry bi@
116 [ [ v*n ] curry ] bi@ bi* v+ ;
118 : leftof ( a b x -- ? ) dup [ lincomb ] dip [ first ] bi@ > ;
120 : ray ( a b x -- ? ) [ between ] [ leftof ] 3bi and ;
122 : raycast ( poly x -- ? )
123 [ dup first suffix [ rest-slice ] [ but-last-slice ] bi ] dip