Point in polygon OR point on polygon using LINQ

As noted in an earlier question, How to Zip enumerable with itself, I am working on some math algorithms based on lists of points. I am currently working on point in polygon. I have the code for how to do that and have found several good references here on SO, such as this link Hit test. So, I can figure out whether or not a point is in a polygon. As part of determining that, I want to determine if the point is actually on the polygon. This I can also do. If I can do all of that, what is my question you might ask?

Can I do it efficiently using LINQ? I can already do something like the following (assuming a Pairwise extension method as described in my earlier question as well as in links to which my question/answers links, and assuming a Position type that has X and Y members). I have not tested much, so the lambda might not be 100% correct. Also, it does not take very small differences into account.

EDIT Note that most recent edit of this code has changed from TakeWhile to TakeWhileInclusive. See the extension method towards the end of this post.

//
// Extend a ray from pt towards the left.  Does this ray intersect the segment p1->p2?
// By definition, the ray extending from pt cannot intersect a horizontal segment, so our
// first check is to see whether or not the segment is horizontal.  
// If it is, there cannot be an intersection.
// If it is not, there could be an intersection.
//
public static bool PointIntersectSegment(Position p1, Position p2, Position pt)
{
  bool intersect = false;
  if (p1.Y != p2.Y)
  {
    // Is pt between (vertically) p1 and p2?  
    // If so, the ray from pt might intersect.  
    // If not, the ray from pt cannot intersect.
    if ((p1.Y >= pt.Y && p2.Y < pt.Y) || (p1.Y < pt.Y && p2.Y >= pt.Y))
    {
      if (p1.X < pt.X && p2.X < pt.X)
      {
        // If the segment is to the left of pt, then the ray extending leftwards from pt will intersect it.
        intersect = true;
      }
      else
        if ((p1.X < pt.X || p2.X < pt.X))
        {
          // If either end of the segment is to the left of pt, calculate intersection (x only) of the 
          // ray from pt and the segment.  If the intersection (x only) is to the left of pt, then
          // the ray extending leftwards from pt will intersect it.
          double inverseSlope = (p1.X - p2.X) / (p1.Y - p2.Y);
          double intersectionX = (pt.Y - p1.Y) * inverseSlope + p1.X;
          if (intersectionX < pt.X)
          {
            intersect = true;
          }
        }
    }
  }
  return intersect;
}

public static bool PointOnSegment(Position p1, Position p2, Position pt)
{
  // Obviously, this is not really going to tell us if pt is on the segment p1->p2.  I am still
  // working on that.  For testing the PointInPolygon algorithm, I can still simulate the "on"
  // case by passing in a pt that is equal to one of the points on the polygon.
  return (pt == p1 || pt == p2);
}

public static PointInPolygonLocation PointInPolygon(IEnumerable<Position> pts, GTPosition pt)
{
  //
  // Implemention of the Jordan Curve theorem to determine if a point is in a polygon.  In essence,
  // this algorithm extends a ray infinitely from pt.  In my implementation I am extending to the left.
  // If the point is inside the polygon, then the ray will intersect the polygon a odd number of times.
  // If the point is outside the polygon, then the ray will intersect the polygon an even number of times.
  // Ideally, we would be able to report not only if the point is inside or outside of the polygon, but
  // also if it is "on" the polygon (i.e. it lies on one of the segments).  Note that "on" and 
  // inside/outside are exlusive.  The point is either on the polygon or it is inside or outside, it
  // cannot be "inside and on" or "outside and on".
  // So, the algorithm is as follows:
  // 1.  Consider the points of the polygon as pairs making up the segments (p1->p2, p2->p3, etc).
  // 2.  For each segment, perform two calculations:  
  //       Does the ray extending left from pt intersect the segment?
  //       Is pt on the segment p[i], p[i+1]?
  // 3.  Count the total number of intersections.  If odd, point is inside.  If even, point is outside.
  // 4.  Here is the tricky part, if the point is on any segment, then we can stop.  If it is on the 
  //     first segment, then there is no need to go through the on/off calculation and the intersection
  //     calculation for the rest of the segments.
  //
  var result = pts.Pairwise( (p1, p2) =>
                         {
                           int intersect = 0;
                           int on = 0;

                           on = PointOnSegment(p1, p2, pt) ? 1 : 0;
                           if (on == 0)
                           {
                             //Don't really need to determine intersection if we already know that pt is on p1->p2.
                             intersect = PointIntersectSegment(p1, p2, pt) ? 1 : 0;
                           }
                           return new { on, intersect };
                         })
                         .TakeWhileInclusive((item) => item.on == 0) //Only consider segments until (or if) pt is on a segment.
                         .Aggregate((curr, next) => new     //Keep a running total of how many intersections.
                                                    { 
                                                      on = curr.on + next.on, 
                                                      intersect = curr.intersect + next.intersect 
                                                    });
  if (result.on != 0)
  {
    return PointInPolygonLocation.On;
  }
  else
  {
    return result.intersect % 2 == 0 ? PointInPolygonLocation.Outside : PointInPolygonLocation.Inside;
  }
}

