Skip to content

Instantly share code, notes, and snippets.

Created September 24, 2020 18:33
Show Gist options
  • Save amirrajan/d92ee728dbbd8a2c556f516c28290d24 to your computer and use it in GitHub Desktop.
Save amirrajan/d92ee728dbbd8a2c556f516c28290d24 to your computer and use it in GitHub Desktop.
Winding number collision detection in Ruby.
# coding: utf-8
module Geometry
class Points
def initialize array
@array = array
def point_inside? point
return false unless point.inside_rect? rect
winding_number(point) == 0
def xs do |p|
def ys do |p|
def min_x
def max_x
def min_y
def max_y
def rect
max_x - min_x,
max_y - min_y]
# 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 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;}
# //===================================================================
# ========================================================================
# ========================================================================
# // 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
# // 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
if cn.mod_zero?(2)
return 0
return 1
# // 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]
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)
if (point_two.y <= p.y)
if (p.ray_test(point_one, point_two) == :left)
return wn
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment