Author Topic: Point is Inside polyline (Now with Bulges)  (Read 22327 times)

0 Members and 1 Guest are viewing this topic.

ymg

  • Guest
Re: Point is Inside polyline
« Reply #30 on: April 25, 2015, 11:53:21 AM »
Lee,

For polylines with bulge, it gets messy.

You can approximate the curve part by adding straight segments.

Code: [Select]
(defun approxpoly (en / px fe)
      (setq px (* 0.75 (acet-geom-pixel-unit))
            fe (acet-list-remove-adjacent-dups
  (acet-geom-object-point-list en (/ px 2.0)))
      )    
   )

But this is bound to create problems when you are on the
boundary.

ymg

GP

  • Newt
  • Posts: 83
  • Vercelli, Italy
Re: Point is Inside polyline
« Reply #31 on: April 25, 2015, 01:58:20 PM »
I don't think this will work with self intersecting polylines.

Also what happen If point P is directly on a vertex of the polyline?

What about polylines with arc segments?


True...  but I think also the other codes posted  :-)

This code should works with self intersecting polylines and also If point P is directly on a vertex/segment of the polyline:

Code - Auto/Visual Lisp: [Select]
  1. (defun inside_p (:p :Lst / far_p cross on)
  2.     (setq far_p (mapcar '+ '(1.0 1.0 0.0) (getvar "extmax")))
  3.     (setq cross 0)
  4.     (if (not (member :p :Lst))
  5.         (mapcar
  6.             '(lambda (a b)
  7.                  (if (inters :p far_p a b) (setq cross (1+ cross)))
  8.                  (if (equal (+ (distance :p a) (distance :p b)) (distance a b) 1e-8) (setq on t))
  9.              )
  10.             :Lst (cdr :Lst)
  11.         )
  12.         (setq on t)
  13.     )
  14.     (or (not (zerop (rem cross 2))) on)
  15. )
« Last Edit: April 25, 2015, 06:29:26 PM by GP »

ymg

  • Guest
Re: Point is Inside polyline
« Reply #32 on: April 25, 2015, 05:03:35 PM »
Gian Paolo,

The winding number method works with self-intersecting polylines
and also if your point is on a vertex.

It does not work with bulges though!

ymg

Quote
Benchmarking .............Elapsed milliseconds / relative speed for 1024 iteration(s):

    (PTINPOLY_P1 P L)..........1279 / 1.52 <fastest>
    (INSIDE_P P L).............1934 / 1.01
    (VK_ISPOINTINSIDE P L).....1950 / 1.00 <slowest>

Images below from:  http://geomalgorithms.com/a03-_inclusion.html
« Last Edit: April 25, 2015, 05:58:59 PM by ymg »

GP

  • Newt
  • Posts: 83
  • Vercelli, Italy
Re: Point is Inside polyline
« Reply #33 on: April 25, 2015, 07:15:57 PM »

Thank you ymg, I realized now, but... seems it does not work always.


roy_043

  • Water Moccasin
  • Posts: 1895
  • BricsCAD 18
Re: Point is Inside polyline
« Reply #34 on: April 26, 2015, 10:00:43 AM »
@ GP:
Further, one must decide whether a point on the polygon's boundary is inside or outside. A standard convention is to say that a point on a left or bottom edge is inside, and a point on a right or top edge is outside. This way, if two distinct polygons share a common boundary segment, then a point on that segment will be in one polygon or the other, but not both at the same time. This avoids a number of problems that might occur, especially in computer graphics displays.

ymg

  • Guest
Re: Point is Inside polyline
« Reply #35 on: April 26, 2015, 10:02:54 AM »
Gian Paolo,

Yes there seems to be some inconsistencies when we are directly
on the line.

Maybe due to some rounding error in the onleft routine.

Thanks,

ymg

Code: [Select]
(defun PtInPoly_p1 (p l / i l p1 p2 v wn y)

   ;; Returns t if point p is strictly to the left of v1->v2   by ymg         ;
   (defun isleft (p v1 v2 / xp yp)
       (setq  xp (car  p)  yp (cadr  p))
       (minusp (- (* (- (cadr v1) yp) (- (car v2) xp)) (* (- (car v1) xp) (- (cadr v2) yp))))
   )
 
   
 
   
   ;; Initialize the loop.                                                    ;
   (setq p1 (car  l)
  y (cadr p)
         wn 0                               
   )
   
   ;; Loop through all edges of polygon.                                      ;
   (foreach p2 (cdr l)
      (if (<= (cadr p1) y)
(if (> (cadr p2) y)                         
    (if (isleft p p1 p2) (setq wn (1+ wn)))

         (if (<= (cadr p2) y)
    (if (not (isleft p p1 p2)) (setq wn (1- wn)))
)
      )
      (setq p1 p2)
   )
   (not (zerop wn));; Returns t if wn not equal to zero.                      ;
)


 ;   // loop through all edges of the polygon                           
 ;   for (int i=0; i<n; i++) {   // edge from V[i] to  V[i+1]           
 ;       if (V[i].y <= P.y) {          // start y <= P.y                 
 ;          if (V[i+1].y  > P.y)      // an upward crossing             
 ;               if (isLeft( V[i], V[i+1], P) > 0)  // P left of  edge   
 ;                    ++wn;            // have  a valid up intersect     
 ;       }                                                               
 ;       else {                        // start y > P.y (no test needed)
 ;           if (V[i+1].y  <= P.y)     // a downward crossing           
 ;                if (isLeft( V[i], V[i+1], P) < 0)  // P right of  edge
 ;                    --wn;            // have  a valid down intersect   
 ;       }                                                               
 ;   }                                                                   
 ;   return wn;                                                         
 ;}                                                                     


