Programming Challenge: 2D Geometry

This blog publishes each week a programming challenge they think would appropriate for a programming job interview. In my opinion most of them aren’t very well suited for a job interview but this isn’t the subject of my post. The aim of my post is to solve this week challenge.

This week question is:

Your input:

1. List of points in a 2D space. If you draw lines between one point to the next one, a closed polygon is created (can be either concave or convex).
2. A single point in a 2D space.

You need to determine whether the given point is inside or outside the given polygon.

[...]

Provide the most efficient and simple algorithm.

The “simplest” solution

The simplest algorithm I have in mind would be:

• draw red point
• draw blue filled polygon
• check point color (red means outside, blue means inside)

Problem solved is a few lines of code . But what happens if coordinates aren’t integers or are really big (> 10000) ? This algorithm is not very accurate, we’d run out of memory and filling up a big polygon would take lots of time.

The geometrical solution

Lets draw a ray starting at our point in a random direction to the infinity. This ray intersects our polygon zero or more times.

Zero times would mean the point is outside the polygon. One time means the point is inside, two times outside and so on… Even number of intersections means the point is outside, odd number means the point is inside. How do I prove that ?

I don’t have a mathematical proof, rather an intuitive one… Lets choose an imaginary point on the ray at infinity. Could this point be inside the polygon ? No, it can’t: because the polygon is closed and finite a point at infinity can’t be inside it. Lets walk backwards from the point to infinity to our reference point. Each time the ray intersects the polygon we switch from inside to outside and vice versa. Even number of switches (intersections) -> the point is outside; odd number of switches -> the point is inside. The algorithm is pretty simple and is O(n) complexity:

• Get a ray originating with the input point
• Transform the input polygon points into segments
• For each segment check if it intersects the ray and count the number of intersections
• If odd number of intersections the point is inside else it is outside

How can we implement that ?! The mathematical representation of a ray is:

Ray(t) = Po + t*Pd with t>=0

Where Po is the origin point and Pd is a point somewhere on the ray witch gives the ray direction.

The mathematical representation of a segment is:

Segment(u) = P1 + u * (P2-P1) with 0<=u<=1

Where P1 and P2 are the head points of the segment. If we replace (P2-P1) by D (as in Delta) we get the same formula as the ray. Only the limits for the argument are different.

To find out if the ray intersects the segment we have to solve the equation

Ray(t) = Segment(u)

and check if t>=0 and 0<=u<=1. But I won’t steal you the pleasure of solving it

Show me the code

First let me warn you that for the sake of simplicity I eliminated input validation code: polygon should have at least three points, segment and ray validation (origin and direction points should be different).

The language used is C# 3.0.

Lets start with some simple geometric classes:

public class Point
{
public Point(float x, float y)
{
X = x;
Y = y;
}

public float X { get; private set; }
public float Y { get; private set; }
}

public class Ray
{
public Ray(Point p1, Point p2)
{
Origin = p1;
Direction = p2;
}

public Point Origin { get; private set; }
public Point Direction { get; private set; }
}

public class Segment
{
public Segment(Point p1, Point p2)
{
Origin = p1;
Delta = new Point(p2.X – p1.X, p2.Y – p1.Y);
}

public Point Origin { get; private set; }
public Point Delta { get; private set; }
}

Now the method witch finds out if a ray and an segment intersect:

///

/// Solves the equation Ray(t) = Segment(u) where
/// Ray(t) = Po + t*Pd with t>=0
/// Segment(u) = P1 + u * (P2-P1) with 0&lt=u&lt=1
///

/// the ray /// the segment /// true if ray and segment intersect
private static bool Intersecting(Ray ray, Segment seg)
{
var a =
seg.Delta.X * ray.Direction.Y -
seg.Delta.Y * ray.Direction.X;

// if ray and segment are parallel they don’t intersect
if (0 == a) return false;

var b =
(ray.Origin.X – seg.Origin.X) * ray.Direction.Y -
(ray.Origin.Y – seg.Origin.Y) * ray.Direction.X;

var u = b / a;

b =
(ray.Origin.X – seg.Origin.X) * seg.Delta.Y -
(ray.Origin.Y – seg.Origin.Y) * seg.Delta.X;

var t = b / a;

return 0 <= t && 0 <= u && 1 > u;
}
} 

Did you notice that something changed ? The constraints on u changed from 0<=u<=1 to 0<=u<1. What happens if the ray contains one of the polygon points ? The intersection would be counted twice ! Hence the change.

To define the polygon the input is a list of points but our algorithm needs segments so we transform the input:

///

/// Transforms a list of polygon points in a list of polygon segments
///

/// list of polygon points /// list of polygon segments
private static IList GetPolySegments(IList polyPoints)
{
IList segments = new List();

for (var i = 0; i < polyPoints.Count - 1; i++)
{
segments.Add(new Segment(polyPoints[i], polyPoints[i + 1]));
}

segments.Add(new Segment(polyPoints[polyPoints.Count - 1], polyPoints[0]));

return segments;
}

To define the ray we need an origin point (the input point) and a direction point. For fun I choosed the direction point to be the center of the polygon.

///

/// Calculates the center of a polygon. (the averege of coordinates)
///

/// collection of polygon /// Center Point
private static Point GetPolyCenter(ICollection polyPoints)
{
return new Point(
polyPoints.Sum(p => p.X) / polyPoints.Count,
polyPoints.Sum(p => p.Y) / polyPoints.Count);
}

Any random direction point would work as long as it is different form the input point.

And now the main algorithm:

///

/// Finds if a point is inside a closed polygon
///

/// point /// list of points definig the polygon /// true if point p is inside the polygon
public static bool IsPointInsidePoly(Point p, IList polyPoints)
{
// define ray
var rayToPolyCenter = new Ray(p, GetPolyCenter(polyPoints));

// count intersections
var intersections =  GetPolySegments(polyPoints).Count(
s => Intersecting(rayToPolyCenter, s));

// Odd number of intersections means inside
return 1 == intersections % 2;
}

So that’s it. There are some ways of improving the code (input validation, test it with very close points, etc.) but I’ll stop here.