Click here to Skip to main content
13,147,956 members (88,608 online)
Rate this:
 
Please Sign up or sign in to vote.
Today's challenge is late but simple one.

Given a set of (x,y) points representing a polygon, determine whether a given (x,y) point is inside or outside the polygon.

For example, the polygon defined by the points (counter-clockwise) { (2,0), (4,1), (4,4), (2,5), (1,2) and (2,0) } contains the point (3,3) but does not contain the point (5,4).

You can make the data type whatever you want. Double precision floating point, integer, cubits - it doesn't matter. You can use whatever language you wish to. You just need to be able to supply the input data and get the output to the question "is this point in the polygon".

There are a couple of articles here that may help (Computational Geometry, C++ and Wykobi[^] and Classes for computational geometry[^])


Last week's challenge was won (and that's a loose term) by Peter Leow[^], with special mentions for the C64 implementation[^] by Kornfeld Eliyahu Peter and the LINQ implementation[^] by Maciej Los. Peter's explanation and over-engineered approach amused the judge immensely. Peter - send us your details and a CodeProject mug is on it's way.

Apologies for the lateness of this week's challenge. Toronto winter has caused mayhem in the office.

What I have tried:

A new snow shovel - but that didn't help much
Posted 16-Dec-16 3:56am
Updated 20-Dec-16 10:41am
v2
Comments
Kornfeld Eliyahu Peter 18-Dec-16 2:40am
   
Try this: https://www.youtube.com/watch?v=IfpDDqyJnSk
ppolymorphe 20-Dec-16 23:03pm
   
Question: a point on the edge of the polygon is considered inside or outside ?
Kornfeld Eliyahu Peter 22-Dec-16 11:52am
   
If a point on the edge is not part of the polygon, than a polygon can not be defined by edge points, but only by all the points are inside the polygon...
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 2

Adding another example using most of what OriginalGriff wrote except no foreach loop.

public class TestPoints {
      GraphicsPath path;
      public TestPoints(Point[] points) {
          path = new GraphicsPath(points);
      }
 
      public bool CheckVisible(int x, int y) {
          if (null == path) {
               return false;
          }
 
          Region region = new Region(path);
          bool visible = region.IsVisible(x, y);
          string output = string.Format("Point ({0}, {1}) is {2} the polygon.", x, y, visible ? "within" : "outside");
          Console.WriteLine(output);
          return visible;
      }
}
 
public class Program {
     public static void Main(string[] args) {
         Point[] points1 = new Point[] {
               new Point(2, 0), new Point(4, 1), new Point(4, 4),
               new Point(2, 5), new Point(1, 2), new Point(2, 0)
         };
 
         Point[] points2 = new Point[] {
               new Point(2, 0), new Point(44, 1), new Point(44, 46),
               new Point(24, 57), new Point(10, 26), new Point(2, 0)
         };
 
         TestPoints test = new TestPoints(points1);
         test.CheckVisible(3, 3);
         test.CheckVisible(5, 4);
 
         test = new TestPoints(points2);
         test.CheckVisible(35, 40);
         test.CheckVisible(45, 43);
     }
}


Output:
Point (3, 3) is within the polygon.
Point (5, 4) is outside the polygon.
Point (35, 40) is within the polygon.
Point (45, 43) is outside the polygon.
  Permalink  
v2
Comments
Kornfeld Eliyahu Peter 18-Dec-16 2:43am
   
Point (35, 40) is within the polygon. - Are you sure???
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 3

Time for some maths! My code uses the ray casting algorithm: it looks at how many times a ray, starting from the given point to whatever direction, intersects with the polygon. If this number is odd, the point is inside the polygon. If it's even, the number is outside the polygon.

/* Usage:
Polygon polygon = new Polygon(new Point(2, 0), new Point(4, 1), new Point(4, 4), new Point(2, 5), new Point(1, 2), new Point(2, 0));
Console.WriteLine(polygon.PointInPolygon(new Point(3, 3))); // True
Console.WriteLine(polygon.PointInPolygon(new Point(5, 4))); // False
*/
 
public class Polygon
{
    Point[] _points;
    List<LineSegment> _segments;
 
    public Polygon(params Point[] points)
    {
        if (points[0] == points[points.Length - 1])
        {
            _points = points;
        }
        else
        {
            _points = points.Concat(new Point[] { points[0] }).ToArray();
            // auto-close the polygon
        }
 
        _segments = new List<LineSegment>();
        for (int i = 0; i < _points.Length; i++)
        {
            if (i != _points.Length - 1)
            {
                _segments.Add(new LineSegment(_points[i], _points[i + 1]));
            }
        }
    }
 