This function, PointInPolygon, takes the input Position, pt, iterates over the input sequence of position values, and uses the Jordan Curve method to determine how many times a ray extended from pt to the left intersects the polygon. The lambda expression will yield, into the "zipped" list, 1 for every segment that is crossed, and 0 for the rest. The sum of these values determines if pt is inside or outside of the polygon (odd == inside, even == outside). So far, so good.

Now, for any consecutive pairs of position values in the sequence (i.e. in any execution of the lambda), we can also determine if pt is ON the segment p1, p2. If that is the case, we can stop the calculation because we have our answer.

Ultimately, my question is this: Can I perform this calculation (maybe using Aggregate?) such that we will only iterate over the sequence no more than 1 time AND can we stop the iteration if we encounter a segment that pt is ON? In other words, if pt is ON the very first segment, there is no need to examine the rest of the segments because we have the answer.

It might very well be that this operation (particularly the requirement/desire to possibly stop the iteration early) does not really lend itself well to the LINQ approach.

It just occurred to me that maybe the lambda expression could yield a tuple, the intersection value (1 or 0 or maybe true or false) and the "on" value (true or false). Maybe then I could use TakeWhile(anontype.PointOnPolygon == false). If I Sum the tuples and if ON == 1, then the point is ON the polygon. Otherwise, the oddness or evenness of the sum of the other part of the tuple tells if the point is inside or outside.

EDIT Cleaned up the code formatting, moved the lambda to a standalone function, added a new function to if the point is on a segment. The new function is really a dummy, it just compares the input point to the start and end points of the segment and returns true if either is equal. So, it should probably be called PointIsOnEndpointOfSegment. There is already a lot of code here and I didn't want to cloud the issue with more fooli开发者_高级运维ng around with x's and y's to do the "real" ON calculation.

With the code as it is now, I can do the following: 1. View sequence of points as a sequence of pairs (segments) (via Pairwise). 2. In the Pairwise result selector, I compute an anonymous type containing the value intersection=1 if the ray extending leftwards from pt intersects the segment and the value on=1 if the point is "on" the segment. 3. Express the intersection/on anonymous types as a sequence.

I thought I could use TakeWhile to take all of the anonymous types until (or if) one of the anonymous types indicates that the pt was ON a segment. Then, using Aggregate, sum the values of intersection and on. If the sum of on != 0, then the point was on one of the segments. Otherwise, the sum of intersection will tell us if the point is inside (odd) or outside (even - or zero).

Unfortunately, I ran into the a similar problem with TakeWhile as was described here:

How to TakeWhile PLUS one more element.

My TakeWhile takes all of the items from the sequence up to, but not including, the item (if any) that indicates an intersection. So, when I aggregate the results, the intersecting item is not there.

