TheSwamp

Code Red => .NET => Topic started by: svenoaks on May 13, 2013, 06:51:20 PM

Title: Algorithm for points inside a polyline inexplicably fails
Post by: svenoaks on May 13, 2013, 06:51:20 PM
Crossposted at http://stackoverflow.com/questions/16531963/algorithm-for-points-inside-a-polyline-inexplicably-fails-autocad (http://stackoverflow.com/questions/16531963/algorithm-for-points-inside-a-polyline-inexplicably-fails-autocad) No answers yet.

I have an algorithm (which I did not write) for determining if a point lies within a polygon (straight lines only).

I have been testing my command on the same 51 polygons repeatedly. 99% it will work perfectly. Every once in a while it will fail on 1 or more of the polygons, returning false for over 2000 points I am creating inside the bounding box of the polyline. I have seen it fail when the polyline is a simple rectangle and all of the points lie distributed in a grid within the polyline. It should have returned true over 2000 times in that case. It will never fail for just 1 of the points, it will fail all of them. I have confirmed that the points are being correctly created where I expect them to be and that the vertices of the polygon are where I expect them to be. When it fails, the angle variable for the last point is at exactly double PI.

I am not doing any multi-threading. The only possibly 'funny' thing I am doing is COM Interop with Excel. This is happening after the transaction has been committed for the part with this algorithm, and I am sure I am cleaning up all my COM objects. I have not been able to reproduce the failure without the COM Interop part but I don't think I've tested it enough yet to have enough absence of evidence.

Any ideas what could be wrong?

Code - C#: [Select]
  1.     bool IsInsidePolygon(Polyline polygon, Point3d pt)
  2.     {
  3.         int n = polygon.NumberOfVertices;
  4.         double angle = 0;
  5.         Point pt1, pt2;
  6.  
  7.         for (int i = 0; i < n; i++)
  8.         {
  9.             pt1.X = polygon.GetPoint2dAt(i).X - pt.X;
  10.             pt1.Y = polygon.GetPoint2dAt(i).Y - pt.Y;
  11.             pt2.X = polygon.GetPoint2dAt((i + 1) % n).X - pt.X;
  12.             pt2.Y = polygon.GetPoint2dAt((i + 1) % n).Y - pt.Y;
  13.             angle += Angle2D(pt1.X, pt1.Y, pt2.X, pt2.Y);
  14.         }
  15.  
  16.         if (Math.Abs(angle) < Math.PI)
  17.             return false;
  18.         else
  19.             return true;
  20.     }
  21.  
  22.     public struct Point
  23.     {
  24.         public double X, Y;
  25.     };
  26.  
  27.     public static double Angle2D(double x1, double y1, double x2, double y2)
  28.     {
  29.         double dtheta, theta1, theta2;
  30.  
  31.         theta1 = Math.Atan2(y1, x1);
  32.         theta2 = Math.Atan2(y2, x2);
  33.         dtheta = theta2 - theta1;
  34.         while (dtheta > Math.PI)
  35.             dtheta -= (Math.PI * 2);
  36.         while (dtheta < -Math.PI)
  37.             dtheta += (Math.PI * 2);
  38.         return (dtheta);
  39.     }
Title: Re: Algorithm for points inside a polyline inexplicably fails
Post by: TheMaster on May 13, 2013, 11:33:41 PM
Going on the assumption that your algorithm is correct, it could be a result of roundoff errors, depending on the location of the coordinates and the vertices (e.g., extremely large coordinates).

You could also try using Regions and the Brep API to do this, which works for any polyline that can be used to create a simple region:

   http://www.theswamp.org/index.php?topic=43572.msg488632#msg488632



Title: Re: Algorithm for points inside a polyline inexplicably fails
Post by: svenoaks on May 14, 2013, 07:25:19 AM
Don't think it's round off errors, my numbers are 'normal' sized. I'll definitely give the Brep API a whirl, looks very promising - thanks.
Title: Re: Algorithm for points inside a polyline inexplicably fails
Post by: WILL HATCH on May 14, 2013, 11:51:57 AM
The algorithm is flawed.  See the breakout:
In this case it is the result of bounding dtheta, the removal of the two while statements would fix it.  This is because by bounding that difference we are getting the smaller angle between the two vectors (which in this case are confusingly labeled pt1 and pt2) and because you are dealing with a simple convex polygon the angle will always be in the range [0,pi] or [-pi,0] depending on rotation direction of points, and the sum will always be either -2pi or 2pi
Title: Re: Algorithm for points inside a polyline inexplicably fails
Post by: svenoaks on May 14, 2013, 05:30:50 PM
Wow, thanks Will - I will try to sit down and understand the next breather I get.