    public bool PointInPolygon(Point point)
    {
        Ray ray = new Ray(point, new Point(point.X + 1, point.Y + 1));
        int counter = 0;
        foreach (LineSegment segment in _segments)
        {
            if (segment.IntersectsWithRay(ray))
            {
                counter++;
            }
        }
        return counter % 2 == 1;
    }
}
public class Point
{
    public double X { get; private set; }
    public double Y { get; private set; }
    public Point(double x, double y)
    {
        X = x;
        Y = y;
    }
}
public class LineSegment
{
    public Point Point1 { get; private set; }
    public Point Point2 { get; private set; }
 
    public LineSegment(Point p1, Point p2)
    {
        Point1 = p1;
        Point2 = p2;
    }
 
    public bool IntersectsWithRay(Ray ray)
    {
        Line raySupport = new Line(ray.Origin, ray.AlsoGoesThrough);
        Line segmentSupport = new Line(Point1, Point2);
 
        Point intersection = raySupport.IntersectionWithOtherLine(segmentSupport);
        if (intersection == null) return false;
 
        bool intersectionOnRay = !((ray.Direction.Item1 == 1 && intersection.Y < ray.Origin.Y) || (ray.Direction.Item1 == -1 && intersection.Y > ray.Origin.Y) ||
            (ray.Direction.Item2 == 1 && intersection.X < ray.Origin.X) || (ray.Direction.Item1 == -1 && intersection.X > ray.Origin.X));
 
        double segmentMinX = Math.Min(Point1.X, Point2.X);
        double segmentMaxX = Math.Max(Point1.X, Point2.X);
        double segmentMinY = Math.Min(Point1.Y, Point2.Y);
        double segmentMaxY = Math.Max(Point1.Y, Point2.Y);
        bool intersectionOnSegment = intersection.X >= segmentMinX && intersection.X <= segmentMaxX && intersection.Y >= segmentMinY && intersection.Y <= segmentMaxY;
 
        return intersectionOnRay && intersectionOnSegment;
     }
}
 
public class Ray
{
    public Point Origin { get; private set; }
    public Point AlsoGoesThrough { get; private set; }
 
    public Tuple<int, int> Direction { get; private set; }
    // Item1 specifies up (1) or down (-1), Item2 specifies right (1) or left (-1)
 
    public Ray(Point origin, Point alsoGoesThrough)
    {
        Origin = origin;
        AlsoGoesThrough = alsoGoesThrough;
 
        int directionUp;
        if (origin.Y == alsoGoesThrough.Y)
        {
            directionUp = 0;
        }
        else if (origin.Y > alsoGoesThrough.Y)
        {
            directionUp = -1;
        }
        else // origin.Y < alsoGoesThrough.Y
        {
            directionUp = 1;
        }
 
        int directionRight;
        if (origin.X == alsoGoesThrough.X)
        {
            directionRight = 0;
        }
        else if (origin.X > alsoGoesThrough.X)
        {
            directionRight = -1;
        }
        else // origin.X < alsoGoesThrough.X
        {
            directionRight = 1;
        }
 
        Direction = new Tuple<int, int>(directionUp, directionRight);
    }
}
 
public class Line
{
    public double Slope { get; private set; }
    public double YOfIntersectionWithYAxis { get; private set; }
 
    public Line(Point goesThrough1, Point goesThrough2)
    {
        Slope = (goesThrough1.Y - goesThrough2.Y) / (goesThrough1.X - goesThrough2.X);
        YOfIntersectionWithYAxis = goesThrough1.Y - Slope * goesThrough1.X;
    }
 
    public Point IntersectionWithOtherLine(Line other)
    {
        if (Slope == other.Slope) return null;
 
        double intersectionX = (other.YOfIntersectionWithYAxis - YOfIntersectionWithYAxis) / (Slope - other.Slope);
        double intersectionY = Slope * intersectionX + YOfIntersectionWithYAxis;
        return new Point(intersectionX, intersectionY);
    }
}
  Permalink  
Comments
CDP1802 20-Dec-16 10:58am
   
And what if the point is actually outside and the ray just touches one point? It's possible.
ProgramFOX 20-Dec-16 15:00pm
   
When does that happen? The most logical response to that question would be when the ray intersects with a vertex of the polygon, but then the number of intersections is still 2! Because such a vertex is the intersection of two line segments, then IntersectsWithRay will return true twice. Am I missing another case?
CDP1802 20-Dec-16 15:16pm
   
Exactly. Geometrically it would touch only one point, but it will touch two lines of the polygon. Ok, it may work.
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 8

For JSOP only :-)
IDENTIFICATION DIVISION.
    PROGRAM-ID. CPCC4-POINT-IN-POLYGON.
