[原创]：射线法及角度法求解点是否在曲线内[角度法在4楼]

fools

今天先贴射线法,明天再贴角度法,均支持各种曲线的测试
[pcode=lisp,true]
(DEFUN InCurve (pt ent / COUNT TMPRAY)
  (IF (EQUAL pt (VLAX-CURVE-GETCLOSESTPOINTTO ent pt) 1E-6)
    (ALERT "***点在线上***")
    (PROGN (SETQ ent        (VLAX-ENAME->VLA-OBJECT ent)
                 TmpRay        (VLAX-ENAME->VLA-OBJECT
                          (ENTMAKEX (LIST '(0 . "RAY")
                                          '(100 . "AcDbEntity")
                                          '(100 . "AcDbRay")
                                          (CONS 10 pt)
                                          (CONS 11 (MAPCAR '+ pt '(1 0)))
                                          ;;相当于(polar pt 0 1)
                                    )
                          )
                        )
                 pt        (VLAX-3D-POINT pt)
                 Count        0
           )
           ;;可根据需要调整扫描角度提高检索速度,本程序采用12度
           (REPEAT 30
             (VLA-ROTATE TmpRay pt (/ PI 15))
             ;;下面这句可以看到扫描过程,实际使用时可以注释掉
             ;;(COMMAND "delay" 50)
             (IF (ZEROP
                   (REM        (LENGTH        (VLAX-INVOKE ent 'INTERSECTWITH TmpRay ACEXTENDNONE)
                        )
                        6
                   )
                 )
               (SETQ Count (1- Count))
               (SETQ Count (1+ Count))
             )
           )
           (VLA-DELETE TmpRay)
           (IF (MINUSP Count)
             (ALERT "***点在线外***")
             (ALERT "***点在线内***")
           )
    )
  )
  (PRINC)
)
[/pcode]

[pcode=lisp,true]
(DEFUN C:test (/ Curve Pt)
  (IF (SETQ Curve (CAR (ENTSEL "\n选择一条曲线:")))
    (WHILE (SETQ Pt (GETPOINT "\n点取测试点:"))
      (InCurve Pt Curve)
    )
  )
  (PRINC)
)
[/pcode]

snoopychen

：）

好像舟自横兄最近也研究了相关的问题。

也找出几个类似的函数继续好好学习。射线法的居多。不好意思，以前测试过，发这篇帖子的时候没有再测试一遍，不知道是否都贴全了。

Luis E的，这篇里面有好些函数
http://groups.google.com/group/a ... nside+curve&rnum=5#
[pcode=lisp,true]
;;; 9/4/03 - 6:36 AM 9/8/2003 by Luis Esquivel
;;; http://www.draftteam.com

(defun point-inside-region-p
       (vla_poly pt aid-pt pts / vla_line lst_value param result)
  (setq
    result
     (cond
       ((and (setq param (vlax-curve-getparamatpoint vla_poly pt))
      (eq (type (- param (fix param))) 'real))
nil)
       ((and (not (vl-position (list (car pt) (cadr pt)) pts))
      (not
        (vl-catch-all-error-p
   (setq lst_value
   (vl-catch-all-apply
     'vlax-safearray->list
     (list (vlax-variant-value
      (vla-intersectwith
        vla_poly
        (setq vla_line
        (vla-addline
          (vla-objectidtoobject
            (vla-get-document
       vla_poly)
            (vla-get-ownerid
       vla_poly))
          (vlax-3d-point pt)
          (vlax-3d-point aid-pt)))
        acextendnone))))))))
(fix (/ (length lst_value) 3.0)))))
  (vl-catch-all-apply 'vla-delete (list vla_line))
  (if (numberp result)
    (not (zerop (rem result 2)))))

(defun C:TEST  (/ ent elst pt vla_poly lc uc)
  (if
    (and (setq ent (car (entsel "\nSelect a polyline: ")))
  (eq (cdadr (setq elst (entget ent))) "LWPOLYLINE")
  (eq (vla-get-closed
        (setq vla_poly (vlax-ename->vla-object ent)))
      :vlax-true)
  (setq pt (getpoint "\nTest point: ")))
     (progn
       (vla-getboundingbox vla_poly 'lc 'uc)
       (point-inside-region-p
  vla_poly
  pt
  (polar pt
  (/ pi 2.0)
  (distance (vlax-safearray->list lc)
     (vlax-safearray->list uc)))
  (mapcar
    'cdr
    (vl-remove-if-not
      (function (lambda (item) (eq (car item) 10)))
      elst))))))
(princ)
[/pcode]

John Uhden的程序
http://groups.google.com/group/a ... =9#11d71de3be9338e2
[pcode=lisp,true]
;; Updated (03-07-03)
;; Given:
;;-------------------------------------------------------
;; Function to find or create an invisible ray VLA-Object
;; on Layer "0" an to create a global symbol $cv_ray.
;; It will leave the Ray object in the drawing to save time.
(defun @cv_ray ( / e ss)
  (if (setq ss (ssget "X" '((0 . "RAY")(8 . "0")(60 . 1))))
    (setq $cv_ray (vlax-ename->vla-object (ssname ss 0)))
    (setq $cv_ray
      (entmakex
        (list
         '(0 . "RAY")'(100 . "AcDbEntity")'(100 . "AcDbRay")
         '(8 . "0")'(60 . 1)'(10 0.0 0.0 0.0)'(11 1.0 0.0 0.0)
        )
      )
      $cv_ray (vlax-ename->vla-object $cv_ray)
    )
  )
)
;; And:
;;-------------------------------------------------
;; Function originated by Ken Alexander (03-05-03)
;; that is 10X faster than my @cv_parse_list,
;; to group data into triplets.
;; Thanks, Ken!
;;
(defun @cv_triple_up (old / new)
  (while
    (setq new (cons (list (car old)(cadr old)(caddr old)) new)
          old (cdddr old)
    )
  )
  (reverse new)
)
;; And:
;;-------------------------------------------------------------------
;; Function to determine if a point <PIQ> is inside a closed polyline
;; based on the number of intersections found between a ray whose
;; basepoint is the PIQ, and that PIQ is not one of the intersection
;; points.
;; Arguments:
;;   PIQ = 3D Point in WCS
;;   Outer = Outer Polyline VLA-Object
;; Returns:
;;   either T (inside) or nil (on or outside)
(defun @cv_inside (PIQ Outer / Points)
  (vl-load-com)
  (and
    (> (vl-list-length PIQ) 1)
    (vl-every 'numberp PIQ)
    (= (type Outer) 'VLA-Object)
    (vlax-property-available-p Outer 'Closed)
    (= (vla-get-closed Outer) :vlax-true)
    (or $cv_ray (@cv_ray))
    (or (vlax-put $cv_ray "Basepoint" PIQ) T)
    (setq Points (vlax-invoke Outer "IntersectWith" $cv_ray acExtendNone))
    (setq Points (@cv_triple_up Points))
    (= (rem (length Points) 2) 1)
    (not (equal PIQ (vlax-curve-getclosestpointto Object PIQ) 1e-11))
  )
)

;; Then:
;; ... define your point and Outer object ...
(setq Inside (@cv_inside PIQ Outer))
[/pcode]

CAB收集的一些。

[pcode=lisp,true]
;; ==========================================================
;;           Collection of methods to determin if a point
;;             is inside a closed poly line
;; ==========================================================

;|Re: Point inside the perimeter of a closed polyline ?
Subject: Re: Point inside the perimeter of a closed polyline ?
From: "Michael Doerr" <mdoerr@cebacus.de>
Date: Tue, 14 Dec 1999 06:47:25 +0000
Newsgroups: autodesk.autocad.customization
Hi Dominique,

the following code calculates the sum of the angles from a given point to
the points of the polyline and decides wether the point is inside or outside
in that way. 'PointInQuestion' means the point which has to be tested and
'Point_List' is the list of Points which represent the polyline. If the
function returns 'T' the given point is inside your polyline.

Greetings, Michael|;

(defun punktinpolylinie        (pointinquestion
                         point_list
                         /
                        )
  (if (equal 0.0 (pipwinkelsumme pointinquestion point_list) 0.0001)
    nil
    t
  )
)

(defun pipwinkelsumme (pointinquestion                 point_list
                       /            count         p1
                       p2            scheitel         winkeleins
                       winkelzwei
                      )
  (setq        winkeleins 0.0
        scheitel   (car point_list)
        count           1
  )
  (while (< 1 (length point_list))
    (setq p1             (car point_list)
          p2             (cadr point_list)
          point_list (cdr point_list)
          winkelzwei (pipwinkelhilfe pointinquestion p1 p2)
          winkelzwei (if (< 180.0 winkelzwei)
                       (- winkelzwei 360.0)
                       winkelzwei
                     )
          winkeleins (+ winkeleins winkelzwei)
    )
    (setq count (1+ count))
  )
  (setq        winkelzwei (pipwinkelhilfe pointinquestion p2 scheitel)
        winkelzwei (if (< 180.0 winkelzwei)
                     (- winkelzwei 360.0)
                     winkelzwei
                   )
  )
  (+ winkeleins winkelzwei)
)

(defun pipwinkelhilfe (pointinquestion p1 p2 / alpha beta)
  (setq        beta  (angle pointinquestion p1)
        alpha (angle pointinquestion p2)
        alpha (- alpha beta)
  )
  (if (< alpha 0)
    (setq alpha (+ (* 2 pi) alpha))
  )
  (* (/ (float alpha) pi) 180.0)
)

;; ==========================================================
;; ==========================================================


;|Here are a set of LISP functions that do the job for old-style polylines
(I did it a few years ago and it may not be perfect LISP). You will have
to modify the function PLINEVERTICES to be able to use new-style plines.
INPOLY does not calculate bulges, but gives correct results even for
self-intersecting plines. I used the fact that the sum of angles from
the checkpoint to all vertices is equal 0 (360) degrees if the
checkpoint is inside and any othre value if it is outside. The FUZZY
value is to avoid rounding errors.
Code is free to use by anybody

Tom
|;

(setq imr_fuzzy 0.0001)

(defun inpoly (pt en /)
  (if (equal 0.0 (anglesum pt (plinevertices en)) imr_fuzzy)
    nil
    t
  )
)


(defun c:inpoly (/ pt en)
  (setvar "cmdecho" 0)
  (if (= "LWPOLYLINE" (cdr (assoc 0 (entget (setq en (car (entsel)))))))
    (while (setq pt (getpoint "\npoint to check: "))
      (if (inpoly pt en)
        (princ "INSIDE  polyline!")
        (princ "OUTSIDE polyline!")
      )
    )
  )
  (prin1)
)

(defun angd (scheitel p1 p2) (r2d (angr scheitel p1 p2)))

(defun r2d (winkel) (* (/ (float winkel) pi) 180.0))

(defun angr (scheitel p1 p2 / alf bet)
  (setq bet (angle scheitel p1)
        alf (angle scheitel p2)
        alf (- alf bet)
  )
  (if (< alf 0)
    (setq alf (+ (* 2 pi) alf))
  )
  alf
)

(defun plistisplanar (plist / isplanar)
  (setq isplanar t
        pl (cddr plist)
  )
  (while (and isplanar (< 3 (length plist)))
    (while (< 1 (length pl))
      (if (inters (car plist) (cadr plist) (car pl) (cadr pl))
        (setq isplanar nil
              pl nil
              plist nil
        )
      )
      (setq pl (cdr pl))
    )
    (if isplanar
      (setq plist (cdr plist)
            pl    (cddr plist)
      )
      (setq plist nil)
    )
  )
  isplanar
)

(defun anglesum (pt plist / as p1 p2 an1 an2 an)
  (setq as    0.0
        stp   (car plist)
        count 1
  )
  (while (< 1 (length plist))
    (setq p1    (car plist)

          p2    (cadr plist)
          plist (cdr plist)
          an    (angd pt p1 p2)
          an    (if (< 180.0 an)
                  (- an 360.0)
                  an
                )
          as    (+ as an)
    )
    (setq count (1+ count))
  )
  (setq
    an (angd pt p2 stp)
    an (if (< 180.0 an)
         (- an 360.0)
         an
       )
  )
  (+ as an)
)

(defun plinevertices (en / el vl)
  (while (= "VERTEX"
            (cdr (assoc 0
                        (setq el (entget (setq en (entnext
                                                    en
                                                  )
                                         )
                                 )
                        )
                 )
            )
         )
    (setq vl (cons (cdr (assoc 10 el)) vl))
  )
  (if (= (car vl) (car (reverse vl)))
    (reverse (cdr (reverse vl)))
    vl
  )
)

;; ==========================================================
;; ==========================================================

;; Updated (03-07-03)
;; Given:
;;-------------------------------------------------------
;; Function to find or create an invisible ray VLA-Object
;; on Layer "0" an to create a global symbol $cv_ray.
;; It will leave the Ray object in the drawing to save time.
(defun @cv_ray ( / e ss)
  (if (setq ss (ssget "X" '((0 . "RAY")(8 . "0")(60 . 1))))
    (setq $cv_ray (vlax-ename->vla-object (ssname ss 0)))
    (setq $cv_ray
      (entmakex
        (list
         '(0 . "RAY")'(100 . "AcDbEntity")'(100 . "AcDbRay")
         '(8 . "0")'(60 . 1)'(10 0.0 0.0 0.0)'(11 1.0 0.0 0.0)
        )
      )
      $cv_ray (vlax-ename->vla-object $cv_ray)
    )
  )
)
;; And:
;;-------------------------------------------------
;; Function originated by Ken Alexander (03-05-03)
;; that is 10X faster than my @cv_parse_list,
;; to group data into triplets.
;; Thanks, Ken!
;;
(defun @cv_triple_up (old / new)
  (while
    (setq new (cons (list (car old)(cadr old)(caddr old)) new)
          old (cdddr old)
    )
  )
  (reverse new)
)
;; And:
;;-------------------------------------------------------------------
;; Function to determine if a point <PIQ> is inside a closed polyline
;; based on the number of intersections found between a ray whose
;; basepoint is the PIQ, and that PIQ is not one of the intersection
;; points.
;; Arguments:
;;   PIQ = 3D Point in WCS
;;   Outer = Outer Polyline entities names
;; Returns:
;;   either T (inside) or nil (on or outside)
(defun @cv_inside (PIQ Outer / Points)
  (vl-load-com)
  (setq Outer (vlax-ename->vla-object Outer))
  (and
    (> (vl-list-length PIQ) 1)
    (vl-every 'numberp PIQ)
    (= (type Outer) 'VLA-Object)
    (vlax-property-available-p Outer 'Closed)
    (= (vla-get-closed Outer) :vlax-true)
    (or $cv_ray (@cv_ray))
    (or (vlax-put $cv_ray "Basepoint" PIQ) T)
    (setq Points (vlax-invoke Outer "IntersectWith" $cv_ray acExtendNone))
    (setq Points (@cv_triple_up Points))
   
    (= (rem (length Points) 2) 1)

(setq tem (vlax-curve-getclosestpointto Outer PIQ))

    (not (equal PIQ tem 1e-11))
  )
)

(defun c:x()
(setq b (getpoint))
(setq temp (vlax-ename->vla-object (car (entsel))))
(@cv_inside b temp)
)

;; Then:
;; ... define your point and Outer object ...
;(setq Inside (@cv_inside PIQ Outer))
;(not (equal b (vlax-curve-getclosestpointto temp b) 1e-11))

;; ==========================================================
;; ==========================================================


;;-----------------------------------------------
;;             LISTING 3
;;-----------------------------------------------
;;
;; Function to determine if a point is inside or
;; outside a set of closed geometry described
;; in a list.
;; Result: -1 if point is inside, +1 if point
;; is outside of geometry list.
;; -------------------
;; Calling Syntax:
;;
;; (inout DataList Point)
;; -------------------
;;
;; DataList is nested list
;;  (((startpoint) (endpoint) "type"
;;    <radius> <(centerpoint)> <direction>)
;;   ...
;;  )
;;  Direction flag is 1 for CCW, 2 for CW
;;  "type" is either "LINE" or "ARC"
;;  <> objects appear only if type is "ARC".
;;
;; Point is normal AutoLISP point list.
;; ------------------
;; Returns:
;;   1 if point is outside of datalist,
;;  -1 if point is inside datalist
;;-----------------------------------------------
(defun InOut (DataList TestPoint
              /         ;;local variable list
              Centroid  ;;point for inters test
              Num       ;;number of elements
              Pts       ;;just the points
              Tmp       ;;temporary usage
              VP1       ;;virtual point
              P1        ;;point list
              P2        ;;point list
              CP        ;;center point list
              RES       ;;temp result value
              I_Cnt     ;;intersection count
              Item      ;;foreach intersection
              TMPANG    ;;temp angle value
              PP3       ;;point list
              EA        ;;end angle
              SA        ;;start angle
              RD        ;;radius
              ALF      ;;temp angle
              )
   ;;
   ;;determine rough estimated centroid
   ;;by averaging all the start points.
   ;;
   (setq Num (length DataList)
         PTS (mapcar 'car DataList)
         Centroid ;;add min and max
           (mapcar '+ (find_lim 'min PTS)
                      (find_lim 'max PTS))
         Centroid
           (mapcar ;;then halve the values
             '(lambda (X)
                (/ X 2.0))
             Centroid)
         VP1 (* 3.0 ;;virtual point
                (distance
                   (find_lim 'min PTS)
                   (find_lim 'max PTS)
                )
             )
   ;;
   ;;vector from Virtual Point to Centroid is
   ;;intersection test vector
   ;;
         TMP (angle TestPoint Centroid)
         VP1 (polar TestPoint TMP VP1)
         I_Cnt nil ;;list of intersections
   )
   ;;
   ;;loop through datalist and find all
   ;;intersections, an even number of
   ;;intersections indicates that the ray is
   ;;starting outside of the object set, and
   ;;odd number means the ray begins inside
   ;;the object set.
   ;;
   (foreach Item DataList
      (cond
        ((= (caddr Item) "LINE")
            (if (setq TMP
                    (inters
                      TestPoint VP1
                      (car Item) (cadr Item)
                      'T))
                (setq I_Cnt (cons TMP I_Cnt)))
        )
        ((= (caddr Item) "ARC")
            (setq P1
                   (if (= (nth 5 Item) 1)
                       (car Item)
                       (cadr Item))
                  P2
                   (if (= (nth 5 Item) 1)
                       (cadr Item)
                       (car Item))
                  CP (nth 4 Item)
                  SA (angle CP P1)
                  EA (angle CP P2)
                  RD (nth 3 Item)
                  ALF (+
                        (/ PI 2.0)
                        (angle TestPoint VP1))
                  PP3 (inters
                         TestPoint VP1
                         CP (polar CP ALF 1.0)
                         nil)
            )
            (if (< EA SA)
               (setq EA (+ EA (* 2.0 PI))))
            (if PP3 (progn
               (setq TMPANG (angle CP PP3))
               (if (< TMPANG SA)
                  (setq TMPANG
                     (+ TMPANG (/ PI 2.0))))
               (if (and
                     (equal
                        (distance CP PP3)
                        RD
                        0.001)
                     (<= SA TMPANG EA))
                   (setq RES nil) ;;tangent
                   (setq RES
                      (line_arc
                         CP
                         (angle CP P1)
                         (angle CP P2)
                         (nth 3 Item)
                         TestPoint
                         VP1
                      )
                   )
               )
            ))
            (if RES (progn
               (setq I_Cnt
                  (cons (car Res) I_Cnt))
               (if (= (length RES) 2)
                  (setq I_Cnt
                     (cons (cadr Res) I_Cnt)))
            ))
        )
      )
   )
   ;;
   ;; remove duplicate intersections
   (if I_Cnt (progn
     (setq TMP (list (car I_Cnt)))
     (foreach Item (cdr I_Cnt)
      (if (apply 'and
            (mapcar
              '(lambda (X)
                (not (equal X Item 0.001)))
              TMP))
         (setq TMP (cons Item TMP)))
     )
     (setq I_Cnt TMP)
   ))
   (if (zerop (rem (length I_Cnt) 2)) 1 -1)
)
;;-----------------------------------------------
;;               LISTING 4
;;-----------------------------------------------
;; FIND_LIM
;; Given a list of points in WHO, function
;; applys the function in WHAT to the X,Y,Z data.
;; Returns a point list - can be used to find
;; the minimum, maximum, sum, and more.
;; Called by INOUT function in listing 3.
;;
(defun Find_Lim (What Who)
   (list (apply What (mapcar 'car Who))
         (apply What (mapcar 'cadr Who))
         (apply What (mapcar 'caddr Who))
   )
)

;; ==========================================================
;; ==========================================================

;
;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/hy ... raphics/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
;

[/pcode]

WHGF

判断是否在封闭多边形内只有这两种方法。射线法比角度法效率高。

fools

多谢snoopychen提供的资料,粗粗看了一下,发现了一些问题,说错的地方请指正:
1.Luis E的用的是line扫描,效果不如ray.
2.John Uhden用的虽然是ray,但只有一条,出错率很高,比如ray通过曲线端点或相切.我的程序也是ray,不过是旋转生成多条ray,然后统计.统计是程序的关键,可以有效规避上述错误,假设只生成7条ray,就可以允许其中三条是通过端点或相切,结果依然正确.
3.CAB收集的程序很长,没全看,好象是角度法,似乎不支持曲线.

就射线法和角度法,我的认识是:
1.射线法效率高,但出错率也高,特别对于自相交的曲线,很难判断.比如五角星,点击中间区域,不同方向的射线交点多数是2点,就好象在曲线外似的.
2.角度法效率低,但只要扫描精度够高,可以保证判断结果正确.角度法可以计算出曲线自相交时,点位于第几层.角度法的关键就是怎样在不出错的情况下提高速度.

我的程序没有附加考虑比如是否选到曲线,曲线是否闭合,UCS/WCS转换等错误判断,只提供了核心部分,这样阅读起来比较直观,实际使用的时候需要增加这些判断.

角度法有其自身的优势,看看下面的程序:
[pcode=lisp,true]
(DEFUN InCurve (pt ent / CURPA DXF0 DXF10 ELST ENDPA KEY NEAPT PTS)
  ;;三维空间曲线的判断未完善
  ;;(SETQ pt (TRANS pt 1 0))                ;To WCS
  (IF (EQUAL pt
             (SETQ NeaPt (VLAX-CURVE-GETCLOSESTPOINTTO ent pt))
             1E-6
      )
    (ALERT "***点在线上***")
    (PROGN (SETQ elst (ENTGET ent))
           (setq EndPa (VLAX-CURVE-GETENDPARAM ent))
           (SETQ dxf0 (CDR (ASSOC 0 elst)))           
           (COND ((EQUAL EndPa  (* 2 PI) 1E-6) ;CIRCLE+ELLIPSE
                  (setq dxf10 (cdr (assoc 10 elst)))
                  (IF (> (DISTANCE pt dxf10) (DISTANCE NeaPt dxf10))
                    (ALERT "***点在线外***")
                    (ALERT "***点在线内***")
                  )
                 )
                 ((= dxf0 "LWPOLYLINE")        ;支持bulge
                  (SETQ        pts (MAPCAR (FUNCTION CDR)
                                    (VL-REMOVE-IF (FUNCTION (LAMBDA (x) (/= 10 (CAR x)))) elst)
                            )
                  )                  
                  (IF (ZEROP (SETQ key (CheckInCurve pt ent (CONS NeaPt pts))))
                    (ALERT "***点在线外***")
                    (ALERT (STRCAT "***点在线内第 " (ITOA key) " 层***"))
                  )
                 )
                 ((= dxf0 "SPLINE")
                  (setq pts (list (vlax-curve-getStartPoint ent)))
                  (setq CurPa 0)
                  ;;可根据需要调整点数
                  (REPEAT 100
                    (SETQ pts (CONS (VLAX-CURVE-GETPOINTATPARAM
                                      ent
                                      (SETQ CurPa (+ CurPa (* EndPa 0.01)))
                                    )
                                    pts
                              )
                    )
                  )                  
                  (IF (ZEROP (SETQ key (CheckInCurve pt ent (CONS NeaPt pts))))
                    (ALERT "***点在线外***")
                    (ALERT (STRCAT "***点在线内第 " (ITOA key) " 层***"))
                  )
                 )
                 (T (ALERT "***暂时不支持的曲线类型***"))
           )
    )
  )
)
(DEFUN CheckInCurve (pt Curve pts)
  ;;添加最近点,也可以用其他方式,
  ;;当曲线存在重复点用vl-sort可能会出错,这里偷了点懒
  (SETQ        pts (VL-SORT pts
                     (FUNCTION (LAMBDA (p1 p2)
                                 (< (VLAX-CURVE-GETPARAMATPOINT Curve p1)
                                    (VLAX-CURVE-GETPARAMATPOINT Curve p2)
                                 )
                               )
                     )
            )
  )
  ;;狂刀兄写的,太精辟了,I Like it!!!
  (FIX
    (+ (/ (ABS
            (APPLY (FUNCTION +)
                   (MAPCAR (FUNCTION
                             (LAMBDA (p1 p2) (REM (- (ANGLE pt p1) (ANGLE pt p2)) PI))
                           )
                           (CONS (LAST pts) pts)
                           pts
                   )
            )
          )
          PI
       )
       0.5
    )
  )
)
[/pcode]

[pcode=lisp,true]
(DEFUN C:test (/ Curve Pt)
  (IF (SETQ Curve (CAR (ENTSEL "\n选择一条曲线:")))
    (WHILE (SETQ Pt (GETPOINT "\n点取测试点:"))
      (InCurve Pt Curve)
    )
  )
  (PRINC)
)
[/pcode]

雨箭风刀

射线法,只要选好了一个角度,一次就可以了
比如角度选 0.148945
想想看,出错的几率是多少?也许要100年才出现一次.你画的曲线(还必须是多义线或直线)刚好有一段角度是0.148945,而你选的点正好在这段曲线或其延伸线上.
不满意? 再来一条,(另外一个角度).
两次,足够了.
如果刚好....... 建议你去买彩票!

另,射线法的好处是支持所有曲线
前面的帖子没有细看,但应该有一个曲线闭合的判断,否则没有意义

如果再好好想一想,是有办法去规避那小小的出错几率的.关键是如何取到一个不可能和曲线边重合的角度.呵呵,话说到这个份上,应该不难做出来了

highflybird

谢谢fools和snoopychen提供的代码，我测试了fools的后面的那段代码，好像对于一般的情况还可以，然而对于自交叉严重的（就是交缠在在一起的，连人都分不出内外的）也不能判断准确，而且有个奇怪的现象，就是在同一区域的不同两点得出的判断也可能相反。

建议不做嵌套判断 ，我总觉得那样没多大实际意义，只做最外层（包括自交叉）的判断。

ScmTools

很好的讨论帖

cxjzxf

好像不带圆弧的才适用。

牢固

[pcode=lisp,true];;功能：判断点在任意封闭曲线内外，本函数自交曲线不适用 By Gu_xl 2012.07.31
;;返回: 点在封闭曲线上或曲线内，返回T，否则返回nil
;;测试: (gxl-PtInCurveP  (car(entsel "\n选择曲线:")) (getpoint))
(defun gxl-PtInCurveP  (POLY PT           /           CP           LW
                                   MINP           MAXP           MINX           MINY
                                   MAXX           MAXY           X           Y
                                   LST           CLOCKWISEP           ENDPARAM
                                   CURVELENGTH           PARAM   DIST
                                   D1           D2           DEV           )
  (cond
    ((equal pt
            (setq cp (vlax-curve-getclosestpointto poly pt))
            1e-8)) ;_ 点在曲线上 T
    ((progn
       (vla-GetBoundingBox
         (setq lw (vlax-ename->vla-object POLY))
         'MinP
         'MaxP)
       (setq MinP (vlax-safearray->list MinP))
       (setq MaxP (vlax-safearray->list MaxP))
       (setq minx (car MinP)
             miny (cadr MinP)
             maxx (car MaxP)
             maxy (cadr MaxP)
             x          (car pt)
             y          (cadr pt)
             )
       (or (< x minx)
           (> x maxx)
           (< y miny)
           (> y maxy)
           )
       )
     NIL ;_ 点在曲线最小包围盒外 nil
     )
    (t
     (setq
       lst (mapcar
             (function
               (lambda (x)
                 (vlax-curve-getParamAtPoint
                   lw
                   (vlax-curve-getClosestPointTo lw x)
                   )
                 )
               )
             (list minp
                   (list minx maxy)
                   MaxP
                   (list maxx miny)
                   )
             )
       ) ;_ 最小包围盒点在曲线上的投影点的参数表
     (setq ClockwiseP
            (if        (or
                  (<= (car lst) (cadr lst) (caddr lst) (cadddr lst))
                  (<= (cadr lst) (caddr lst) (cadddr lst) (car lst))
                  (<= (caddr lst) (cadddr lst) (car lst) (cadr lst))
                  (<= (cadddr lst) (car lst) (cadr lst) (caddr lst))
                  ) ;_  or
              t
              ) ;_  if
           ) ;_ 判断曲线是否为顺时针，顺时针 = T
     (setq endparam    (vlax-curve-getendparam poly)
           curvelength (vlax-curve-getDistAtParam poly endparam) ;_ 曲线长度
           )
     (setq param (vlax-curve-getparamatpoint poly cp)
           dist         (vlax-curve-getDistAtParam poly param)
           )
     (if (equal param (fix param) 1e-8)
       (progn
         (setq d1 (- dist 1e-8))
         (if (minusp d1)
           (setq d1 (+ curvelength d1))
           )
         (setq d2 (+ dist 1e-8))
         (if (> d2 curvelength)
           (setq d2 (- d2 curvelength)))
         (if (<        (distance pt (vlax-curve-getpointatdist poly d1))
                (distance pt (vlax-curve-getpointatdist poly d2))
                )
           (setq param (vlax-curve-getparamatdist poly d1))
           (setq param (vlax-curve-getparamatdist poly d2))
           )
         )
       )
     (setq dev (vlax-curve-getFirstDeriv poly param)
           cp  (vlax-curve-getpointatparam poly param)
           )
     (=        ClockwiseP
        (
         (lambda (p1 p2 p3)
           (<
             (* (- (car p2) (car p1)) (- (cadr p3) (cadr p1)))
             (* (- (cadr p2) (cadr p1)) (- (car p3) (car p1)))
             )
           )
          pt
          cp
          (mapcar '+ cp dev)
          )
        )
       )
    )
  )[/pcode]这个函数虽然对自交曲线不适用，但是效率比较高（一般在实际使用中需要对自交曲线进行内外判断的情况很少吧），这个算法比起射线法来说，效率要高不少！射线法虽然简单，可是因为要考虑极端情况，比如射线经过顶点、切点的情况，代码要复杂许多，导致效率低下！

梦断江南

发一个以前写的，有注解。修改一下可以作为函数
[pcode=lisp,true];|功能:测试点是否在曲线内------------------梁雄啸.2003.12
编程思路 :最近点构造xline 求点在交点表中的位置,偶数为在外,奇数在内;
v 1.0 完成;
|;
(defun c:ptin ()
  (princ "测试点是否在曲线内------------------梁雄啸.2003.12")
  (defun *error* (msg) (princ "\n出错! 未选中曲线")(setq *error* nil)) ;;出错处理;
  (vl-load-com)
  (setq pt    (getpoint "\n测试点:")
        objcur (vlax-ename->vla-object (car(entsel "\n选择封闭曲线:")))
        ptnea (vlax-curve-getclosestpointto objcur pt nil))  ;求最近点,不扩展;
  (cond
    ((equal pt ptnea) (princ "\n点在曲线上:"))
    (ptnea (setq ptnea (list (+ 0.019287 (car ptnea)) (cadr ptnea) (caddr ptnea)) ;使xline不水平或垂直,方便排序;
                 objxl   (vla-addxline mspace (vlax-3D-point pt) (vlax-3D-point ptnea)) ;构造xline,点需转换为ActiveX 兼容的（变体）三维点;
                 intlst  (x_intlst objxl objcur 0)      ;自定义函数:求两曲线交点表;
                 intlsti (vl-sort-i (cons pt intlst) '(lambda (x y) (< (car x) (car y))));pt加入表头,按x轴排序得序号表;
                 i       (vl-position 0 intlsti))       ;pt (i=0) 在排序的索引表中的位置;
           (entdel (entlast))                           ;删除xline;
           (cond ((= nil (vlax-curve-isClosed objcur)) (princ "\n曲线未封闭"))
                 ((= 0 (rem i 2))  (princ "\n点在曲线外"));奇偶判断,不能用(/ i 2);
                 (T (princ "\n点在曲线内")))
    )            
    ((= nil ptnea) (princ "\n未选中曲线"));此句无效,用出错处理;
  )
  (princ)
)[/pcode]

牢固

梦断江南 发表于 2013-5-5 01:47
发一个以前写的，有注解。修改一下可以作为函数
;|功能:测试点是否在曲线内------------------梁雄啸.2003 ...

这个还是对极端情况未处理，诸如xline穿过顶点，切线的情况！

GTJ116600

经典帖子，对初学者学习判断点和多义线（或多段线）非常有益。
再补充上老大提供的 http://www.xdcad.net/forum/forum.php?mod=redirect&goto=findpost&ptid=667514&pid=3446158帖子和http://bbs.xdcad.net/forum.php?mod=viewthread&tid=628915&page=1#pid3458036帖子的论文，可以深入细致的学习了

redcat2012

牢固 发表于 2013-5-5 08:08
这个还是对极端情况未处理，诸如xline穿过顶点，切线的情况！

判断曲线与射线各个交点是否为曲线的顶点或者切点，如果是将交点数量+1否则+0，最后累计数是否被整除判断

张和平

超级顶。。。。速度收藏

kkt123

超级顶。。。。速度收藏:(:(