Author Topic: Point in poly  (Read 4109 times)

0 Members and 1 Guest are viewing this topic.

Bryco

  • Water Moccasin
  • Posts: 1883
Point in poly
« on: July 17, 2007, 12:47:10 AM »
I've been working on converting the vba version I had to C#, I'm getting there.
The winding method used by http://softsurfer.com is far superior to the usual counting ray crossing method, and does work for selfintersecting polylines. Unfortunately it doesn't work for curves.
Although I have added some math that works fine for bumps, it's pretty expensive and I don't understand math well enough to make it leaner. Perhaps someone can help. I hope the actual code is not too offensive to those who know C#. I gotta say that I totally miss the point with enums.
Code: [Select]
using Autodesk.AutoCAD.ApplicationServices;
using Autodesk.AutoCAD.DatabaseServices;
using Autodesk.AutoCAD.EditorInput;
using Autodesk.AutoCAD.Geometry;
using Autodesk.AutoCAD.Runtime;
using System.Collections;
using System;
using cad = Autodesk.AutoCAD.ApplicationServices.Application;
using Util;

namespace Getit
{

    public class PtinPoly
    {
       
        ////Uses http://softsurfer.com but also makes it work for curves

        //07-16-07 CCP
       [Flags]
       public enum isInside
        {Inside,Outside,Selfintersecting}

        public static isInside PtInsidePoly(Point3d p, Polyline oPline)
        {
            //1) Use the winding count to see if the point is inside a straigth sided poly
            //2) Then checks for arcs
            Editor ed = cad.DocumentManager.MdiActiveDocument.Editor;           
            Point2d P1, P2;
            double dBulge;
            int Testbulge, wn = 0; //winding count
            int Cnt = oPline.NumberOfVertices;
            Point2dCollection PtS=new Point2dCollection();

            for (int i = 0; i < Cnt; i++){
                Point2d Pt = oPline.GetPoint2dAt(i);
                PtS.Add(Pt);
            }
            PtS.Add(PtS[0]);
            PtS.TrimToSize();

            for (int i = 0; i < Cnt; i++){ //loop through all edges of the polygon
                P1 = PtS[i];
                P2 = PtS[i + 1];
                if (P1.Y <= p.Y){      //start y <= P.y
                    if (P2.Y > p.Y){   //an upward crossing
                        if (MyMath.isLeft(P1, P2, p) > 0){  //P left of edge Cw 
                            ++wn; }   //have a valid up intersect
                    }
                }
                else{                 //start y > P.y (no test needed)
                    if (P2.Y <= p.Y){      //a downward crossing
                        if (MyMath.isLeft(P1, P2, p) < 0){   //P right of edge Cww
                            --wn; //have a valid down intersect
                        }                           
                    }
                }
            }

            if ((oPline.IsOnlyLines)){    //Lines only   
                if (wn == 0){
                    return isInside.Outside;
                }
                else{
                    ed.WriteMessage("\nWinding={0}", wn);
                    return isInside.Inside;
                }
            }

            //Now for the arcs
            //if( wn = 0 ) //Point is inside verticies using straight lines
            //if( bumps go out and the point is in one ) point is in the poly
            //The opposite applies if( wn<>0)

            if (SelfIntersectingCheck(oPline)){  //This is hardly necessary as it often gives the correct answer
                ed.WriteMessage("\nThis poly is self intersecting.");
                return isInside.Selfintersecting;
            }

            bool Cw = MyMath.IsClockwise(oPline);
            for (int i = 0; i < Cnt; i++){
                dBulge = oPline.GetBulgeAt(i);
                if (dBulge == 0) goto skip;  // continue             
                //Now we only have curves
                Testbulge = Math.Sign(dBulge);
                if (wn == 0) { Testbulge = -Testbulge; }

                if (Cw){    //poly is clockwise
                    if (Testbulge < 0) goto skip; //inward bumps dont matter
                }
                else{    //poly is counterclockwise
                    if (Testbulge > 0) goto skip;  //inward bumps dont matter
                }

                if (PtinSeg(p, PtS[i], PtS[i + 1], dBulge) == true){
                    if (wn == 0){
                        return isInside.Inside;
                    }
                    else{
                        return isInside.Outside;
                    }
                }
            skip: ;
            }
            if (wn != 0){
                return isInside.Inside;
            }
            else { return isInside.Outside;
            }
        }




