Topic: check the supplied point is inside a polygon

[sub]

Hello all,

I want to check the supplied co-ordinate is inside a polygon or not.
polygon is not rectangle or square.

Any Ideas

[Anonymous]

Hi sub,
One way is to radiate from your pt to each vertices of your polygon. If the sum of the subtended angles adds up to 360; the point is strictly in. Another way is to put a point at the co-ordinate to be tested, then use the list of vertex points and with the help of ssget "WP" check if it is picked! Entdel the punctual object when over. Also, using c:Bpoly to see if pt in, returns T and creates a polyline. Others might have hints or ready-made routines! HTH

[Serge J. Gianolla]

$hit or should i be polite and say merde, the loggin' bizo plays up again, anyway sub, above 3 methods should take you somewhere.

[SMadsen]

I prefer the method of shooting out rays and count the number of times they hit the edges. If even, the point is outside, otherwise it's inside. Drawbacks:
1. If point is located on an edge or coincides with a vertex the result can not be trusted.
2. Rays may hit vertices instead of edges but averaging hits from numeral rays can give a pretty good certainty.

Serge, won't any point in- or outside a polygon return 360 from the sum of angles? For example:

Code: [Select]

(defun PtPolyTest (/ plst)
  (setq pl  (car (entsel))
        pll (entget pl)
  )
  (foreach n pll
    (and (= (car n) 10) (setq plst (cons (cdr n) plst)))
  )
  (setq pt   (getpoint "\nPoint: ")
        ang  0.0
        plst (reverse plst)
        plst (append plst (list (car plst)))
  )
  (while (cdr plst)
    (setq ang  (+ ang (- (angle pt (cadr plst)) (angle pt (car plst))))
          plst (cdr plst)
    )
  )
  (terpri)(princ (angtos ang))(terpri)
  (equal ang 0.0 1e-6)
)

[CAB]

This stuff is over my head but I collected this some time ago.

Re: Point inside the perimeter of a closed polyline ?
Subject: Re: Point inside the perimeter of a closed polyline ?
From: kboutora@francemel.com (Kamal Boutora)
Date: Thu, 16 Dec 1999 15:40:45 +0000
Newsgroups: autodesk.autocad.customization
In article <838ik5$kje17@adesknews2.autodesk.com>,
not.robertb@mwengineers.com says...
> Could you post code, or a link to some?
>

  From the faq of comp.graphics.algorithms (a very interesting faq that
  every programmer should have IMHO).

    you can look at the latest HTML version at either of these two sites:

http://www.exaflop.org/docs/cgafaq
http://www.cis.ohio-state.edu/hypertext/faq/usenet/graphics/algorithms-
faq/faq.html

    The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>

    int pnpoly(int npol, float *xp, float *yp, float x, float y)
    {
      int i, j, c = 0;
      for (i = 0, j = npol-1; i < npol; j = i++) {
        if ((((yp<=y) && (y<yp[j])) ||
             ((yp[j]<=y) && (y<yp))) &&
            (x < (xp[j] - xp) * (y - yp) / (yp[j] - yp) + xp))

          c = !c;
      }
      return c;
    }




Subject: Re: About Polyline
Author: Horst Kraemer <horst.kraemer@snafu.de>
Date Posted: Oct 2 1999 9:15:21:000AM


On Sat, 02 Oct 1999 09:34:42 GMT, philipma@ms5.hinet.net (Philip Ma)
wrote:

> Hi,
>
> Maybe the questions are stupid, but I need the answers. Thank you in
> advance!
>
> There are some questions about 2-D closed polyline as following:
>   1. How to know a polyline is crisscross (intersection itself) ?
>   2. If the polyline is not crisscross, how to know the polyline is
> clockwise or counterclockwise?
>   3. If the polyline is not crisscross, how to know a point is inside,
> outside or on the polyline?
>   4. If the polyline is not crisscross, how to know the area of the
> polyline?

In the order of ascending computational cost:

If  P0=(x0,y0),P1=(x1,y1),....,PN=(xN,yN)=P0

is your closed polygon consisting of a sequence of N distinct points
P1,..,PN then the _oriented_ area is 1/2 of the sum of

   x_i*y_{i+1}-y_i*x_{i+1}

for i from 0 to N-1, i.e. for all adjacent pairs of vertices of the
closed polygon. This area may be positive or negative, according to
the orientation of the sequence relative to the coordinate system. In
any case its magnitude is the surface area.

If the coordinate system is oriented in a way that turning the
positive x-axis counterclockwise by 90deg will turn it into the
positive y-axis, then a positive area means that the polygon is
oriented counterclockwise. In fact "positive" doesn't mean
"counterclockwise". It just means "the same orientation as the
coordinate system". In fact in 2D geometry "clockwise" has no meaning.
If you look at a plane from the other side "clockwise" will turn into
"counterclockwise". You have to define additionally which side of the
plane is the "upside" and this definition can only be done by "looking
at it" from 3D. (Sorry for the side-step...)

There are many methods to test if a point of the plane is inside a
simple polygon. None of them is technical trivial.

The method I prefer is this one. It has been published by the german
mathematician Ahlemeyer, although he may not be the first who invented
it.  Let's assume for simplicity that the point to be tested is the
origin.

Divide the plane into 3 regions.

S0 consists only of the origin (0,0)

S+ is the "upper sector", consisting of all points with y>0 + all
points (x,0) with x>0

S- is the "lower sector", consisting of all points with y<0 + all
points (x,0) with x<0.

Every point of the plane belongs to exactly one sector.

Initialize a counter with 0.

Now make a round trip as before through all pairs of vertex points. If
any vertex is in S0 then (0,0) is a a vertex of the polygon and you
make your decision depending on the detail if an edge/vertex is
"inside" or "outside" the polygon.

Otherwise
  if two adjacent vertices (x_i,y_i),(x_{i+1},y_{i+1}) are in
  the same sector, do nothing.
  else if they are in two distinct sectors then calculate

      d =  x_i*y_{i+1) - y_i*x_{i+1)

  if d=0 then the origin is situated on the edge and you
  make your decision
  else if d>0 then add 1 to the counter
  else (d<0) subtract 1 from the counter

If no final decision could be taken in the middle then the point is
outside if the counter is ZERO, otherwise the point is outside. (Then
the counter is either 2 or -2).


Testing if a polygon is simple (not crossing itself) is probably the
most difficult task in terms of computation time. The crude and
straightforward approach is to test for each pair of edges if they
intersect or not. This will only be a problem for _very_ many
vertices.

Some time ago we hat a long thread about this problem in
sci.math.research, but nobody could come up with a _working_ algorithm
which was not O(n^2), i.e. which did need significantly less than
n^2/2 distinct tests. There were a lot of "convincing" proposals but
none of them was waterproof.


Regards
Horst

[CAB]

If it code you want look here.
http://www.theswamp.org/lilly.pond/CAB/Pline%20Point%20Inside.lsp

[SMadsen]

Excellent links. Thanks CAB

The non-trivial suggestions by Horst Kraemer may be even more non-trivial when one considers polylines with bulged segments.
Which is why I prefer the ray-edge intersection test, despite a few drawbacks. Manual intersection with lines/arcs isn't that hard to do but  IntersectWith already handles bulges and spits out a ready-to-use list of intersections.

[CarlB]

You can get a routine named "GE_PtInPoly" at http://www.4d-technologies.com/techcenter/index.htm  It uses the "ray crossing athe boundary" method.

[SMadsen]

CarlB, very good routine for pline with straight segments. Nice touch with the decision of how to treat points located on segments.

It doesn't work for bulged plines, though.

Here's a very quick outline of a "count IntersectWith" method. Ok, so it's obviously a reinvention of the wheel, as the good stuff CAB posted also shows an IntersectWith routine .. but hey, my routine is animated!  
Adjust the length of the "ray" to fit a possible test .. here's it's just set to 100.0  

Code: [Select]

(defun PtInPoly (/ clst apt delta ints oLine oPline pl pline someAngle)
  (vl-load-com)
  (cond
    ((setq pline (car (entsel "\nPolyline: ")))
     ;; NO error checking .. assume a lwpoly is selected!
     (setq oPline       (vlax-ename->vla-object pline)
           ;; pick point to investigate
           aPt       (getpoint "\nPick point: ")
           ;; initialize starting angle for "ray"
           someAngle 0.0
           ;; set delta angle for rotation of "ray"
           delta     (/ (* 2 pi) 8)
     )
     (cond
       ;; start by making a "ray" shooting out from aPt
       ;; (should use vla-add but what the heck)
       ;; (should also use a real RAY and manipulate
       ;;  unit vector .. but what the heck)
       ((entmake (list '(0 . "LINE") (cons 10 aPt)
                       (cons 11 (polar aPt 0.0 100.0))
                 )
        )
        (setq oline (vlax-ename->vla-object (entlast)))
        (repeat 8
          (cond
            ;; get intersection points with pline and "ray"
            ((setq ints (vla-intersectwith oPline oLine acExtendNone))
             ;; should be using vlax-...-u-bounds and all that to check
             ;; safearray (but what the heck, this is quicker to write)
             (if (not (vl-catch-all-error-p (setq ints (vl-catch-all-apply
                                     'vlax-safearray->list
                                     (list (vlax-variant-value ints)))
                        )
                      )
                 )
               ;; just put nil for uneven number of hits and T for even
               (setq clst (cons (not (zerop (rem (length ints) 2.0))) clst))
             )
            )
          )
          ;; pause to see the "ray" move around
          (while (not (grread nil 10)))
          ;; move endpoint of "ray"
          (vla-put-endpoint oLine (vlax-3D-point
              (polar aPt (setq someAngle (+ delta someAngle)) 100.0)
            )
          )
        )
        (vla-delete oLine)
        (vlax-release-object oLine)
        (vlax-release-object oPline)
       )
     )
    )
  )
  ;; .. and a lazy decoding of result:
  (cond ((not (member 'nil clst))(princ "Inside"))
        ((not (member 'T clst))(princ "Outside"))
        ((princ "Probably on an edge or vertex")))
  (terpri)
  clst
)

[Serge J. Gianolla]

Yep Stig you got me on the Sigma of angles so I believe it is my turn to report a flaw in the intersect method . When working with circles or ellipses, one can count if object is intersected once, point is in; if intersecting twice point is outside. But with odd shape plines no such luck. Draw a pline in the outline of letter E, then pick point to evaluate in middle branch, with your method it will sometimes cut a number of even edges and other times odd number of sides. Not full proof! How 'bout ssget "WP"?

[CarlB]

Serge you are correct about sum of the angles of all vertices as a test for point in poly, you just might not have stated the correct angle.  Various strategies for testing if point is in a polygon is discussed at http://www.acm.org/pubs/tog/editors/erich/ptinpoly/. Yeah Stig, probably still doesn't address bulges.

[SMadsen]

Without having seen the discussion on the link CarlB posted I'll take a wild guess that the angle method is good if "sign" is considered (i.e. direction)?

Serge, I don't doubt the intersect method will fail for some shapes. However, I can't find where it'll fail in the case you mention unless the ray hits a vertice or coincides with an edge. I'm sure it's just me who can't see it and I take your word for it. But can you make a small sketch to illustrate?

[Anonymous]

Hey Stig, you made me doubt myself about sigma of angles adding to 360! You cheeky boy Obviously was not 100% sure of myself Thanks CarlB to believe in me! Seriously Stig, let us have a rectangle and label the bottom right corner 1, the top right 2, TL=3 and BL=4 and the radiating point is A and located below the bottom base of rectangle. Now, we have to accept the precept that the angle cannot be greater than 180 deg. It can be 180, meaning that A is precisely on the line connecting the points or less than 180! Each angle is to be tested if greater than 180. So, we have angle 1A2, 2A3, 3A4 and 4A1. The last one being larger than 180, we use alpha=360-4A1 where alpha is now clockwise [as AutoCAD works] while the others are CCW! If A is inside, angles add up to 360 if outside resulting in close to 0 for rounding off errors. I am not at the office yet but when there i'll post something with bpoly.
Since there was no more input from sub, it is hard to gear in the right direction! Is the polygon drawn using 3DPOLY, UCS with world, are there any object in the way that could confuse the ray method...
Serge

[Serge J. Gianolla]

Sub was not specific, a polygon does not necessarily mean drawn with pline! This method will help user by picking a point and creates a pline if a bounded perimeter is found. It can be massaged to have the point already preset as opposed to picking. This is not gonna' work if there are other objects obstructing in the sought polygon!

Code: [Select]

(defun upoint (bit kwd msg def bpt / inp)
  (if def ;check for a default
    (setq pts (strcat (rtos (car def))
     ","
     (rtos (cadr def)) ;formats X,Y 2D pt as string
     (if ;formats 3D ,Z if supplied and FLATLAND off
(and (caddr def) (= 0 (getvar "FLATLAND")))
(strcat "," (rtos (caddr def)))
""
     ) ;_ end of if
     ) ;if&strcat
 msg (strcat "\n" msg " <" pts ">: ") ;string them with default
 bit (* 2 (fix (/ bit 2))) ;a default and no null bit code conflict so
    ) ;this reduces bit by 1 if odd, to allow null
    (setq msg (strcat "\n" msg ": ")) ;or without
  ) ;if a default was specified
  (initget bit kwd)
  (setq inp (if bpt ;check for base point
     (getpoint msg bpt) ;and use it
     (getpoint msg) ;but not if nil
   ) ;_ end of if
  ) ;setq&if
  (if inp
    inp
    def
  ) ;evaluate results and return proper value
) ;defun
(defun InorOut (/ BDRY P BPSET)
  (if (setq BfSET (ssget "i"))
    (princ "\nUsing current existing boundary set.")
  ) ;_ end of if
   (setq pStarter (upoint 6 "" "Select internal point" nil nil) )
     (cond
       ((setq BDRY (BPOLY pStarter BfSET))
;; (command "_.ERASE" BDRY "")
       )
       (T (princ "\nValid boundary not found."))
     ) ;endcond
  (princ)
)