extra/rosetta-code/raycasting/raycasting.factor

   1 ! Copyright (c) 2012 Anonymous
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: kernel math math.vectors sequences ;
   4 IN: rosetta-code.raycasting
   5
   6
   7 ! http://rosettacode.org/wiki/Ray-casting_algorithm
   8
   9 ! Given a point and a polygon, check if the point is inside or
  10 ! outside the polygon using the ray-casting algorithm.
  11
  12 ! A pseudocode can be simply:
  13
  14 ! count ← 0
  15 ! foreach side in polygon:
  16 !   if ray_intersects_segment(P,side) then
  17 !     count ← count + 1
  18 ! if is_odd(count) then
  19 !   return inside
  20 ! else
  21 !   return outside
  22
  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.
  26
  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).
  37
  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
  40 ! possible ways.
  41
  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).
  47
  48 ! So the problematic points are those inside the white area (the
  49 ! box delimited by the points A and B), like P4.
  50
  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).
  63
  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.
  67
  68 ! An algorithm for the previous speech could be (if P is a
  69 ! point, Px is its x coordinate):
  70
  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
  78 !   Py ← Py + ε
  79 ! end if
  80 ! if Py < Ay or Py > By then
  81 !   return false
  82 ! else if Px > max(Ax, Bx) then
  83 !   return false
  84 ! else
  85 !   if Px < min(Ax, Bx) then
  86 !     return true
  87 !   else
  88 !     if Ax ≠ Bx then
  89 !       m_red ← (By - Ay)/(Bx - Ax)
  90 !     else
  91 !       m_red ← ∞
  92 !     end if
  93 !     if Ax ≠ Px then
  94 !       m_blue ← (Py - Ay)/(Px - Ax)
  95 !     else
  96 !       m_blue ← ∞
  97 !     end if
  98 !     if m_blue ≥ m_red then
  99 !       return true
 100 !     else
 101 !       return false
 102 !     end if
 103 !   end if
 104 ! end if
 105
 106 ! (To avoid the "ray on vertex" problem, the point is moved
 107 ! upward of a small quantity ε)
 108
 109 : between ( a b x -- ? ) [ last ] tri@ [ < ] curry bi@ xor ;
 110
 111 : lincomb ( a b x -- w )
 112     3dup [ last ] tri@
 113     [ - ] curry bi@
 114     nipd
 115     neg 2dup + [ / ] curry bi@
 116     [ [ v*n ] curry ] bi@ bi*  v+ ;
 117
 118 : leftof ( a b x -- ? ) dup [ lincomb ] dip [ first ] bi@ > ;
 119
 120 : ray ( a b x -- ? ) [ between ] [ leftof ] 3bi and ;
 121
 122 : raycast ( poly x -- ? )
 123     [ dup first suffix [ rest-slice ] [ but-last-slice ] bi ] dip
 124     [ ray ] curry 2map
 125     f [ xor ] reduce ;