DATA DIVISION.
    WORKING-STORAGE SECTION.
        01 T1 PIC S9(6) VALUE ZERO.
        01 T2 PIC S9(6) VALUE ZERO.
        01 T3 PIC S9(6) VALUE ZERO.
        01 T4 PIC S9(6) VALUE ZERO.
    LOCAL-STORAGE SECTION.
        01 POLYGON.
            02 POINT OCCURS 6 TIMES.
                03 X PIC 9(3) VALUE ZERO.
                03 Y PIC 9(3) VALUE ZERO.
        01 CHECK-POINT.
            02 CP-X PIC 9(3) VALUE ZERO.
            02 CP-Y PIC 9(3) VALUE ZERO.
        01 CNT PIC 9(1) VALUE 1.
        01 WN PIC 9(3) VALUE ZERO.
        01 SIDE PIC S9(3) VALUE ZERO.
PROCEDURE DIVISION.
    PERFORM READ-CHECK-POINT.
    PERFORM INITIALIZE-POLYGON.
    PERFORM COMPUTE-WINDING-NUMBER VARYING CNT FROM 1 BY 1 UNTIL CNT = 6.
    IF WN = 0 THEN
        DISPLAY 'THE POINT (' CP-X ',' CP-Y ') IS OUTSIDE THE POLYGON.',
    ELSE
        DISPLAY 'THE POINT (' CP-X ',' CP-Y ') IS INSIDE THE POLYGON.',
    END-IF.
    STOP RUN.
    
    COMPUTE-WINDING-NUMBER.
        PERFORM CHECK-SIDE.
 
        IF Y(CNT) <= CP-Y THEN
			IF Y(CNT + 1) > CP-Y AND SIDE > 0 THEN
			    ADD 1 TO WN
			 END-IF
		ELSE
		    IF Y(CNT + 1) <= CP-Y AND SIDE < 0 THEN
		        SUBTRACT 1 FROM WN
		    END-IF
        END-IF.
 
    CHECK-SIDE.
        SUBTRACT X(CNT) FROM X(CNT + 1) GIVING T1. 
        SUBTRACT Y(CNT) FROM CP-Y GIVING T2.
        MULTIPLY T1 BY T2 GIVING T1.
        
        SUBTRACT Y(CNT) FROM Y(CNT + 1) GIVING T3.
        SUBTRACT X(CNT) FROM CP-X GIVING T4.
        MULTIPLY T3 BY T4 GIVING T3.
        
        SUBTRACT T3 FROM T1 GIVING SIDE.
    
    READ-CHECK-POINT.
        DISPLAY 'ENTER X COORDINATE:'.
        ACCEPT CP-X.
        DISPLAY 'ENTER Y COORDINATE:'.
        ACCEPT CP-Y.
 
    INITIALIZE-POLYGON.
        MOVE 2 TO X(1).
        MOVE 0 TO Y(1).
 
        MOVE 4 TO X(2).
        MOVE 1 TO Y(2).
 
        MOVE 4 TO X(3).
        MOVE 4 TO Y(3).
 
        MOVE 2 TO X(4).
        MOVE 5 TO Y(4).
 
        MOVE 1 TO X(5).
        MOVE 2 TO Y(5).
 
        MOVE 2 TO X(6).
        MOVE 0 TO Y(6).

And a bit for the math...
It is a combination of the winding and the crossing method... From the crossing method I borrowed the idea of the ray, while from the winding method the idea of direction...
1. Uses an imaginary horizontal line that contains the point we looking for as ray
2. Checks only edges crossed by that line
3. Uses edges as vectors to see if we go CW or CCW
4. Checks on which side the point of the vector (edge) - and that side is relative to the direction...
5. Increments (left side) or decrements (right side) the winding number
6. Zero is outside
  Permalink  
v3
Comments
grgran 19-Dec-16 10:44am
   
Thumbs up for doing this in COBOL :-) !
Chris Maunder 6-Jan-17 17:50pm
   
This impressed me so much I added COBOL syntax colouring (missed LOCAL-STORAGE - will add next time)
Kornfeld Eliyahu Peter 7-Jan-17 12:48pm
   
:-) But don't be disappointed if that will not cause a immense flow of COBOL code samples...
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 11

First, I'd like to explain my methodology. Second, I'll explain my code.

So first, my technique is similar to the ray-casting algorithm. I double cast a ray in both directions and count the intersections individually then OR them together which allows me to catch points on sides like (4,1),(4,4). My technique for doing so is home-brewed. So let's get to it:

A line is represented by $y = mx + b$. So to calculate the lines that can intersect we test versus the y-values of the test point and the test lines. We test both as $> y$ in order to avoid multiple tests on connecting lines. Now we're posed with an interesting question. We can either A) calculate b by plugging in one of the existing points, or B) nullify that calculation with an axis transformation.

I went with option B. To explain: In order to calculate b you need to create a boundary for the problem. In the standard xy plane this boundary can be relative necessitating the calculation of b manually. But if you transform the points by rotating them 90 degrees, you're able to bound the problem by the y-values (which are already bounded by the ray-cast). This means you can set b equal to one of the y-bounds already determined. Crude MS Paint demonstration.[^]

