amirrajan/winding_number.rb

## winding_number.rb
# coding: utf-8
module Geometry
  class Points
    def initialize array
      @array = array
    end

    def point_inside? point
      return false unless point.inside_rect? rect
      winding_number(point) == 0
    end

    def xs
      @array.map do |p|
        p.x
      end
    end

    def ys
      @array.map do |p|
        p.y
      end
    end

    def min_x
      xs.min
    end

    def max_x
      xs.max
    end

    def min_y
      ys.min
    end

    def max_y
      ys.max
    end

    def rect
      [min_x,
       min_y,
       max_x - min_x,
       max_y - min_y]
    end

    # The Winding Number
    # On the other hand, the winding number accurately determines if a
    # point is inside a nonsimple closed polygon. It does this by
    # computing how many times the polygon winds around the point. A point
    # is outside only when the polygon doesn't wind around the point at
    # all which is when the winding number wn = 0. More generally, one can
    # define the winding number wn(P,C) of any closed continuous curve C
    # around a point P in the 2D plane. Let the continuous 2D curve C be
    # defined by the points C(u)=(x(u),y(u)), for u.ge-0.le-1 with
    # C(0)=C(1). And let P be a point not on C. Then, define the vector
    # c(P,u) = C(u) – P from P to C(u), and the unit vector w(P,u) =
    # c(P,u) / |c(P,u)| which gives a continuous function W(P)=C-map-S1
    # mapping the point C(u) on C to the point w(P,u) on the unit circle
    # S1={(x,y).st.x2+y2=1}. This map can be represented in polar
    # coordinates as W(P)(u)=(cos-theta(u),sin-theta(u) where theta(u) is
    # a positive counterclockwise angle in radians. The winding number
    # wn(P,C) is then equal to the integer number of times that W(P) wraps
    # C around S1. This corresponds to a homotopy class of S1, and can be
    # computed by the integral:

    # wn_integral

    # When the curve C is a polygon with vertices V0,V1,...,Vn = V0, this
    # integral reduces to the sum of the (signed) angles that each edge
    # ViVi+1 subtends with the point P. So, if theta-i=angle(PV-i,PV-i+1),
    # we have:


    # wn1
    # wn_sum


    # This formula is clearly not very efficient since it uses a
    # computationally expensive arccos() trig function. But, a simple
    # observation lets us replace this formula by a more efficient
    # one. Pick any point Q on S1. Then, as the curve W(P) wraps around
    # S1, it passes Q a certain number of times. If we count (+1) when it
    # passes Q counterclockwise, and (–1) when it passes clockwise, then
    # the accumulated sum is exactly the total number of times that W(P)
    # wraps around S1, and is equal to the winding number
    # wn(P,C). Further, if we take an infinite ray R starting at P and
    # extending in the direction of the vector Q, then intersections where
    # R crosses the curve C correspond to the points where W(P) passes
    # Q. To do the math, we have to distinguish between positive and
    # negative crossings where C crosses R from right-to-left or
    # left-to-right. This can be determined by the sign of the dot product
    # between a normal vector to C and the direction vector q = Q [Foley
    # et al, 1996, p-965], and when the curve C is a polygon, one just
    # needs to make this determination once for each edge. For a
    # horizontal ray R from P, testing whether an edge's endpoints are
    # above and below the ray suffices. If the edge crosses the positive
    # ray from below to above, the crossing is positive (+1); but if it
    # crosses from above to below, the crossing is negative (–1). One then
    # simply adds all crossing values to get wn(P,C). For example:

    # wn2

    # Additionally, one can avoid computing the actual edge-ray
    # intersection point by using the isLeft() attribute; however, it
    # needs to be applied differently for ascending and descending
    # edges. If an upward edge crosses the ray to the right of P, then P
    # is on the left side of the edge since the triangle ViVi+1P is
    # oriented counterclockwise. On the other hand, a downward edge
    # crossing the positive ray would have P on the right side since the
    # triangle ViVi+1P would then be oriented clockwise.

    # Pseudo-Code: Winding Number Inclusion
    # This results in the following wn algorithm which is an adaptation of
    # the cn algorithm and uses the same edge crossing rules as before to
    # handle special cases.

    # typedef struct {int x, y;} Point;

    # wn_PnPoly( Point P, Point V[], int n )
    # {
    #     int    wn = 0;    // the  winding number counter

    #     // loop through all edges of the polygon
    #     for (each edge E[i]:V[i]V[i+1] of the polygon) {
    #         if (E[i] crosses upward ala Rule #1)  {
    #             if (P is  strictly left of E[i])    // Rule #4
    #                  ++wn;   // a valid up intersect right of P.x
    #         }
    #         else
    #         if (E[i] crosses downward ala Rule  #2) {
    #             if (P is  strictly right of E[i])   // Rule #4
    #                  --wn;   // a valid down intersect right of P.x
    #         }
    #     }
    #     return wn;    // =0 <=> P is outside the polygon

    # }

    # Clearly, this winding number algorithm has the same efficiency
    # as the analogous crossing number algorithm. Thus, since it is
    # more accurate in general, the winding number algorithm should
    # always be the preferred method to determine inclusion of a point
    # in an arbitrary polygon.

    # The wn algorithm's efficiency can be improved further by
    # rearranging the crossing comparison tests. This is shown in the
    # detailed implementation of wn_PnPoly() given below. In that
    # code, all edges that are totally above or totally below P get
    # rejected after only two (2) inequality tests. However, currently
    # popular implementations of the cn algorithm ([Franklin, 2000], [
    # Haines, 1994], [O'Rourke, 1998]) use at least three (3)
    # inequality tests for each rejected edge. Since most of the edges
    # in a large polygon get rejected in practical applications, there
    # is about a 33% (or more) reduction in the number of comparisons
    # done. In runtime tests using very large (1,000,000 edge) random
    # polygons (with edge length < 1/10 the polygon diameter) and 1000
    # random test points (inside the polygon's bounding box), we
    # measured a 20% increase in efficiency overall.


    # Enhancements
    # There are some enhancements to point in polygon algorithms
    # [Haines, 1994] that software developers should be aware of. We
    # mention a few that pertain to ray crossing algorithms. However,
    # there are other techniques that give better performance in
    # special cases such as testing inclusion in small convex polygons
    # like triangles. These are discussed in [Haines, 1994].


    # Bounding Box or Ball
    # It is efficient to first test that a point P is inside the
    # bounding box or ball of a large polygon before testing all edges
    # for ray crossings. If a point is outside the bounding box or
    # ball, it is also outside the polygon, and no further testing is
    # needed. But, one must precompute the bounding box (the max and
    # min for vertex x and y coordinates) or the bounding ball (center
    # and minimum radius) and store it for future use. This is worth
    # doing if more than a few points are going to be tested for
    # inclusion, which is generally the case. Further information
    # about computing bounding containers can be found in Algorithm
    # 8.

    # 3D Planar Polygons
    # In 3D applications, one sometimes wants to test a point and
    # polygon that are in the same plane. For example, one may have
    # the intersection point of a ray with the plane of a polyhedron's
    # face, and want to test if it is inside the face. Or one may want
    # to know if the base of a 3D perpendicular dropped from a point
    # is inside a planar polygon.

    # 3D inclusion is easily determined by projecting the point and
    # polygon into 2D. To do this, one simply ignores one of the 3D
    # coordinates and uses the other two. To optimally select the
    # coordinate to ignore, compute a normal vector to the plane, and
    # select the coordinate with the largest absolute value [Snyder &
    # Barr, 1987]. This gives the projection of the polygon with
    # maximum area, and results in robust computations.

    # Convex or Monotone Polygons
    # When a polygon is known to be convex, many computations can be
    # speeded up. For example, the bounding box can be computed in
    # O(log n) time [O'Rourke, 1998] instead of O(n) as for nonconvex
    # polygons. Also, point inclusion can be tested in O(log n)
    # time. More generally, the following algorithm works for convex
    # or monotone polygons.

    # A convex or y-monotone polygon can be split into two polylines,
    # one with edges increasing in y, and one with edges decreasing in
    # y. The split occurs at two vertices with max y and min y
    # coordinates. Note that these vertices are at the top and bottom
    # of the polygon's bounding box, and if they are both above or
    # both below P’s y-coordinate, then the test point is outside the
    # polygon. Otherwise, each of these two polylines intersects a ray
    # paralell to the x-axis once, and each potential crossing edge
    # can be found by doing a binary search on the vertices of each
    # polyline. This results in a practical O(log n) (preprocessing
    # and runtime) algorithm for convex and montone polygon inclusion
    # testing. For example, in the following diagram, only three
    # binary search vertices have to be tested on each polyline (A1,
    # A2, A3 ascending; and B1, B2, B3 descending):

    # This method also works for polygons that are monotone in an
    # arbitrary direction. One then uses a ray from P that is
    # perpendicular to the direction of monotonicity, and the
    # algorithm is easily adapted. A little more computation is needed
    # to determine which side of the ray a vertex is on, but for a
    # large enough polygon, the O(log n) performance more than makes
    # up for the overhead.

    # // Copyright 2000 softSurfer, 2012 Dan Sunday
    # // This code may be freely used and modified for any purpose
    # // providing that this copyright notice is included with it.
    # // SoftSurfer makes no warranty for this code, and cannot be held
    # // liable for any real or imagined damage resulting from its use.
    # // Users of this code must verify correctness for their application.


    # // a Point is defined by its coordinates {int x, y;}
    # //===================================================================


    # ========================================================================
    # DONE
    # ========================================================================
    # // isLeft(): tests if a point is Left|On|Right of an infinite line.
    # //    Input:  three points P0, P1, and P2
    # //    Return: >0 for P2 left of the line through P0 and P1
    # //            =0 for P2  on the line
    # //            <0 for P2  right of the line
    # //    See: Algorithm 1 "Area of Triangles and Polygons"
    # inline int
    # isLeft( Point P0, Point P1, Point P2 )
    # {
    #     return ( (P1.x - P0.x) * (P2.y - P0.y)
    #             - (P2.x -  P0.x) * (P1.y - P0.y) );
    # }
    # //===================================================================
    def ray_test p0, p1, p2
      p2.ray_test p0, p1
    end


    # // cn_PnPoly(): crossing number test for a point in a polygon
    # //      Input:   P = a point,
    # //               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
    # //      Return:  0 = outside, 1 = inside
    # // This code is patterned after [Franklin, 2000]
    # int
    # cn_PnPoly( Point P, Point* V, int n )
    # {
    #     int    cn = 0;    // the  crossing number counter

    #     // loop through all edges of the polygon
    #     for (int i=0; i<n; i++) {    // edge from V[i]  to V[i+1]
    #        if (((V[i].y <= P.y) && (V[i+1].y > P.y))     // an upward crossing
    #         || ((V[i].y > P.y) && (V[i+1].y <=  P.y))) { // a downward crossing
    #             // compute  the actual edge-ray intersect x-coordinate
    #             float vt = (float)(P.y  - V[i].y) / (V[i+1].y - V[i].y);
    #             if (P.x <  V[i].x + vt * (V[i+1].x - V[i].x)) // P.x < intersect
    #                  ++cn;   // a valid crossing of y=P.y right of P.x
    #         }
    #     }
    #     return (cn&1);    // 0 if even (out), and 1 if  odd (in)

    # }
    # //===================================================================
    def crossing_number point
      cn = 0
      v  = @array
      p  = point
      v.each_with_index do |p, i|
        if i != v.length - 1
          if (((v[i].y <= p.y) && (v[i+1].y > p.y)) ||
              ((v[i].y > p.y) && (v[i+1].y <= p.y)))
            vt = (p.y - v[i].y).fdiv(v[i+1].y - v[i].y)
            if (p.x < v[i].x + vt * (v[i+1].x - v[i].x))
              cn += 1
            end
          end
        end
      end

      if cn.mod_zero?(2)
        return 0
      else
        return 1
      end
    end

    # // wn_PnPoly(): winding number test for a point in a polygon
    # //      Input:   P = a point,
    # //               V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
    # //      Return:  wn = the winding number (=0 only when P is outside)
    # int
    # wn_PnPoly( Point P, Point* V, int n )
    # {
    #     int    wn = 0;    // the  winding number counter

    #     // loop through all edges of the polygon
    #     for (int i=0; i<n; i++) {   // edge from V[i] to  V[i+1]
    #         if (V[i].y <= P.y) {          // start y <= P.y
    #             if (V[i+1].y  > P.y)      // an upward crossing
    #                  if (isLeft( V[i], V[i+1], P) > 0)  // P left of  edge
    #                      ++wn;            // have  a valid up intersect
    #         }
    #         else {                        // start y > P.y (no test needed)
    #             if (V[i+1].y  <= P.y)     // a downward crossing
    #                  if (isLeft( V[i], V[i+1], P) < 0)  // P right of  edge
    #                      --wn;            // have  a valid down intersect
    #         }
    #     }
    #     return wn;
    # }
    # //===================================================================
    def winding_number point
      wn = 0
      v = @array
      p = point
      v.each_with_index do |_, i|
        point_one = v[i]
        point_two = v[i+1]
        if !point_two
          point_two = v[0]
        end

        positive_winding = :right
        negative_winding = :left
        if (point_one.y <= p.y)
          if (point_two.y > p.y)
            if (p.ray_test(point_one, point_two) == :right)
              wn+=1
            end
          end
        else
          if (point_two.y <= p.y)
            if (p.ray_test(point_one, point_two) == :left)
              wn-=1
            end
          end
        end
      end

      return wn
    end
  end
