T Pt
   Float x, y

   F (x, y)
      .x = x
      .y = y

   F String()
      R ‘Pt(x=#., y=#.)’.format(.x, .y)

T Edge
   Pt a, b

   F (a, b)
      .a = a
      .b = b

   F String()
      R ‘Edge(a=#., b=#.)’.format(.a, .b)

T Poly
   String name
   [Edge] edges

   F (name, edges)
      .name = name
      .edges = edges

V _eps = 0.00001
V _huge = 1e+100
V _tiny = 1e-100

F rayintersectseg(=p, edge)
   V a = edge.a
   V b = edge.b
   I a.y > b.y
      swap(&a, &b)
   I p.y == a.y | p.y == b.y
      p = Pt(p.x, p.y + :_eps)

   V intersect = 0B

   I (p.y > b.y | p.y < a.y) | (p.x > max(a.x, b.x))
      R 0B

   I p.x < min(a.x, b.x)
      intersect = 1B
   E
      Float m_red, m_blue
      I abs(a.x - b.x) > :_tiny
         m_red = (b.y - a.y) / Float(b.x - a.x)
      E
         m_red = :_huge
      I abs(a.x - p.x) > :_tiny
         m_blue = (p.y - a.y) / Float(p.x - a.x)
      E
         m_blue = :_huge
      intersect = m_blue >= m_red
   R intersect

F ispointinside(p, poly)
   R sum(poly.edges.map(edge -> Int(rayintersectseg(@p, edge)))) % 2 == 1

F polypp(poly)
   print("\n  Polygon(name='#.', edges=(".format(poly.name))
   print(‘    ’(poly.edges.map(e -> String(e)).join(",\n    ")"\n    ))"))

V polys = [
  Poly(name' ‘square’, edges' [Edge(Pt(0, 0), Pt(10, 0)), Edge(Pt(10, 0), Pt(10, 10)), Edge(Pt(10, 10), Pt(0, 10)), Edge(Pt(0, 10), Pt(0, 0))]),
  Poly(name' ‘square_hole’, edges' [Edge(Pt(0, 0), Pt(10, 0)), Edge(Pt(10, 0), Pt(10, 10)), Edge(Pt(10, 10), Pt(0, 10)), Edge(Pt(0, 10), Pt(0, 0)), Edge(Pt(2.5, 2.5), Pt(7.5, 2.5)), Edge(Pt(7.5, 2.5), Pt(7.5, 7.5)), Edge(Pt(7.5, 7.5), Pt(2.5, 7.5)), Edge(Pt(2.5, 7.5), Pt(2.5, 2.5))]),
  Poly(name' ‘strange’, edges' [Edge(Pt(0, 0), Pt(2.5, 2.5)), Edge(Pt(2.5, 2.5), Pt(0, 10)), Edge(Pt(0, 10), Pt(2.5, 7.5)), Edge(Pt(2.5, 7.5), Pt(7.5, 7.5)), Edge(Pt(7.5, 7.5), Pt(10, 10)), Edge(Pt(10, 10), Pt(10, 0)), Edge(Pt(10, 0), Pt(2.5, 2.5))]),
  Poly(name' ‘exagon’, edges' [Edge(Pt(3, 0), Pt(7, 0)), Edge(Pt(7, 0), Pt(10, 5)), Edge(Pt(10, 5), Pt(7, 10)), Edge(Pt(7, 10), Pt(3, 10)), Edge(Pt(3, 10), Pt(0, 5)), Edge(Pt(0, 5), Pt(3, 0))])]

V testpoints = [Pt(5, 5), Pt(5, 8),
                Pt(-10, 5), Pt(0, 5),
                Pt(10, 5), Pt(8, 5),
                Pt(10, 10)]

print("\n TESTING WHETHER POINTS ARE WITHIN POLYGONS")
L(poly) polys
   polypp(poly)
   print(‘    ’testpoints[0.<3].map(p -> ‘#.: #.’.format(p, I ispointinside(p, @poly) {‘True’} E ‘False’)).join("\t"))
   print(‘    ’testpoints[3.<6].map(p -> ‘#.: #.’.format(p, I ispointinside(p, @poly) {‘True’} E ‘False’)).join("\t"))
   print(‘    ’testpoints[6.. ].map(p -> ‘#.: #.’.format(p, I ispointinside(p, @poly) {‘True’} E ‘False’)).join("\t"))