This preserves the relative position of the points while enabling our calculations to be easier. So now we determine intersection with a line based on half-planes[^]. A ray casted to the right calculates the left half-plane and a ray casted to the left the right half-plane. This can be represented as $x < ((x_1 - x_2)/(y_1 - y_2))*(y - y_2 ) + x_2$ and $x > ((x_1 - x_2)/(y_1 - y_2))*(y - y_2 ) + x_2$ respectively.

My solution, organized and commented, is as follows:
>>,>>>,>>,,<<<,>>,,>>>>>>,                                      Read test point and first two polygon points
[                                                               Until end of input
  [<<<<+>+>>>-]<<<[>>>+<<<-]                                    t0 = y2
  <<<<[>>>>+>+<<<<<-]>>>>>[<<<<<+>>>>>-]                        t1 = y
  <<[>>+<[>[-]>+<<-]>>[<<+>>-]<[<<[-]+<<<<<+>>>>>>>-]<-<-]>[-]  s2 = y2 greater than y
  <<<<<[>>>>+>+<<<<<-]>>>>>[<<<<<+>>>>>-]                       t0 = y1                                               
  <<<<[>>>>+>+<<<<<-]>>>>>[<<<<<+>>>>>-]                        t1 = y
  <<[>>+<[>[-]>+<<-]>>[<<+>>-]<[<<[-]+>>>>>+<<<-]<-<-]>[-]      s1 = y1 greater than y
  >>>>[<<<<+>>>>-]                                              t1 = s1
  <<<<<<<<<<[>>>>>>-<+<<<<<-]>>>>>[<<<<<+>>>>>-]>[>>>>+<<<<[-]] s1 = s1 != s2
  >>>>[-                                                        if (s1)
    <<<<<<[>+>+<<-]>>[<<+>>-]                                   t0 = x1
    <<<[>>>+>+<<<<-]>>>>[<<<<+>>>>-]                            t1 = x2 
    <<[>>+<[>[-]>+<<-]>>[<<+>>-]<[<<[-]+>>>>>+<<<-]<-<-]>[-]    s1 = x1 greater than x2 
    +>>>>[-<<<<-                                                if (s1)
      <<<[>>+<-<-]>>[<<+>>-]                                    x1 = x1 minus x2
    >>>>>]<<<<[-                                                else
      <<<[>>+>+<<<-]>>>[<<<+>>>-]<<[>-<-]>[<+>-]                x1 = x2 minus x1
    >]
    +<<<<<<[->>>>>>-                                            if (s2)
      >>>[<<<<+>+>>>-]<<<[>>>+<<<-]                             t0 = y2
      <<<<[>>>>+<-<<<-]>>>>[<<<<+>>>>-]                         t0 = y2 minus y
    <<<<<<]>>>>>>[-                                             else
      <<<<[>>>+>+<<<<-]>>>>[<<<<+>>>>-]                         t0 = y
      >>>[<<<+<->>>>-]<<<[>>>+<<<-]                             t0 = y minus y2
    ]
    <<[>>>+<<<-]                                                t2 = x1
    >>>[<<[>+<<+>-]>[<+>-]>-]<<[-]                              x1 = x1 * t0
    <<<<[>>>>+>+<<<<<-]>>>>>[<<<<<+>>>>>-]                      t0 = y1
    >>>[<<<+>+>>-]<<[>>+<<-]                                    t1 = y2
    <<[>>+<[>[-]>+<<-]>>[<<+>>-]<[<<[-]+>>>>>+<<<-]<-<-]>[-]    s1 = y1 greater than y2
    +>>>>[-<<<<-                                                if (s1)
      >>>[<<+<<<<<<->>>>>>>>-]<<[>>+<<-]                        y1 = y1 minus y2
    >>>]<<<<[-                                                  else
       >>>[<<<+<+>>>>-]                                         t0 = y2
       <<<<<<<<[>>>>-<<<<-]>>>>[<<<<+>>>>-]>[>>>+<<<-]          y1 = y2 minus y1
    ]
    <<[>+<-]                                                    t0 = x1
    >[<<<<[>>>>>+>+<<<<<<-]>>>>>>[<<<<<<+>>>>>>-]
    <[>+<<-[>>[-]>+<<<-]>>>[<<<+>>>-]<[<-[<<->>[-]]+>-]<-]<<+>] x1 = x1 / y1
    <<[>+>+<<-]>>[<<+>>-]                                       x1 = x1 plus x2
    <<<<<<[>>>>>>+>+<<<<<<<-]>>>>>>>[<<<<<<<+>>>>>>>-]          t0 = x
    <<[>>>+>+<<<<-]>>>[<<<+>>>-]                                t3 = x1  
    <<[>+<-]>>>[<<->+>-]<[>+<-]<[<+>[-]]                        t0 = x != x1
    +<[->-                                                      if (t0)
      <<<<<<<[>>>>>>+>+<<<<<<<-]>>>>>>>[<<<<<<<+>>>>>>>-]       t0 = x
      <[>>+>[<[-]<+>>-]<<[>>+<<-]>[<<[-]+<<<<<<<<[>>>+       
      <<<-]+>>>[<<<->>>-]>>>>>>>-]>-<<<-]>>>[-]<<<              c1 = x greater than x1 then !c2
      <<<<<<[>>>>>>+>+<<<<<<<-]>>>>>>>[<<<<<<<+>>>>>>>-]        t0 = x
      <<[>>>+<<[>>[-]<+<-]>[<+>-]>[<<<[-]+<<<<<<[>>+<<-]+         
      >>[<<->>-]>>>>>>>-]<<-<-]>[-]                             c2 = x less than x1 then !c2
    ]>[-                                                        else
      <<<<<<<<[>>+<<-]+>>[<<->>-]                               C2 = !C2
    >>>>>>]
  >>>>]
  <<<<<<[-]<<<[-]<[-]
  >>>[>+<-],,>>>>>>[<<<<<<<<+>>>>>>>>-],                        Copy (x2 y2) to (x1 y1) read (x2 y2)
]              
<<<<<<<<<<<<[>+<-]>[>>+<<[-]]>>[<<+>>-]                         c2 = c2 || c1
++++++++[<<++++++>>-]<<.                                        Print c