end
	# coding: utf-8
	module Geometry
	class Points
	def initialize array
	@array = array
	end

	def point_inside? point
	return false unless point.inside_rect? rect
	winding_number(point) == 0
	end

	def xs
	@array.map do \|p\|
	p.x
	end
	end

	def ys
	@array.map do \|p\|
	p.y
	end
	end

	def min_x
	xs.min
	end

	def max_x
	xs.max
	end

	def min_y
	ys.min
	end

	def max_y
	ys.max
	end

	def rect
	[min_x,
	min_y,
	max_x - min_x,
	max_y - min_y]
	end

	# The Winding Number
	# On the other hand, the winding number accurately determines if a
	# point is inside a nonsimple closed polygon. It does this by
	# computing how many times the polygon winds around the point. A point
	# is outside only when the polygon doesn't wind around the point at
	# all which is when the winding number wn = 0. More generally, one can
	# define the winding number wn(P,C) of any closed continuous curve C
	# around a point P in the 2D plane. Let the continuous 2D curve C be
	# defined by the points C(u)=(x(u),y(u)), for u.ge-0.le-1 with
	# C(0)=C(1). And let P be a point not on C. Then, define the vector
	# c(P,u) = C(u) – P from P to C(u), and the unit vector w(P,u) =
	# c(P,u) / \|c(P,u)\| which gives a continuous function W(P)=C-map-S1
	# mapping the point C(u) on C to the point w(P,u) on the unit circle
	# S1={(x,y).st.x2+y2=1}. This map can be represented in polar
	# coordinates as W(P)(u)=(cos-theta(u),sin-theta(u) where theta(u) is
	# a positive counterclockwise angle in radians. The winding number
	# wn(P,C) is then equal to the integer number of times that W(P) wraps
	# C around S1. This corresponds to a homotopy class of S1, and can be
	# computed by the integral:

	# wn_integral

	# When the curve C is a polygon with vertices V0,V1,...,Vn = V0, this
	# integral reduces to the sum of the (signed) angles that each edge
	# ViVi+1 subtends with the point P. So, if theta-i=angle(PV-i,PV-i+1),
	# we have:


	# wn1
	# wn_sum


	# This formula is clearly not very efficient since it uses a
	# computationally expensive arccos() trig function. But, a simple
	# observation lets us replace this formula by a more efficient
	# one. Pick any point Q on S1. Then, as the curve W(P) wraps around
	# S1, it passes Q a certain number of times. If we count (+1) when it
	# passes Q counterclockwise, and (–1) when it passes clockwise, then
	# the accumulated sum is exactly the total number of times that W(P)
	# wraps around S1, and is equal to the winding number
	# wn(P,C). Further, if we take an infinite ray R starting at P and
	# extending in the direction of the vector Q, then intersections where
	# R crosses the curve C correspond to the points where W(P) passes
	# Q. To do the math, we have to distinguish between positive and
	# negative crossings where C crosses R from right-to-left or
	# left-to-right. This can be determined by the sign of the dot product
	# between a normal vector to C and the direction vector q = Q [Foley
	# et al, 1996, p-965], and when the curve C is a polygon, one just
	# needs to make this determination once for each edge. For a
	# horizontal ray R from P, testing whether an edge's endpoints are
	# above and below the ray suffices. If the edge crosses the positive
	# ray from below to above, the crossing is positive (+1); but if it
	# crosses from above to below, the crossing is negative (–1). One then
	# simply adds all crossing values to get wn(P,C). For example:

	# wn2

	# Additionally, one can avoid computing the actual edge-ray
	# intersection point by using the isLeft() attribute; however, it
	# needs to be applied differently for ascending and descending
	# edges. If an upward edge crosses the ray to the right of P, then P
	# is on the left side of the edge since the triangle ViVi+1P is
	# oriented counterclockwise. On the other hand, a downward edge
	# crossing the positive ray would have P on the right side since the
	# triangle ViVi+1P would then be oriented clockwise.

	# Pseudo-Code: Winding Number Inclusion
	# This results in the following wn algorithm which is an adaptation of
	# the cn algorithm and uses the same edge crossing rules as before to
	# handle special cases.

	# typedef struct {int x, y;} Point;

	# wn_PnPoly( Point P, Point V[], int n )
	# {
	# int wn = 0; // the winding number counter

	# // loop through all edges of the polygon
	# for (each edge E[i]:V[i]V[i+1] of the polygon) {
	# if (E[i] crosses upward ala Rule #1) {
	# if (P is strictly left of E[i]) // Rule #4
	# ++wn; // a valid up intersect right of P.x
	# }
	# else
	# if (E[i] crosses downward ala Rule #2) {
	# if (P is strictly right of E[i]) // Rule #4
	# --wn; // a valid down intersect right of P.x
	# }
	# }
	# return wn; // =0 <=> P is outside the polygon

	# }

	# Clearly, this winding number algorithm has the same efficiency
	# as the analogous crossing number algorithm. Thus, since it is
	# more accurate in general, the winding number algorithm should
	# always be the preferred method to determine inclusion of a point
	# in an arbitrary polygon.

	# The wn algorithm's efficiency can be improved further by
	# rearranging the crossing comparison tests. This is shown in the
	# detailed implementation of wn_PnPoly() given below. In that
	# code, all edges that are totally above or totally below P get
	# rejected after only two (2) inequality tests. However, currently
	# popular implementations of the cn algorithm ([Franklin, 2000], [
	# Haines, 1994], [O'Rourke, 1998]) use at least three (3)
	# inequality tests for each rejected edge. Since most of the edges
	# in a large polygon get rejected in practical applications, there
	# is about a 33% (or more) reduction in the number of comparisons
	# done. In runtime tests using very large (1,000,000 edge) random
	# polygons (with edge length < 1/10 the polygon diameter) and 1000
	# random test points (inside the polygon's bounding box), we
	# measured a 20% increase in efficiency overall.


	# Enhancements
	# There are some enhancements to point in polygon algorithms
	# [Haines, 1994] that software developers should be aware of. We
	# mention a few that pertain to ray crossing algorithms. However,
	# there are other techniques that give better performance in
	# special cases such as testing inclusion in small convex polygons
	# like triangles. These are discussed in [Haines, 1994].


	# Bounding Box or Ball
	# It is efficient to first test that a point P is inside the
	# bounding box or ball of a large polygon before testing all edges
	# for ray crossings. If a point is outside the bounding box or
	# ball, it is also outside the polygon, and no further testing is
	# needed. But, one must precompute the bounding box (the max and
	# min for vertex x and y coordinates) or the bounding ball (center
	# and minimum radius) and store it for future use. This is worth
	# doing if more than a few points are going to be tested for
	# inclusion, which is generally the case. Further information
	# about computing bounding containers can be found in Algorithm
	# 8.

	# 3D Planar Polygons
	# In 3D applications, one sometimes wants to test a point and
	# polygon that are in the same plane. For example, one may have
	# the intersection point of a ray with the plane of a polyhedron's
	# face, and want to test if it is inside the face. Or one may want
	# to know if the base of a 3D perpendicular dropped from a point
	# is inside a planar polygon.

	# 3D inclusion is easily determined by projecting the point and
	# polygon into 2D. To do this, one simply ignores one of the 3D
	# coordinates and uses the other two. To optimally select the
	# coordinate to ignore, compute a normal vector to the plane, and
	# select the coordinate with the largest absolute value [Snyder &
	# Barr, 1987]. This gives the projection of the polygon with
	# maximum area, and results in robust computations.

	# Convex or Monotone Polygons
	# When a polygon is known to be convex, many computations can be
	# speeded up. For example, the bounding box can be computed in
	# O(log n) time [O'Rourke, 1998] instead of O(n) as for nonconvex
	# polygons. Also, point inclusion can be tested in O(log n)
	# time. More generally, the following algorithm works for convex
	# or monotone polygons.

	# A convex or y-monotone polygon can be split into two polylines,
	# one with edges increasing in y, and one with edges decreasing in
	# y. The split occurs at two vertices with max y and min y
	# coordinates. Note that these vertices are at the top and bottom
	# of the polygon's bounding box, and if they are both above or
	# both below P’s y-coordinate, then the test point is outside the
	# polygon. Otherwise, each of these two polylines intersects a ray
	# paralell to the x-axis once, and each potential crossing edge
	# can be found by doing a binary search on the vertices of each
	# polyline. This results in a practical O(log n) (preprocessing
	# and runtime) algorithm for convex and montone polygon inclusion
	# testing. For example, in the following diagram, only three
	# binary search vertices have to be tested on each polyline (A1,
	# A2, A3 ascending; and B1, B2, B3 descending):

	# This method also works for polygons that are monotone in an
	# arbitrary direction. One then uses a ray from P that is
	# perpendicular to the direction of monotonicity, and the
	# algorithm is easily adapted. A little more computation is needed
	# to determine which side of the ray a vertex is on, but for a
	# large enough polygon, the O(log n) performance more than makes
	# up for the overhead.

	# // Copyright 2000 softSurfer, 2012 Dan Sunday
	# // This code may be freely used and modified for any purpose
	# // providing that this copyright notice is included with it.
	# // SoftSurfer makes no warranty for this code, and cannot be held
	# // liable for any real or imagined damage resulting from its use.
	# // Users of this code must verify correctness for their application.


	# // a Point is defined by its coordinates {int x, y;}
	# //===================================================================


	# ========================================================================
	# DONE
	# ========================================================================
	# // isLeft(): tests if a point is Left\|On\|Right of an infinite line.
	# // Input: three points P0, P1, and P2
	# // Return: >0 for P2 left of the line through P0 and P1
	# // =0 for P2 on the line
	# // <0 for P2 right of the line
	# // See: Algorithm 1 "Area of Triangles and Polygons"
	# inline int
	# isLeft( Point P0, Point P1, Point P2 )
	# {
	# return ( (P1.x - P0.x) * (P2.y - P0.y)
	# - (P2.x - P0.x) * (P1.y - P0.y) );
	# }
	# //===================================================================
	def ray_test p0, p1, p2
	p2.ray_test p0, p1
	end


	# // cn_PnPoly(): crossing number test for a point in a polygon
	# // Input: P = a point,
	# // V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
	# // Return: 0 = outside, 1 = inside
	# // This code is patterned after [Franklin, 2000]
	# int
	# cn_PnPoly( Point P, Point* V, int n )
	# {
	# int cn = 0; // the crossing number counter

	# // loop through all edges of the polygon
	# for (int i=0; i<n; i++) { // edge from V[i] to V[i+1]
	# if (((V[i].y <= P.y) && (V[i+1].y > P.y)) // an upward crossing
	# \|\| ((V[i].y > P.y) && (V[i+1].y <= P.y))) { // a downward crossing
	# // compute the actual edge-ray intersect x-coordinate
	# float vt = (float)(P.y - V[i].y) / (V[i+1].y - V[i].y);
	# if (P.x < V[i].x + vt * (V[i+1].x - V[i].x)) // P.x < intersect
	# ++cn; // a valid crossing of y=P.y right of P.x
	# }
	# }
	# return (cn&1); // 0 if even (out), and 1 if odd (in)

	# }
	# //===================================================================
	def crossing_number point
	cn = 0
	v = @array
	p = point
	v.each_with_index do \|p, i\|
	if i != v.length - 1
	if (((v[i].y <= p.y) && (v[i+1].y > p.y)) \|\|
	((v[i].y > p.y) && (v[i+1].y <= p.y)))
	vt = (p.y - v[i].y).fdiv(v[i+1].y - v[i].y)
	if (p.x < v[i].x + vt * (v[i+1].x - v[i].x))
	cn += 1
	end
	end
	end
	end

	if cn.mod_zero?(2)
	return 0
	else
	return 1
	end
	end

	# // wn_PnPoly(): winding number test for a point in a polygon
	# // Input: P = a point,
	# // V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
	# // Return: wn = the winding number (=0 only when P is outside)
	# int
	# wn_PnPoly( Point P, Point* V, int n )
	# {
	# int wn = 0; // the winding number counter

	# // loop through all edges of the polygon
	# for (int i=0; i<n; i++) { // edge from V[i] to V[i+1]
	# if (V[i].y <= P.y) { // start y <= P.y
	# if (V[i+1].y > P.y) // an upward crossing
	# if (isLeft( V[i], V[i+1], P) > 0) // P left of edge
	# ++wn; // have a valid up intersect
	# }
	# else { // start y > P.y (no test needed)
	# if (V[i+1].y <= P.y) // a downward crossing
	# if (isLeft( V[i], V[i+1], P) < 0) // P right of edge
	# --wn; // have a valid down intersect
	# }
	# }
	# return wn;
	# }
	# //===================================================================
	def winding_number point
	wn = 0
	v = @array
	p = point
	v.each_with_index do \|_, i\|
	point_one = v[i]
	point_two = v[i+1]
	if !point_two
	point_two = v[0]
	end

	positive_winding = :right
	negative_winding = :left
	if (point_one.y <= p.y)
	if (point_two.y > p.y)
	if (p.ray_test(point_one, point_two) == :right)
	wn+=1
	end
	end
	else
	if (point_two.y <= p.y)
	if (p.ray_test(point_one, point_two) == :left)
	wn-=1
	end
	end
	end
	end

	return wn
	end
	end
	end