It looks like one way that people have handled this situation before has been to write a TakeUntil extension that is like TakeWhile, but it includes the first item that fails the predicate. Would such a TakeUntil extension force an evaluation of the entire sequence?

This has mainly been an exercise to see if this algorithm could be implemented using Linq subject to the following requirements: 1. Entire sequence of input points is iterated ONLY if the input point is NOT on any of the segments. 2. If the input point is on one of the segments, the iteration over the input points should stop. If the input point is on the first segment of a 1000 point polygon, there is no need to do the intersection and "on" calculations for the rest the segments.

We have implemented this algorithm previously in C++ using a for loop over the points, breaking out if the point is ever on a segment. I could certainly implement the same way, I just wanted to see if I could LINQify if without incurring more interations than the "old" for loop way.

I will probably try the TakeUntil approach just for kicks and see what happens.

Edit

With this code (TakeWhileInclusive IEnumerable that is virtually identical to the TakeUntil described in the link about TakeWhile PLUS one more element):

// Mimic TakeWhile's overload which takes an index as a parameter to the predicate.
public static IEnumerable<T> TakeWhileInclusive<T>(this IEnumerable<T> seq, Func<T, int, bool> predicate)
{
  int i = 0;
  foreach( T e in seq)
  {
    if (!predicate(e,i))
    {
      // If here, then this is first item to fail predicate.  Yield this item and then break.
      yield return e;
      yield break;
    }
    // yield each item from the input sequence until end of sequence or first failure (see above).
    yield return e;
    i++;
  }
}

//
// First saw this here: http://blog.foxxtrot.net/2009/06/inclusively-take-elements-using-linq-and-custom-extensions.html
// TakeWhileInclusive - IEnumerable extension to mimic TakeWhile, but also to return the first
// item that failed the predicate.
// e.g.  seq = 1 2 3 4
// TakeWhile(p => p != 3) will yield 1 2
// TakeWhileInclusive(p => p != 3) will yield 1 2 3
// e.g.  seq = 0 0 0 1 0
// TakeWhile(p => p == 0) will yield 0 0 0
// TakeWhileInclusive(p => p == 0) will yield 0 0 0 1
// Similar to TakeUntil from here: https://stackoverflow.com/questions/2242318/how-could-i-take-1-more-item-from-linqs-takewhile
//
public static IEnumerable<T> TakeWhileInclusive<T>(this IEnumerable<T> seq, Func<T, bool> predicate)
{
  return seq.Select((x, i) => new { Item = x, Index = i })
            .TakeWhileInclusive((x, i) => predicate(x.Item))
            .Select(x => x.Item);
}

And the corresponding change in my original code (replace TakeWhile with TakeWhileInclusive), I can get the answer (either the point is on the polygon or it is inside or outside) and the iteration stops (via TakeWhileInclusive) after it hits the segment of the polygon that the point is ON. I will note again that my "PointOnSegment" code is bogus, but is fine for testing as long as I am testing a point that is equal to one of the polygon points.

I will leave this for now in case anyone else wants to comment. Right now I am inclined to accept my own answer as it does what I intended to do. Whether or not this particular approach is ultimately a good one, I don't know yet.

I don't quite understand the algorithm itself, but I think I can help:
What you basically need, assuming I understood that much, is to check if any of the tuples fit a certain predicate. Then you can put each pair of Positions in a tuple, and pass it along as an IEnumerable<Tuple<Position, Position>> on which you can do .Any(predicate goes here) to see if any of them fit the description (or, if you need the actual value for which the predicate was true, use .FirstOrDefault() and check null)
So your code will look something like this:

   public static PointInPolygonLocation PointInPolygon(IEnumerable<Position> pts, Position pt)
    {
        bool isIn = pts.Pairwise((p1, p2) => Tuple.Create(p1, p2)).Any(tuple => tuple.Item1.Y != tuple.Item2.Y &&
        ((tuple.Item1.Y >= pt.Y && tuple.Item2.Y < pt.Y) || (tuple.Item1.Y < pt.Y && tuple.Item2.Y >= pt.Y))
        && ((tuple.Item1.X < tuple.Item1.X && tuple.Item2.X < pt.X) || ((pt.Y - tuple.Item1.Y) * ((tuple.Item1.X - tuple.Item2.X) / (tuple.Item1.Y - tuple.Item2.Y)) * tuple.Item1.X) < pt.X));
        return isIn ? PointInPolygonLocation.Inside : PointInPolygonLocation.Outside;
    }