        private static bool SelfIntersectingCheck(Polyline P)
        {
            try
            {
                DBObjectCollection curves = new DBObjectCollection();
                curves.Add(P);
                DBObjectCollection Rs = Region.CreateFromCurves(curves);
            }
            catch (System.Exception)
            {
                return true;
            }
            return false;
        }





        //07-16-07 CCP
        public static bool PtinSeg(Point3d pt, Point2d p1, Point2d p2, double dBulge)
        {
            //Find the center of the arc
            //Center is on the perp. biscector of p1,p2
            //Check if length from center is > radius
            //Check if pt
            if (dBulge == 0) return false;

            double    PDist;
            double m,  mV, dLength;
            double Y, X;
            int iLeft;
            double deltaX = p2.X - p1.X;
            double deltaY = p2.Y - p1.Y;
            if (Math.Abs(deltaX) < 1.0e-8) { deltaX = 0; }//Check for close to zero
            if (Math.Abs(deltaY) < 1.0e-8) { deltaY = 0; }
            double halfseg = 0.5 * Util.MyMath.Length(p1, p2);
            double dRad = halfseg / Math.Sin(Math.Atan(dBulge) * 2);
            double Dist = dRad - (dBulge * halfseg); //tan=opp/adj
            if (deltaX > 0 && deltaY >= 0) { Dist = -Dist; }
            if (deltaX < 0 && deltaY > 0) { Dist = -Dist; }

            if (deltaY == 0){  //Line p1,p2 is horizontal
                Y = p1.Y;
                X = p1.X + 0.5 * deltaX;
                Y = Y - Dist;
            }

            else if (deltaX == 0){ //Line p1,p2 is vertical
                Y = p1.Y + 0.5 * deltaY;
                X = p1.X;
                if (deltaY < 0) Dist = -Dist;
                X = X - Dist;
            }

            else{
                m = deltaY / deltaX;
                mV = -1 / m;  //Invert slope for perpendicular bisector
                X = p1.X + 0.5 * deltaX;  //Find the midpoint
                Y = p1.Y + 0.5 * deltaY;
                //X = X1 + distance / dLength * DX
                //proving for   DeltaX=1 in slope direction
                //slope=DeltaY / 1 =>  DeltaY=slope
                dLength = Math.Sqrt((mV*mV) + 1);
                X = X + (Dist / dLength) * 1;
                Y = Y + (Dist / dLength) * mV;
            }
            Point3d CenPt = new Point3d(X, Y, 0);//To do Sort out 3d
            //Quicktest
            PDist = Util.MyMath.Length(pt, CenPt);
            if (PDist > Math.Abs(dRad))
            { return false; }
            else
            {
                iLeft = MyMath.isLeft(p1, p2, pt);
                if (dBulge > 0){
                    if (iLeft == -1)
                        return true;
                }
                else{ //if (dBulge < 0)
                    if (iLeft == 1)
                        return true;
                }
            }
            return false;

        }



    public bool BBCheck(Point3d P,Entity  Ent){

        Extents3d Ext=Ent.GeomExtents;
        Point3d max=Ext.MaxPoint;
        Point3d min=Ext.MinPoint;
       
        if (P.X > max.X) {return false;}
        if (P.X < min.X) {return false;}
        if (P.Y > max.Y ) {return false;}
        if (P.Y < min.Y ) {return false;}
            return true;
        }



        public isInside PointinPoly(Point3d p, Polyline oPline)
        {
            if (BBCheck(p, oPline) == false){
                return isInside.Outside;
            }
            return(PtInsidePoly(p, oPline));
        }