Memory: Pointer movement optimization prioritized.
-------------------------------------------------------------------
C1 | C2 | X | S2 | Y1 | Y | X2 | X1 | T0 | T1 | T2 | T3 | Y2 | S1 |
-------------------------------------------------------------------

Test input:
33,20,41,44,25,12,20

This solution runs into similar problems to all other general solutions due to finite precision mathematics. I've prioritized catching out-of-bounds values over on-the-bounds values. Any value outside of the polygon will return a 0, while any value inside of the polygon will return a 1. Any value on the polygon on a Z-most vertex (top/left/right/bottom) will be indeterminate similar to ray-casting. This algorithm does, however, catch inner-vertices. If you wish to prioritize values on the polygon, you can remove the "if (t0)....else...." section, adjusting memory locations respectively. This will, however, allow (1,4) in the polygon given in the original post to register as true due to finite precision math as described above.

This solution allows a polygon with points in the range of (0-9, 0-9). You can run the program at Le Brainfuck[^].
  Permalink  
v11
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 1

Well ... you could just use standard framework methods:
Point[] points = new Point[] { new Point(2, 0), new Point(4, 1), new Point(4, 4), new Point(2, 5), new Point(1, 2), new Point(2, 0 ) };
GraphicsPath path = new GraphicsPath();
Point last = points[0];
foreach (Point p in points)
    {
    if (p != last)
        {
        path.AddLine(last, p);
        last = p;
        }
    }
Region polygon = new Region(path);
if (polygon.IsVisible(new Point(3, 3)))
    {
    Console.WriteLine("In");
    }
else
    {
    Console.WriteLine("Out");
    }
if (polygon.IsVisible(new Point(5, 4)))
    {
    Console.WriteLine("In");
    }
else
    {
    Console.WriteLine("Out");
    }
  Permalink  
v2
Comments
Chris Maunder 16-Dec-16 10:33am
   
+1 for brevity and understandability. -5 for the overabundance of common sense.
OriginalGriff 16-Dec-16 10:43am
   
I knew that would be my undoing! :laugh:
   
Maybe Marc can do this in COBOL.
Ramza360 16-Dec-16 11:06am
   
This works well however you do not need the foreach, use the GraphicsPath.AddLines(points) does the same thing. Quick addition to your solution i added below.

And yes it is simple in c#
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 4

If I recall you take an imaginary line from your point (x,y) to (xMax,y) and check if each of the vertices crosses it.

If the number that crosses is even you are outside the polygon (inclusive of 0) otherwise you are inside the polygon.
  Permalink  
Comments
ppolymorphe 20-Dec-16 16:46pm
   
Is it considered a solution ?
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 6

If you allow me an extension method so I can get fluent coding:

public static class GPExt
{
    public static GraphicsPath AddPolygon2(this GraphicsPath gp, PointF[] points)
    {
        gp.AddPolygon(points);
        return gp;
    }
}


I can write this as a one line method, parsing an input string of coordinates (I hate all those new Point calls, haha):

static bool Contains(string points, PointF p)
{
    return new GraphicsPath()
    .AddPolygon2(
        new Regex(@"(?<coord>(\((?<x>([-+]?[\d]+(\.\d+)*)),\s*(?<y>([-+]?[\d]+(\.\d+)*))\),*\s*))")
            .Matches(points)
            .Cast<Match>().Select(m => new PointF(Convert.ToSingle(m.Groups["x"].Value), Convert.ToSingle(m.Groups["y"].Value))).ToArray())
    .IsVisible(p);
}