(defun listpol (en / i p l) 
  (setq i (if (vlax-curve-IsClosed en)
       (vlax-curve-getEndParam en)
     (+ (vlax-curve-getEndParam en) 1)
  )
  )      
  (while (setq p (vlax-curve-getPointAtParam en (setq i (1- i))))
      (setq l (cons (trans p 0 1 ) l))
  )
)


(defun c:test (/ en l p)
   (princ "\nSelect Polyline Under Test: ")
   (setq en (car (entsel))
  l (listpol en)
  l (cons (last l) l)
   )
   (setvar 'PDSIZE -1.5)
   (while (setq p (getpoint "\nPick a Point: "))
      (if (ptinpoly_p1 p l)
         (entmakex (list (cons 0 "POINT") (cons 10 p) (cons 62 1)))
(entmakex (list (cons 0 "POINT") (cons 10 p) (cons 62 7)))
      )
   )
   (princ)

« Last Edit: April 26, 2015, 10:15:22 AM by ymg »

ymg

  • Guest
Re: Point is Inside polyline
« Reply #36 on: April 26, 2015, 03:14:46 PM »
Roy,

Here I've modified on left by removing the minusp test,
and modified ptinpoly_p to test exactly like Sunday's
routine.

Works more or less OK as by his description.

I do run into some roundoff error on the left side
of the test poly.

Code: [Select]
(defun PtInPoly_p5 (p l / i l p1 p2 v wn y)

   ;; Returns t if point p is strictly to the left of v1->v2   by ymg         ;
   (defun isleft (p v1 v2 / xp yp)
       (setq  xp (car  p)  yp (cadr  p))
       (- (* (- (cadr v1) yp) (- (car v2) xp)) (* (- (car v1) xp) (- (cadr v2) yp)))
   )
 
   
 
   
   ;; Initialize the loop.                                                    ;
   (setq p1 (car  l)
  y (cadr p)
         wn 0                               
   )
   
   ;; Loop through all edges of polygon.                                      ;
   (foreach p2 (cdr l)
      (if (<= (cadr p1) y)
(if (> (cadr p2) y)                         
    (if (> (isleft p p1 p2) 0) (setq wn (1+ wn)))

         (if (<= (cadr p2) y)
    (if (< (isleft p p1 p2) 0) (setq wn (1- wn)))
)
      )
      (setq p1 p2)
   )
   (not (zerop wn));; Returns t if wn not equal to zero.                      ;
)


roy_043

  • Water Moccasin
  • Posts: 1895
  • BricsCAD 18
Re: Point is Inside polyline
« Reply #37 on: April 27, 2015, 09:23:57 AM »
@ ymg:
I was trying to replicate your experiment and while doing so came across some strange results:
Code: [Select]
(setq ptLst '((0.0 0.0 0.0) (104.0 0.0 0.0) (104.0 72.0 0.0) (0.0 72.0 0.0)))

(GEO_GEOM_POINTINSIDEPOINTLIST_P '(1.0 1.0 0.0) ptLst) => T
(PtInPoly_p1                     '(1.0 1.0 0.0) ptLst) => T
(PtInPoly_p2                     '(1.0 1.0 0.0) ptLst) => T
(PtInPoly_p3                     '(1.0 1.0 0.0) ptLst) => T
(PtInPoly_p5                     '(1.0 1.0 0.0) ptLst) => nil

(GEO_GEOM_POINTINSIDEPOINTLIST_P '(-1.0 0.0 0.0) ptLst) => nil
(PtInPoly_p1                     '(-1.0 0.0 0.0) ptLst) => T
(PtInPoly_p2                     '(-1.0 0.0 0.0) ptLst) => T
(PtInPoly_p3                     '(-1.0 0.0 0.0) ptLst) => T
(PtInPoly_p5                     '(-1.0 0.0 0.0) ptLst) => nil

roy_043

  • Water Moccasin
  • Posts: 1895
  • BricsCAD 18
Re: Point is Inside polyline
« Reply #38 on: April 27, 2015, 09:40:30 AM »
It is interesting to note that where http://geomalgorithms.com/a03-_inclusion.html mentions that, by convention, points on the left and bottom edges are considered inside the polyline, ymg's code favors the right and bottom edges instead and the code I have 'ported' favors the left and top edges.

ymg

  • Guest
Re: Point is Inside polyline
« Reply #39 on: April 27, 2015, 11:13:27 AM »
Roy,

This is a finicky! beast.

It is somewhat unsastisfying for example that a point of the polyline
would be considered outside.

ymg

roy_043

  • Water Moccasin
  • Posts: 1895
  • BricsCAD 18
Re: Point is Inside polyline
« Reply #40 on: April 28, 2015, 05:39:36 AM »
@ ymg:
I now see that your code requires a different point list than mine. For your code the last point in the list must be equal to the first point. The strange results reported in post #37 are all due to this oversight. My apologies for that.

The same test with correct point lists:
Code: [Select]
(setq ptLst    '((0.0 0.0 0.0) (104.0 0.0 0.0) (104.0 72.0 0.0) (0.0 72.0 0.0)))
(setq ptLstAlt '((0.0 0.0 0.0) (104.0 0.0 0.0) (104.0 72.0 0.0) (0.0 72.0 0.0) (0.0 0.0 0.0)))

(GEO_GEOM_POINTINSIDEPOINTLIST_P '(1.0 1.0 0.0) ptLst)    => T
(PtInPoly_p1                     '(1.0 1.0 0.0) ptLstAlt) => T
(PtInPoly_p2                     '(1.0 1.0 0.0) ptLstAlt) => T
(PtInPoly_p3                     '(1.0 1.0 0.0) ptLstAlt) => T
(PtInPoly_p5                     '(1.0 1.0 0.0) ptLstAlt) => T
(vk_IsPointInside                '(1.0 1.0 0.0) ptLstAlt) => T

(GEO_GEOM_POINTINSIDEPOINTLIST_P '(-1.0 0.0 0.0) ptLst)    => nil
(PtInPoly_p1                     '(-1.0 0.0 0.0) ptLstAlt) => nil
(PtInPoly_p2                     '(-1.0 0.0 0.0) ptLstAlt) => nil
(PtInPoly_p3                     '(-1.0 0.0 0.0) ptLstAlt) => nil
(PtInPoly_p5                     '(-1.0 0.0 0.0) ptLstAlt) => nil
(vk_IsPointInside                '(-1.0 0.0 0.0) ptLstAlt) => nil

ymg

  • Guest
Re: Point is Inside polyline
« Reply #41 on: April 28, 2015, 08:44:32 AM »
Roy,

No apologies necessary,

I attach a lot of value to your Input.


ymg

ymg

  • Guest
Re: Point is Inside polyline (Now with Bulges)
« Reply #42 on: April 29, 2015, 09:50:59 AM »
Here I've modified approxpoly to work with a tolerance
instead of the zoom level of your screen.

If we use it to get the point list of a polyline with bulges,
ptinpoly_p returns valid results to within the tolerance
given above.

I attached a test program and a drawings with the
polylines I used for testing.

Would appreciate to know if approxpoly can be ported
to briscad.

ymg

Code - Auto/Visual Lisp: [Select]
  1. ;; Routine to Test if a Point is inside a polylines                           ;
  2. ;; Will work with bulges and Self-Crossing and Overlapping Poly                ;
  3. ;;                                                                            ;
  4.  
  5. (defun c:test (/ *AcadDoc* *error* en l p varl listpol approxpoly ptinpoly_p)
  6.     ;;; Error Handler by ElpanovEvgenyi                                       ;
  7.     (defun *error* (msg)
  8.         (mapcar 'eval varl)
  9.         (if (and msg (not (wcmatch (strcase msg) "*BREAK*,*CANCEL*,*EXIT*")))
  10.            (princ (strcat "\nError: " msg))
  11.         )
  12.         (and *AcadDoc* (vla-endundomark *AcadDoc*))
  13.         (princ)
  14.     )
  15.      
  16.     (setq varl '("CLAYER" "OSMODE" "CMDECHO" "PDMODE" "PDSIZE")
  17.           varl (mapcar (function (lambda (a) (list 'setvar a (getvar a)))) varl)
  18.     )    
  19.      
  20.    (or *AcadDoc* (setq *AcadDoc* (vla-get-activedocument (vlax-get-acad-object))))
  21.  
  22.    ;;                                                                         ;
  23.    ;; listpol     by ymg    (Simplified a Routine by Gile Chanteau)           ;
  24.    ;;                                                                         ;
  25.    ;; Parameter:  en,  Entity Name or Object Name of Any Type of Polyline     ;
  26.    ;;                                                                         ;
  27.    ;; Returns:    List of Points in Current UCS                               ;
  28.    ;;                                                                         ;
  29.    ;; Notes:      On Closed Polyline the Last Vertex is Same as First)        ;
  30.    ;;                                                                         ;
  31.  
  32.    (defun listpol (en / i l)
  33.       (repeat (setq i (fix (1+ (vlax-curve-getEndParam en))))
  34.          (setq l (cons (trans (vlax-curve-getPointAtParam en (setq i (1- i))) 0 1) l))
  35.       )
  36.    )
  37.  
  38.    ;;                                                                         ;
  39.    ;; approxpoly        by ymg       (Derived from Fast Select)               ;
  40.    ;;                                                                         ;
  41.    ;; Parameters: en   Entity Name of a Polyline (Any types)                  ;
  42.    ;;             tol  Tolerance Allowed for Lines Segments.                  ;
  43.    ;;                                                                         ;
  44.    ;; Returns:   A List of Points Replacing Bulges With Straight Segments.    ;
  45.    ;;                                                                         ;
  46.  
  47.    (defun approxpoly (en tol)
  48.       (acet-list-remove-adjacent-dups
  49.            (acet-geom-object-point-list en tol)
  50.       )    
  51.    )
  52.  
  53.    ;;=========================================================================;
  54.    ;;                                                                         ;
  55.    ;;      PtInPoly_p(): Winding number test for a point in a polygon         ;
  56.    ;;                                                                         ;
  57.    ;;      Input:    p = a point                                              ;
  58.    ;;               en = ename or vla-object of a closed polyline.            ;
  59.    ;;                                                                         ;
  60.    ;;      Return:  t, if point is in polyline. (wn is not 0)                 ;
  61.    ;;                                                                         ;
  62.    ;; Original code in C++ by Dan Sunday, Translated by ymg                   ;
  63.    ;; See: http://geomalgorithms.com/a03-_inclusion.html                      ;
  64.    ;;=========================================================================;
  65.  
  66.    (defun ptinpoly_p (p l / i p1 p2 v wn y isleft)
  67.  
  68.       ;; Returns t if point p is strictly to the left of v1->v2   by ymg      ;
  69.       (defun isleft (p v1 v2 / xp yp)
  70.           (setq  xp (car  p)  yp (cadr  p))
  71.           (- (* (- (cadr v1) yp) (- (car v2) xp)) (* (- (car v1) xp) (- (cadr v2) yp)))
  72.       )
  73.    
  74.       ;; Initialize the loop.                                                 ;
  75.       (setq p1 (car  l)
  76.              y (cadr p)
  77.             wn 0                              
  78.       )
  79.    
  80.       ;; Loop through all edges of polygon.                                   ;
  81.       (foreach p2 (cdr l)
  82.          (if (<= (cadr p1) y)
  83.             (if (> (cadr p2) y)                          
  84.                (if (> (isleft p p1 p2) 0) (setq wn (1+ wn)))
  85.             )
  86.             (if (<= (cadr p2) y)
  87.                (if (< (isleft p p1 p2) 0) (setq wn (1- wn)))
  88.             )  
  89.          )
  90.          (setq p1 p2)
  91.       )
  92.       (not (zerop wn)) ;; Returns t if wn not equal to zero.                  ;
  93.    )
  94.  
  95.  
  96.    (princ "\nSelect Polyline Under Test: ")
  97.    (setq en (car (entsel))
  98.          ;l (listpol en)           ;   Use this if poly has no bulges         ;
  99.           l (approxpoly en 0.001)  ;   Use this if poly has bulges.           ;
  100.    )
  101.    
  102.    (setvar 'PDMODE 34) (setvar 'PDSIZE -3.5)
  103.    (while (setq p (getpoint "\nPick a Point: "))
  104.       (if (ptinpoly_p p l)
  105.          (entmakex (list (cons 0 "POINT") (cons 10 p) (cons 62 1)))
  106.          (entmakex (list (cons 0 "POINT") (cons 10 p) (cons 62 7)))
  107.       )
  108.    )
  109.    (*error* nil)
  110. )  
  111.  
« Last Edit: April 29, 2015, 10:01:51 AM by ymg »

ronjonp

  • Needs a day job
  • Posts: 7526
Re: Point is Inside polyline (Now with Bulges)
« Reply #43 on: April 29, 2015, 10:30:33 AM »
THIS could be used to trace the polyline rather than relying on ET functions.
Here is a quick test with your code above.

Windows 11 x64 - AutoCAD /C3D 2023

Custom Build PC

ymg

  • Guest
Re: Point is Inside polyline (Now with Bulges)
« Reply #44 on: April 29, 2015, 10:39:11 AM »
Ronjomp,

Looks good!! to me.

I'll do some test to see what happen when
there is no or little bulges.

ymg