I would seriously consider adding some comments to that monstrous predicate. I kind of zipped it up. nested ifs turned into &&s, etc. You should check it, ofcourse, but I'm 90% sure it's the same as what you did. To quote knuth : Beware of bugs in the above code; I have only proved it correct, not tried it.

Edit: Much more readable now, wouldn't you say?

   public static PointInPolygonLocation PointInPolygon(IEnumerable<Position> pts, Position pt)
    {
        bool isIn = pts.Pairwise((p1, p2) => Tuple.Create(p1, p2)).Any(tuple => ReadablePredicate(tuple.Item1, tuple.Item2, pt));
        return isIn ? PointInPolygonLocation.Inside : PointInPolygonLocation.Outside;
    }
    public static bool ReadablePredicate(Position p1, Position p2, Position pt)
    {
        if (p1.Y == p2.Y)
            return false;
        if (!((p1.Y >= pt.Y && p2.Y < pt.Y) || (p1.Y < pt.Y && p2.Y >= pt.Y)))
            return false;
        if (p1.X < pt.X && p2.X < pt.X) // Originally was (p1.X < p1.X && ...). that makes no sense, so I assumed you meant p1.X < pt.X
            return true;
        if ((p1.X < pt.X || p2.X < pt.X) && ((pt.Y - p1.Y) * ((p1.X - p2.X) / (p1.Y - p2.Y)) * p1.X) < pt.X)
            return true;
        return false;
    }

I have been through a few iterations (ha ha) of this problem. Ultimately there were two keys:

1. My result selector function computes two values per segment of the polygon. Does the ray extending leftwards from the point intersect the segment? (1 or 0 to facilitate summing via Aggregate later on). Is the point ON the segment? (1 or 0 to facilitate summing via Aggregate later on). If the point on segment value is 1 (point is on segment), then the iteration should stop.

2. Since I want to iterate over as few items as necessary, I need to limit the sequence in the case of the point being on a segment. If the point is not on any segment, I have not choice but to iterate over all segments. If the point is on a segment, I want to quit, but I also need to Aggregate that one last value (that caused the "iterate until the point is on a segment" condition to fail). I though that TakeWhile would work as I would take the { intersect, on } pairs until on == 1 (i.e. until there the point is on a segment). Unfortunately, TakeWhile detects this condition and stops the iteration, but my output (which feeds to Aggregate) does not contain the "failure" condition, so the results of Aggregate to did not reflect that the point was on a segment.

I thought that, if I were to accomplish this via LINQ, I would have to iterate over the entire sequence and then check the results (of Aggregate) to determine if the point had landed on any of the segments over which I just iterated. Using the TakeWhileInclusive extension method inspired/copied from the post linked above and from a non-SO website noted in the comment in the code, I am able to get the TakeWhile + 1 behavior that I needed.

Maybe there is a better way to do this? I don't know. I am really just getting started. The Pairwise pattern works well for many geometry calculations that involve sequences of points. Given a sequence of points and an algorithm expressed in terms of pairs of the points (p1->p2, p2->p3, p3->p4, etc), it can be "simply" a matter of expressing the algorithm in a lambda (easy ones like Distance) or in a standalone function (like the one in this post, Centroid, etc).

I will leave this for now in case anyone else want to comment or make suggestions or come up with alternate answers. Ultimately, I will probably accept this answer (as expressed in the additional code and comments in the original question.