And you call it like this:

Console.WriteLine(Contains("(2,0),(4,1),(4,4),(2,5),(1,2)", new PointF(3, 3)));
Console.WriteLine(Contains("(2,0),(4,1),(4,4),(2,5),(1,2)", new PointF(5, 4)));


Results being:

True
False

Rather a wild combination using regex, Linq, and a user friendly input of the polygon.

Marc
  Permalink  
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 7

Picked out this new mission while waiting for my flight home. One look and this is going to involve geometry, that means Cartesian coordinate system, vectors, edges, and angles sort of things. To kill time, I started drawing out my thought on the phone...

There are three ways to solve this problem:
1. The easiest one is to pluck and use one of the many Python packages for manipulation and analysis of geometric objects on the Cartesian plane.
2. The challenging one will be to find your own algorithm (excluding known ones) to determine that a point is inside or outside a polygon.
3. The last and extreme one will be to explore certain machine learning techniques such as neural network to learn and classify a point to be either inside or outside of a polygon. This is time-consuming, as it will need a big data set of polygons and points with known results for the computer to learn. I will not be attempting it here.

Let's begin...

The Easiest One

I am using the Shapely 1.6b2[^] for this demo. The code below is commented with URLs that point to the respective online references:
"""
IsPointInPolygon_1.py
The easier way using shapely package
by Peter Leow the point seeker
 
"""
 
from shapely.geometry.polygon import Polygon
from shapely.geometry.point import Point
 
def pointWithinPolygon(point, polygon):
    # http://toblerity.org/shapely/manual.html#object.within
    return point.within(polygon)
 
def polygonContainsPoint(point, polygon):
    # http://toblerity.org/shapely/manual.html#object.contains
    return polygon.contains(point)
 
def main():
 
    # http://toblerity.org/shapely/manual.html#polygons
    polygon = Polygon([(2,0), (4,1), (4,4), (2,5), (1,2)])
    
    while True:
        # http://toblerity.org/shapely/manual.html#points
        coords = input("Input the x y coordinates of a point separated by space: ").split()
        if len(coords) != 2:
            break
 
        # Numeric validation omitted
        point= Point(float(coords[0]), float(coords[1]))
 
        resultWithin = 'Is {} within {}? {}.'.format(point, polygon, pointWithinPolygon(point, polygon))
 
        resultContains = 'Does {} contain {}? {}.'.format(polygon, point, polygonContainsPoint(point, polygon))
 
        print(resultWithin)
 
        print(resultContains)
  
main()

An example output:
Input the x y coordinates of a point separated by space: 3 3
Is POINT (3 3) within POLYGON ((2 0, 4 1, 4 4, 2 5, 1 2, 2 0))? True.
Does POLYGON ((2 0, 4 1, 4 4, 2 5, 1 2, 2 0)) contain POINT (3 3)? True.
Input the x y coordinates of a point separated by space: 5 4
Is POINT (5 4) within POLYGON ((2 0, 4 1, 4 4, 2 5, 1 2, 2 0))? False.
Does POLYGON ((2 0, 4 1, 4 4, 2 5, 1 2, 2 0)) contain POINT (5 4)? False.
Input the x y coordinates of a point separated by space: 4 4
Is POINT (4 4) within POLYGON ((2 0, 4 1, 4 4, 2 5, 1 2, 2 0))? False.
Does POLYGON ((2 0, 4 1, 4 4, 2 5, 1 2, 2 0)) contain POINT (4 4)? False.
Input the x y coordinates of a point separated by space: 

You will have to install the shapely package from Shapely 1.6b2 : Python Package Index[^] prior to running this program.

The Challenging One

As I was exploring the different polygons and points on the phone, I observed one trait that seems to prevail only when a point is inside a polygon--The sum of all angles from the point to adjacent vertices of the encompassing polygon is always 2pi in radians (or 360 degrees). I wrote the following code to implement this assumption:
"""
IsPointInPolygon_2.py
The pain-in-the-a** way
by Peter Leow the point seeker
 
"""
 
import math
 
def offsetCoords(polygon, newOrigin):
    offsetPolygon = []
    for v in polygon:
        offsetPolygon.append((v[0]-newOrigin[0], v[1]-newOrigin[1]))
    return offsetPolygon
 
def dotProduct(v1, v2):
    return v1[0]*v2[0]+v1[1]*v2[1]
 
def vectorLen(v):
    return (v[0]**2+v[1]**2)**.5
 
def angleBetweenVectors(v1, v2):
    cosine = dotProduct(v1, v2) / (vectorLen(v1) * vectorLen(v2))
    return math.acos(cosine)
 