        [CommandMethod("PTP")]
        public void CreateEmployee()
        {
            Document doc = cad.DocumentManager.MdiActiveDocument;
            Editor ed = doc.Editor;
            Database db = doc.Database;
            PromptEntityResult res = ed.GetEntity("\nSelect the polyline: ");
            PromptPointResult Pr = ed.GetPoint("\nPick a point: ");
            Point3d P = Pr.Value;
            if (res.Status != PromptStatus.OK) return;
            if (Pr.Status != PromptStatus.OK) return;
            using (Transaction tr = db.TransactionManager.StartTransaction()){
                Polyline poly = tr.GetObject(res.ObjectId, OpenMode.ForRead, false) as Polyline;
                if (poly == null) return;
                ed.WriteMessage("\nIs clockwise={0}", MyMath.IsClockwise(poly));
                isInside isP = PointinPoly(P, poly);
                ed.WriteMessage("\nPoint is inside poly={0}", isP);
            }
        }

        [CommandMethod("cw")]
        public void iscw()
        {
            Document doc = cad.DocumentManager.MdiActiveDocument;
            Editor ed = doc.Editor;
            Database db = doc.Database;
            PromptEntityResult res = ed.GetEntity("\nSelect the polyline: ");
            DBObjectCollection curves=new DBObjectCollection();

            if (res.Status != PromptStatus.OK) return;
            Transaction tr = db.TransactionManager.StartTransaction();
            try
            {
                Polyline poly = tr.GetObject(res.ObjectId, OpenMode.ForRead, false) as Polyline;
                curves.Add(poly);
                DBObjectCollection Rs = Region.CreateFromCurves(curves);
                if (poly == null) return;
                ed.WriteMessage("\nIs clockwise={0}", MyMath.IsClockwise(poly));
            }
            catch
            {ed.WriteMessage("\nThe poly is selfcrossing.");}
            finally{tr.Dispose();}
        }

    }
}
and a couple of utilities
Code: [Select]
    class MyMath
    {


        public const int Horizontal = 0;
        public const int Vertical = 1;
        public const int Sloped = 2;

        public static double Abs(Double Num)
        {
            if (Math.Abs(Num) < 1.0e-8)
            {
                Num = 0;
                return 0;
            }
            else
            {
                if (Num > 0)
                { return Num; }
                else
                { return -Num; }
            }
        }


        public static int isLeft(Point2d Startpoint, Point2d Endpoint, Point2d P)
        {
            double Ans = ((Endpoint.X - Startpoint.X) * (P.Y - Startpoint.Y)) -
              ((P.X - Startpoint.X) * (Endpoint.Y - Startpoint.Y));
            if (Math.Abs(Ans) < 1.0e-8)
            { return 0; } //P is on the line
            else
            {
                if (Ans > 0)
                { return 1; } //P is left of the line (CW)
                else
                { return -1; } //P is right of the line (CCW)
            }
        }




        public static int isLeft(Point2d Startpoint, Point2d Endpoint, Point3d P)
        {
            double Ans = (Endpoint.X - Startpoint.X) * (P.Y - Startpoint.Y) -
              (P.X - Startpoint.X) * (Endpoint.Y - Startpoint.Y);
            if (Math.Abs(Ans) < 1.0e-8)
            { return 0; } //P is on the line
            else
            {
                if (Ans > 0)
                { return 1; } //P is left of the line (CW)
                else
                { return -1; } //P is right of the line (CCW)
            }
        }



        public static int isLeft(Point3d Startpoint, Point3d Endpoint, Point3d P)
        {
            double Ans = (Endpoint.X - Startpoint.X) * (P.Y - Startpoint.Y -
              (P.X - Startpoint.X) * (Endpoint.Y - Startpoint.Y));
            if (Math.Abs(Ans) < 1.0e-8)
            { return 0; } //P is on the line
            else
            {
                if (Ans > 0)
                { return 1; } //P is left of the line (CW)
                else
                { return -1; } //P is right of the line (CCW)
            }
        }

        //Length

        public static double Length(Point2d P1)
            { return Math.Sqrt((P1.X * P1.X) + (P1.Y * P1.Y)); }

            public static double Length(Point2d P1, Point2d P2)
            {
                double dX = P2.X - P1.X;
                double dY = P2.Y - P1.Y;
                return Math.Sqrt((dX * dX) + (dY * dY));
            }

            public static double Length(Point3d P1)
            {
                return Math.Sqrt((P1.X * P1.X) + (P1.Y * P1.Y) + (P1.Z * P1.Z));
            }


        public static double Length(Point3d P1, Point3d P2)
            {
                double dX = P2.X - P1.X;
                double dY = P2.Y - P1.Y;
                double dZ = P2.Z - P1.Z;
                return Math.Sqrt((dX * dX) + (dY * dY) + (dZ * dZ));
            }




        //http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm#orientation2D_polygon()
       
   public static bool     IsClockwise(Polyline oPline)
        {
            // first find rightmost lowest vertex of the polygon
            int rmin = 0;
            int n = oPline.NumberOfVertices;
            if (n == 2) { return false; }

            Point2dCollection V = new Point2dCollection();
            for(int i=0;i<n;i++)
                {V.Add(oPline.GetPoint2dAt(i));}
                V.Add (V[0]);
                V.TrimToSize();
            double xmin = V[0].X;
            double ymin = V[0].Y;
           
           
            for (int i = 1; i < n; i++)
            {
                if (V[i].Y > ymin)
                    continue;
                if (V[i].Y == ymin)
                {    // just as low
                    if (V[i].X < xmin)   // and to left
                        continue;
                }
                rmin = i;          // a new rightmost lowest vertex
                xmin = V[i].X;
                ymin = V[i].Y;
            }

            // test orientation at this rmin vertex
            // ccw <=> the edge leaving is left of the entering edge
         
            if (rmin == 0)
            {
                if (MyMath.isLeft(   V[n - 1], V[1]  , V[0])>0)
                { 
                return true; }
            }
            else
            {
                if (MyMath.isLeft(V[rmin - 1], V[rmin + 1] , V[rmin]) > 0)
                {
                    return true;}
                }
       //If we get to here it's false
            return false;
        }





        public static bool isZero(double num)
        {
            //if( (Math.Abs(num) < 0.00000001)
            if (num > -0.00000001)
            {
                if (num < 0.00000001)
                {
                    num = 0;
                    return true;
                }
            }
            return false;

        }
        public static bool Zero(double num)
        {
            if (Math.Abs(num) < 0.00000001)
            {
                num = 0;
                return true;
            }
            return false;
        }

So Vba is verbose ha?

Bryco

  • Water Moccasin
  • Posts: 1883
Re: Point in poly
« Reply #1 on: July 17, 2007, 01:04:22 AM »
the isclockwise function from softsurfer is also worthy of note as it is so fast. Of course if you do have a crossover at the lowest rightest vertex then the answer will be wrong.

LE

  • Guest
Re: Point in poly
« Reply #2 on: July 17, 2007, 09:59:50 AM »
hi Bryco;

looks like you have been doing a lot of homework.

so, you just want to have a function to test if a point is inside of a polyline/or a curve closed object no?

i posted in the show your stuff the function isPointInside() in c# - the only condition i did not include was to ignore self intersection, i'll look at that...

http://www.theswamp.org/index.php?topic=17505.0
« Last Edit: July 17, 2007, 10:01:56 AM by LE »

LE

  • Guest
Re: Point in poly
« Reply #3 on: July 17, 2007, 11:12:20 AM »

Bryco

  • Water Moccasin
  • Posts: 1883
Re: Point in poly
« Reply #4 on: July 17, 2007, 11:34:02 AM »
I'll mess with that tonight, (More likely a week) thanks Luis.
If you draw an hour glass closed poly with straight lines with the intersection in the middle, cad will give the poly a zero area, makes me think it is using similar techniques.

LE

  • Guest
Re: Point in poly
« Reply #5 on: July 17, 2007, 11:39:52 AM »
If you draw an hour glass closed poly with straight lines with the intersection in the middle, cad will give the poly a zero area, makes me think it is using similar techniques.

?

by using testPointInside() it will stop and show a message that it is a self intersection curve.

Bryco

  • Water Moccasin
  • Posts: 1883
Re: Point in poly
« Reply #6 on: July 18, 2007, 12:47:25 AM »
Luis I'm having trouble getting that to work.

LE

  • Guest
Re: Point in poly
« Reply #7 on: July 18, 2007, 01:06:43 PM »
Luis I'm having trouble getting that to work.

Is something wrong with the code?

BTW, I use A2007.