def isWithin(point, polygon):
    # weed out the obvious
    if point in polygon:
        return False
    # Find angle of each adjacent vectors
    sumAngles = 0
    for i in range(len(polygon) - 1):
        sumAngles += angleBetweenVectors(polygon[i], polygon[i+1])
    if math.isclose(sumAngles, math.pi*2):
        return True
    else:
        return False
 
def main():
 
    # A list of coords tuples
    polygon = [(2,0), (4,1), (4,4), (2,5), (1,2), (2,0)]
    
    while True:
        coords = input("Input the x y coordinates of a point separated by space: ").split()
        if len(coords) != 2:    # Exit on Enter key
            break
 
        # Numeric validation omitted
        point = (float(coords[0]), float(coords[1]))
 
        '''
        Set this point as the origin (0 0) of the cartesian coordinate plane
        by offsetting all vertices of the polygon against it
        '''
        offsetPolygon = offsetCoords(polygon, point)
        offsetPoint = (0, 0)
 
        result = 'Is POINT {} within POLYGON {}? {}.'.format(point, polygon, isWithin(offsetPoint, offsetPolygon))
 
        print(result)
  
main()

An example output is shown below:
Input the x y coordinates of a point separated by space: 3 3
Is POINT (3.0, 3.0) within POLYGON [(2, 0), (4, 1), (4, 4), (2, 5), (1, 2), (2, 0)]? True.
Input the x y coordinates of a point separated by space: 5 4
Is POINT (5.0, 4.0) within POLYGON [(2, 0), (4, 1), (4, 4), (2, 5), (1, 2), (2, 0)]? False.
Input the x y coordinates of a point separated by space: 4 4
Is POINT (4.0, 4.0) within POLYGON [(2, 0), (4, 1), (4, 4), (2, 5), (1, 2), (2, 0)]? False.
Input the x y coordinates of a point separated by space: 

So far so good, my assumption seems to hold true. Anyone have time to spare? Please help to test it out.
  Permalink  
v9
Comments
Jon McKee 18-Dec-16 15:05pm
   
Just wanted to point out that the 2pi solution is called the "Winding Number Algorithm" if you were curious =)
Peter Leow 18-Dec-16 21:15pm
   
Thank you for enlightening me. Will read up on it, and see if I can come out with a better code with time permitting. This code challenge has become the hotbed for knowledge sharing. Cheer!
grgran 19-Dec-16 10:51am
   
It would seem to me that particularly irregular polygons could produce sums greater than 2pi, because "wedges" (triangles formed by the test point and edge points) could overlap. Imagine a shape where a "left" edge point connected to a point into the polygon inside and near the "right" edge. I wish I could draw here :-)
Peter Leow 19-Dec-16 21:56pm
   
You are right. The overlapping results in sum of angle > 2pi. So my assumption is false. Thank you.
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 9

Here is my submission.
Looks like I took a little original approach. I guess the method is known but I have no clue about it. The method works with clockwise and non clockwise polygons, convex or not.
First I move the polygon as if the testing point was {0,0}.
Then I compute the angle covered by a segment using arc tangent function and sum it for every segments.
If sum is 0, point is outside of polygon or on border.
If sum is 2Pi (radians), point is inside of polygon.

*   CCCP Code Challenge Code Project
    * Point in polygon
    clear
    ? "CCCP Code Challenge Code Project"
    ? "Point in polygon"
    /*
         2
     x  x
      1
    x    x
      x
    
    */
    poly= {{2,0}, {4,1}, {4,4}, {2,5}, {1,2}, {2,0}}
    points= {{3,3}, {5,4}, {4,4}}
    FOR scan= 1 to len(points)
        rep= inpoly(points[scan], poly)
        ? rep
        ?
    next
    return
    
    function inpoly(point, poly)
        local scan
        somme= 0
        for scan=1 to len(poly)- 1
            // move segment coordinates as if testing point was {0,0}
            x1= poly[scan,1]- point[1]
            y1= poly[scan,2]- point[2]
            x2= poly[scan+1,1]- point[1]
            y2= poly[scan+1,2]- point[2]
            // calc covering angle of segment
            tmp= atn2(x2,y2)- atn2(x1,y1)
            if tmp > pi()
                tmp -= 2*pi()
            endif
            if tmp < -pi()
                tmp += 2*pi()
            endif
            // sum it
            somme += tmp
        next
        return abs(somme) > 3
  Permalink  
v7
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 10

Abusing Chris' generosity regarding data type, I made a solution which has rather nice runtime complexity :p
#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;
 
class point
{
public:
	bool x, y;
	point(bool _x, bool _y) : x(_x), y(_y) {}
	bool operator==(const point& o) const { return o.x == x && o.y == y; }
};
 
typedef vector<point> polygon;
 
bool is_point_in_polygon(const point& pt, const polygon& area)
{
	return find(area.begin(), area.end(), pt) != area.end();
}
 
int main()
{
	polygon area;
	area.push_back(point(true, true));
	area.push_back(point(true, false));
	area.push_back(point(false, true));
 
	point pt1(false, false);
	cout << (is_point_in_polygon(pt1, area) ? "Inside :)" : "Outside :(") << endl;
 
	point pt2(true, false);
	cout << (is_point_in_polygon(pt2, area) ? "Inside :)" : "Outside :(") << endl;
 
    return 0;
}
giving the output
Outside :(
Inside :)
The boolean 2-D space is a dreamworld when working with geometric algorithms since it only contains 4 points. A coordinate cannot lie truly within a polygon, but must be on one vertex or it's outside.
Although the boolean 2-D space is of low practical use, the pure theoretical concept is almost as bad, if not worse.
Next week I'm hoping for something more challenging, like computing the convex hull or doing polygon union. :)
  Permalink  
v2
Rate this: bad
 
good
Please Sign up or sign in to vote.

Solution 12

using System;
using System.Collections.Generic;
using System.Windows.Media.Media3D;
 
class Program
{
    public static bool TestPointAgainstPolygon(List<Vector3D> Polygon,
                                               Vector3D TestPoint)
    {
        Vector3D Normal;
        int i;
 

        if (Polygon == null && Polygon.Count < 3)
        {
           return false;
        }
 
        for (i = 0; i < (Polygon.Count - 1); i++)
        {
            Normal = Vector3D.CrossProduct(Polygon[i + 1] - Polygon[i],
                                           TestPoint - Polygon[i]);
            if (Normal.Z < 0)
            {
                return false;
            }
        }
 
        return true;
    }
 
    static void Main(string[] args)
    {
        bool Result;
 

        List<Vector3D> Polygon = new List<Vector3D> 
                                           { 
                                                new Vector3D(2,0,0),
                                                new Vector3D(4,1,0),
                                                new Vector3D(4,4,0),
                                                new Vector3D(2,5,0),
                                                new Vector3D(1,2,0),
                                                new Vector3D(2,0,0)
                                           };
 
        Result = TestPointAgainstPolygon(Polygon, new Vector3D(3, 3, 0));
        if (Result == true)
        {
            Console.WriteLine("Point (3,3) is inside the polygon.");
        }
        else
        {
            Console.WriteLine("Point (3,3) is outside the polygon.");
        }
 
        Result = TestPointAgainstPolygon(Polygon, new Vector3D(5, 4, 0));
        if (Result == true)
        {
            Console.WriteLine("Point (5,4) is inside the polygon.");
        }
        else
        {
           Console.WriteLine("Point (5,4) is outside the polygon.");
        }
 
        Console.Read();
    }
}


Output:
Point (3,3) is inside the polygon.
Point (5,4) is outside the polygon.


This approach is a little different. The vectors are expanded to 3D, with the Z components set to zero. This enables us to use the cross product to calculate the normal between the to the point and each line of the polygon.

If the polygon is defined counterclockwise, then each normal's Z component must be positive if the point is inside the polygon. It is outside if even a single normal's Z component is negative.

Note: We don't actually need to calculate the full cross product. It is sufficient to determine the sign of the Z component. Also, we can drop the requirement that the polygon is defined counterclockwise. If all normals have the same sign (no matter if positive or negative), then the point is inside the polygon. If only one sign is different, then the point is outside.

Note 2: This was a question in one of my exams. The professor also gave us the polygon in 2D and expected us to expand it to 3D and use the cross product.
  Permalink  
v3
Comments
ppolymorphe 20-Dec-16 16:52pm
   
Does it work if the polygon is not convex ?
CDP1802 20-Dec-16 16:58pm
   
That's the weak point. Fortunately you can divide every concave polygon up into convex polygons. The other nasty case would be if the polygon could contain another polygon as a 'hole'.
ppolymorphe 20-Dec-16 18:15pm
   
I know Delaunay triangulation. I fear no solution here can handle holes in polygon.
Have a look at:
href="http://www.cs.cmu.edu/~quake/robust.html">Fast Robust Predicates for Computational Geometry[^]
In this page, there a 2D cross product.
CDP1802 20-Dec-16 18:24pm
   
The orientation test looks familiar.
ppolymorphe 20-Dec-16 18:33pm
   
I found the page while studying Voronoi/Delaunay.
The interest is to avoid resorting to 3D points.

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

  Print Answers RSS
Top Experts
Last 24hrsThis month


Advertise | Privacy |
Web01 | 2.8.170915.1 | Last Updated 6 Jan 2017
Copyright © CodeProject, 1999-2017
All Rights Reserved. Terms of Service
Layout: fixed | fluid

CodeProject, 503-250 Ferrand Drive Toronto Ontario, M3C 3G8 Canada +1 416-849-8900 x 100