是源码，就学习下，谢谢楼主。
算法绝对经典
支持分享、不错啊。。
学习中........什么时候能达到呢..
感谢分享，好好学习~
21楼的问题，我也同问。期待highflybird 兄有解决办法。
GP方法本身存在2个问题：
1、当第1步STEP1的数量不是足够大时，可能得不到正确结果；
2、当存在多个最大圆时得不到正确结果；
Saften方法存在2个问题：
1、自交曲线存在点在曲线内误判的可能。
2、运行效率低下。
以下GP模式改进版，支持自交曲线；（有待进一步改进得到最更多最大圆）
;;; maximum circle inscribed in a closed polyline
;;; Gian Paolo Cattaneo
(defun C:TesT (/POLY POLY_vl Dx Dy
LpList_vert_poly list_p_int
P_centerdist step1 step2 t1
t2
)
(prompt "\nSelect Polyline: ")
(if (setq POLY (ssname (ssget ":S" '((0 . "LWPOLYLINE"))) 0))
(progn
(setq i1
t1 (getvar "MilliSecs")
)
(setq step1 40) ;--> grid_1
(setq step2 20) ;--> grid_2
(setq list_vert_poly (LWPoly->List POLY 10))
(grid_1)
(Point_int)
(grid+)
(Point_center)
(repeat 2
(grid_2)
(Point_center)
)
(entmake
(list
(cons 0 "CIRCLE")
(cons 10 P_center)
(cons 40 dist)
)
)
(setq t2 (getvar "MilliSecs"))
(princ (strcat "time = " (rtos (- t2 t1) 2 0) " MilliSecs"))
(princ)
)
)
)
;; Restituisce una griglia di punti all'interno del getboundingbox della poly selezionata
;; Returns a grid of points within the BoundingBox of the selected poly
(defun grid_1 (/ p1 p2 X1 Y1 l1)
(vla-getboundingbox (vlax-ename->vla-object POLY) 'p1 'p2)
(setq p1 (vlax-safearray->list p1)
p2 (vlax-safearray->list p2)
p1 (list (car p1) (cadr p1))
p2 (list (car p2) (cadr p2))
)
(setq Dx (/ (- (car p2) (car p1)) step1))
(setq Dy (/ (- (cadr p2) (cadr p1)) step1))
(setq Lp (list p1)
X1 (car p1)
Y1 (cadr p1)
)
(repeat step1
(setq Lp (cons (list (setq X1 (+ X1 Dx)) Y1) Lp))
)
(setq Lp (list Lp))
(repeat step1
(setqLp (cons (mapcar (function (lambda (x)
(list (car x) (+ (cadr x) Dy))
)
)
(car lp)
)Lp)
)
)
(setq Lp (apply (function append) Lp))
)
;; Restituisce una griglia di punti intorno al punto centrale (provvisorio)
;; Returns a grid of points around the center point (provisional)
(defun grid_2 (/ P1_ P> n)
(setq list_p_int nil)
(setq P1_ (list (- (car P_center) (* Dx 2))
(- (cadr P_center) (* Dy 2))
)
)
(setq Dx (/ (* 4 Dx) step2))
(setq Dy (/ (* 4 Dy) step2))
(setq n 0)
(setq P> P1_)
(setq list_p_int (list P1_))
(repeat (* (1+ step2) step2)
(setq P> (list (+ (car P>) Dx) (cadr P>)))
(setq list_p_int (cons P> list_p_int))
(setq n (1+ n))
(if (= n step2)
(progn
(setq n 0)
(setq P1_ (list (car P1_) (+ (cadr P1_) Dy)))
(setq P> P1_)
(setq list_p_int (cons P> list_p_int))
)
)
)
)
;; restituisce la lista dei punti interni ad un poligono
;; dati:- lista coordinate dei punti -> Lp
;; - lista coordinate vertici poligono -> list_vert_poly
;; Returns the list of points inside the polyline
(defun Point_int (/ n Pr cont attr p# Pa Pa_ Pb)
(setq list_p_int nil)
(foreach Pr Lp
(if (> (Point-in-ClosedCurve-p list_vert_poly Pr) 0)
(setq list_p_int (cons Pr list_p_int))
)
)
)
;; Infittisce la griglia inserendo altri punti
;; nel centro delle diagonali tra i punti interni
;; Insert points (interior) to increase the density of the grid
(defun grid+ (/ G+)
(setq G+
(mapcar '(lambda (x)
(list (+ (car x) (/ Dx 2)) (+ (cadr x) (/ Dy 2)))
)
list_p_int
)
)
(setq list_p_int (append G+ list_p_int))
)
;; Da una lista di punti restituisce quello più lontano da un oggetto
;; dati:- lista dei punti -> list_p_int
;; - oggetto -> POLY_vl
;; Returns the farthest point from the polyline
(defun Point_center (/ Pa n Pvic)
(setq Dist 1e-7)
(setq P_center nil)
(foreach Pa list_p_int
(setq Pvic (vlax-curve-getClosestPointTo Poly Pa))
(if (> (distance Pa Pvic) Dist)
(progn
(setq P_center Pa)
(setq Dist (distance Pa Pvic))
)
)
)
)
;;
(defun LWPoly->List (en acc / a b vetex bu p1 p2 l r ang an1 N)
;;Acc --- 0 ~ 99
(setq ent (entget en))
(while (setq ent (member (assoc 10 ent) ent))
(setq b (cons (cdar ent) b)
ent (member (assoc 42 ent) ent)
b (cons (cdar ent) b)
ent (cdr ent)
vetex (cons b vetex)
b nil
)
)
(while vetex
(setq a (car vetex)
vetex (cdr vetex)
bu (car a)
p1 (cadr a)
)
(if l
(setq p2 (car l))
(setq p2 (cadr (last vetex))
l(cons p2 l)
)
)
(if (equal bu 0 1e-6)
(setq l (cons p1 l))
(progn
(setq ang (* 2 (atan bu))
r (/ (distance p1 p2)
(* 2 (sin ang))
)
c (polar p1
(+ (angle p1 p2) (- (/ pi 2) ang))
r
)
r (abs r)
ang (abs (* ang 2.0))
N (abs (fix (/ ang 0.0174532925199433)))
N (min N (1+ Acc))
)
(if (= N 0)
(setq l (cons p1 l))
(progn
(setq an1 (/ ang N)
ang (angle c p2)
)
(if (not (minusp bu))
(setq an1 (- an1))
)
(repeat (1- N)
(setq ang (+ ang an1))
(setq l (cons (polar c ang r) l))
)
(setq l (cons p1 l))
)
)
)
)
)
l
)
;;
;; This method suggest by Lee Mac from http://en.wikipedia.org/wiki/Winding_number
;; function : determin the point position with the closed curve by widding-number method
;; l---- point set of a Closed Curve , First item must same as Last item .
;; pt ---- a given point to determin position with the Closed Curve
;;; return a num
;;; -----1pt out of curve
;;; ---- 0pt at curve
;;; ---- 1pt in curve
;; by GSLS(SS) 2012-08-02
(defun Point-in-ClosedCurve-p (l pt / ang p1 p2 n r at)
(setq ang 0.0
atnil
)
(while (and (cadr l) (not at))
(setq p1 (car l)
p2 (cadr l)
l(cdr l)
)
(if (equal p1 p2 1e-6)
(setq an1 0.0)
(setq an1
((lambda (/ a b c d e f g)
(setq b (distance p1 pt)
c (distance p2 pt)
a (distance p1 p2)
d (- (* (- (car p1) (car pt)) (- (cadr p2) (cadr pt)))
(* (- (car p2) (car pt)) (- (cadr p1) (cadr pt)))
)
)
(if (and (equal d 0.0 1e-4) (setq at T))
pi
(progn
(setq
e (+ (* b b) (* c c) (* -1 a a))
f (abs ((lambda (x)
(cond ((equal x 0.0 1e-6)(* pi 0.5))
((equal x 1.0 1e-6)0.0)
((atan (/ (sqrt (- 1 (* x x)))
x
)
))
))
(/ e 2. b c)
)
)
g (if (> d 0)1-1)
)
(if (< e 0)
(* g (- pi f))
(* g f)
)
)
)
)
)
)
)
(setq ang (+ ang an1))
)
;;deal widding number
(if at
0
(progn
(setq n (/ ang 2. pi))
(if (equal (fix n) n 1e-4)
(setq n (fix n))
(if (and (> n 0) (equal (1+ (fix n)) n 1e-4))
(setq n (1+ (fix n)))
(if (and (< n 0) (equal (1- (fix n)) n 1e-4))
(setq n (1- (fix n)))
(setq n (fix n))
)
)
)
(if (= (rem n 2) 0)
-1
1
)
)
)
)
;|
(defun c:t1 ( / en l n)
(setq en (car (entsel))
l (LWPoly->List en 10))
(while (setq pt (getpoint ))
(setq n (Point-in-ClosedCurve-p l pt))
(cond ((> n 0)
(alert "IN"))
((= n 0)
(alert "AT"))
(t
(alert "OUT")))))
|;
对于椭圆及SPLine可以用下面函数取点:;; get point set of a closed curve by order
;; this function you improve by yourself acordding your need .
(defun get-closed-curve-pts (en / ent et)
;;by GSLS(SS)
(setq
ent (entget en)
et (cdr (assoc 0 ent))
)
(cond
((= et "LWPOLYLINE")
((lambda (/ a b vetex bu p1 p2 l r ang an1 N)
(while (setq ent (member (assoc 10 ent) ent))
(setq b (cons (cdar ent) b)
ent (member (assoc 42 ent) ent)
b (cons (cdar ent) b)
ent (cdr ent)
vetex (cons b vetex)
b nil
)
)
(while vetex
(setq a (car vetex)
vetex (cdr vetex)
bu (car a)
p1 (cadr a)
)
(if l
(setq p2 (car l))
(setq p2 (cadr (last vetex))
l(cons p2 l)
)
)
(if (equal bu 0 1e-6)
(setq l (cons p1 l))
(progn
(setq ang (* 2 (atan bu))
r (/ (distance p1 p2)
(* 2 (sin ang))
)
c (polar p1
(+ (angle p1 p2) (- (/ pi 2) ang))
r
)
r (abs r)
ang (abs (* ang 2.0))
N (abs (fix (/ ang 0.0174532925199433)))
)
(if (= N 0)
(setq l (cons p1 l))
(progn
(setq an1 (/ ang N)
ang (angle c p2)
)
(if (not (minusp bu))
(setq an1 (- an1))
)
(repeat (1- N)
(setq ang (+ ang an1))
(setq l (cons (polar c ang r) l))
)
(setq l (cons p1 l))
)
)
)
)
)
l
)
)
)
((= et "CIRCLE")
((lambda (c R / sa ptl)
(setq sa 0.0)
(repeat 180
(setq ptl (cons (polar c sa R) ptl)
sa(+ sa 0.0174532925199433)
)
)
(setq ptl (reverse ptl))
(append
ptl
(mapcar (function
(lambda (p)
(mapcar (function +) c (mapcar (function -) c p))
)
)
ptl
)
)
)
(cdr (assoc 10 ent))
(cdr (assoc 40 ent))
)
)
((= et "SPLINE")
((lambda (/ r l _oce)
(setq _oce (getvar "CMDECHO"))
(setvar "CMDECHO" 0)
(if (vl-catch-all-apply
(function vl-cmdf)
(list "_PEDIT"
(vlax-vla-object->ename
(vla-copy (vlax-ename->vla-object en))
)
""
10
""
)
)
(progn
(setq l (ss-assoc 10 (entget (setq r (entlast)))))
(if (vlax-curve-isClosed r)
(setq l (append l (list (car l))))
)
(entdel r)
)
)
(setvar "CMDECHO" _oce)
(append l (list (car l)))
)
)
)
((= et "ELLIPSE")
((lambda (/ e l _os)
(setq _os (getvar "OSMODE"))
(setvar "OSMODE" 0)
(vl-catch-all-apply
(function vla-offset)
(list (vlax-ename->vla-object en) 0.1)
)
(setq e (entlast))
(vl-catch-all-apply
(function vla-offset)
(list (vlax-ename->vla-object (entlast)) -0.1)
)
(entdel e)
(setq e (entlast))
(setq l (ss-assoc 10 (entget e)))
(entdel e)
(setvar "OSMODE" _os)
(append l (list (car l)))
)
)
)
)
)
GP再改进：
1.为了适应任意管状封闭曲线, 将第一步网格数的确定方式改为自适应方式，即使用周长与面积的开根号的比例来控制，这样可以避免陷入非最优解的困扰；
2.强制第2步Step2的网格细分数为10，10等分已经足够让半径的精度提高1个等级了；
3.改进Grid_2的计算区域，减少为原来的1/4，即将原来为上次网格划分的16个区格改为4个区格；
;;old grid_2 area
+ + + + +
+ + + + +
+ + o + +
+ + + + +
+ + + + +
;;new grid_2 area
+ + +
+ o +
+ + +4.改进点在曲线点集内的计算函数，仅使用angle函数，效率比原来提高了非常多倍；
5.使用控制精度的while循环来替代原本由次数控制的repeat循环，既提高了精度又可能减少大量运算，一般情况下都比原来快些；
6.限制循环次数，避免存在条状最大圆的情况陷入死循环。
新代码如下：
;;; maximum circle inscribed in a closed polyline
;;; writed by Gian Paolo Cattaneo
;;; edited by GSLS(SS) 2012-8-5
(defun C:TesT (/ POLY POLY_vl Dx Dy
Lp List_vert_poly list_p_int
P_center dist step1 step2 t1
t2R0area len i
)
(gc)
(prompt "\nSelect Polyline: ")
(if (setq POLY (ssname (ssget ":S" '((0 . "LWPOLYLINE"))) 0))
(progn
(setq i1
t1 (getvar "MilliSecs")
)
(setq area (vlax-curve-getArea poly)
len (vlax-curve-getDistAtParam poly (vlax-curve-getEndParam poly)))
(setq step1 (max 10 (fix (/ len 0.4 (sqrt area)))));_--> grid_1
(setq step2 10);_--> grid_2
(setq list_vert_poly (LWPoly->List POLY 10))
(grid_1)
(Point_int)
(grid+)
(Point_center)
(setq i 0)
(while (and (> (- Dist R0) 1e-8)(< i 10))
(grid_2)
(Point_center)
(setq i (1+ i))
)
(entmake
(list
(cons 0 "CIRCLE")
(cons 10 P_center)
(cons 40 dist)
)
)
(setq t2 (getvar "MilliSecs"))
(princ (strcat "time = " (rtos (- t2 t1) 2 0) " MilliSecs"))
(princ)
)
)
)
;; Restituisce una griglia di punti all'interno del getboundingbox della poly selezionata
;; Returns a grid of points within the BoundingBox of the selected poly
(defun grid_1 (/ p1 p2 X1 Y1 l1)
(vla-getboundingbox (vlax-ename->vla-object POLY) 'p1 'p2)
(setq p1 (vlax-safearray->list p1)
p2 (vlax-safearray->list p2)
p1 (list (car p1) (cadr p1))
p2 (list (car p2) (cadr p2))
)
(setq Dx (/ (- (car p2) (car p1)) step1))
(setq Dy (/ (- (cadr p2) (cadr p1)) step1))
(setq Lp (list p1)
X1 (car p1)
Y1 (cadr p1)
)
(repeat step1
(setq Lp (cons (list (setq X1 (+ X1 Dx)) Y1) Lp))
)
(setq Lp (list Lp))
(repeat step1
(setq Lp (cons (mapcar (function (lambda (x)
(list (car x) (+ (cadr x) Dy))
)
)
(car lp)
)
Lp
)
)
)
(setq Lp (apply (function append) Lp))
)
;; Restituisce una griglia di punti intorno al punto centrale (provvisorio)
;; Returns a grid of points around the center point (provisional)
(defun grid_2 (/ X1 Y1 P1)
(setq list_p_int nil
X1 (- (car P_center) Dx)
Y1 (- (cadr P_center) Dy)
P1 (list X1 Y1)
Dx (/ (* 2 Dx) step2)
Dy (/ (* 2 Dy) step2)
)
(setq list_p_int (list P1))
(repeat step2
(setq list_p_int (cons (list (setq X1 (+ X1 Dx)) Y1) list_p_int))
)
(setq list_p_int (list list_p_int))
(repeat step2
(setq list_p_int
(cons (mapcar (function (lambda (x)
(list (car x) (+ (cadr x) Dy))
)
)
(car list_p_int)
)
list_p_int
)
)
)
(setq list_p_int (apply (function append) list_p_int))
)
;; restituisce la lista dei punti interni ad un poligono
;; dati:- lista coordinate dei punti -> Lp
;; - lista coordinate vertici poligono -> list_vert_poly
;; Returns the list of points inside the polyline
(defun Point_int (/ n Pr cont attr p# Pa Pa_ Pb)
(setq list_p_int nil)
(foreach Pr Lp
(if (> (Point-in-ClosedCurve-p list_vert_poly Pr) 0)
(setq list_p_int (cons Pr list_p_int))
)
)
)
;; Infittisce la griglia inserendo altri punti
;; nel centro delle diagonali tra i punti interni
;; Insert points (interior) to increase the density of the grid
(defun grid+ (/ G+)
(setq G+
(mapcar '(lambda (x)
(list (+ (car x) (/ Dx 2)) (+ (cadr x) (/ Dy 2)))
)
list_p_int
)
)
(setq list_p_int (append G+ list_p_int))
)
;; Da una lista di punti restituisce quello più lontano da un oggetto
;; dati:- lista dei punti -> list_p_int
;; - oggetto -> POLY_vl
;; Returns the farthest point from the polyline
(defun Point_center (/ Pa Pvic)
(foreach Pa list_p_int
(setq Pvic (vlax-curve-getClosestPointTo Poly Pa))
(if (> (distance Pa Pvic) Dist)
(setq P_center Pa
R0 Dist
Dist (distance Pa Pvic))
)
)
)
;;
(defun LWPoly->List (en acc / a b vetex bu p1 p2 l r ang an1 N)
;;Acc --- 0 ~ 99
(setq ent (entget en))
(while (setq ent (member (assoc 10 ent) ent))
(setq b (cons (cdar ent) b)
ent (member (assoc 42 ent) ent)
b (cons (cdar ent) b)
ent (cdr ent)
vetex (cons b vetex)
b nil
)
)
(while vetex
(setq a (car vetex)
vetex (cdr vetex)
bu (car a)
p1 (cadr a)
)
(if l
(setq p2 (car l))
(setq p2 (cadr (last vetex))
l(cons p2 l)
)
)
(if (equal bu 0 1e-6)
(setq l (cons p1 l))
(progn
(setq ang (* 2 (atan bu))
r (/ (distance p1 p2)
(* 2 (sin ang))
)
c (polar p1
(+ (angle p1 p2) (- (/ pi 2) ang))
r
)
r (abs r)
ang (abs (* ang 2.0))
N (abs (fix (/ ang 0.0174532925199433)))
N (min N (1+ Acc))
)
(if (= N 0)
(setq l (cons p1 l))
(progn
(setq an1 (/ ang N)
ang (angle c p2)
)
(if (not (minusp bu))
(setq an1 (- an1))
)
(repeat (1- N)
(setq ang (+ ang an1)
l (cons (polar c ang r) l)
)
)
(setq l (cons p1 l))
)
)
)
)
)
l
)
;;
;; This method suggest by Lee Mac from http://en.wikipedia.org/wiki/Winding_number
;; function : determin the point position with the closed curve by widding-number method
;; l---- point set of a Closed Curve , First item must same as Last item .
;; pt ---- a given point to determin position with the Closed Curve
;; return a num
;; -----1pt out of curve
;; ---- 0pt at curve
;; ---- 1pt in curve
;; by GSLS(SS) 2012-08-02
;; Edited : Change widding number's Acc 1e-4 into 1e-2 ,
;; for checking point in or out it's a integer
;; Edited 2012-8-5
;; Improved vector angle calculation , only use angle function .
;;
(defun Point-in-ClosedCurve-p (l pt / ang p1 p2 n r at)
(setq ang 0.0)
(while (and (cadr l) (not at))
(setq p1(car l)
p2(cadr l)
l (cdr l)
an1 (- (angle pt p2) (angle pt p1))
)
(if (< an1 (- pi))
(setq an1 (+ an1 pi pi))
(if (> an1 pi)
(setq an1 (- an1 pi pi))
)
)
(if (equal (abs an1) pi 1e-14);_If it's just used to solve the maximum radius of the circle,
;_here'sprecision 1e-14 can be set lower , such as 1e-1 ,
;_this will exclude the point of the curve edge .
;_but for ultra-narrow four-point rectangle will generate an error .
(setq at T)
)
(setq ang (+ ang an1))
)
;;deal widding number
(if at
0
(progn
(setq n (/ ang 2. pi))
(if (and (> n 0) (equal (1+ (fix n)) n 1e-2))
(setq n (1+ (fix n)))
(if (and (< n 0) (equal (1- (fix n)) n 1e-2))
(setq n (1- (fix n)))
(setq n (fix n))
)
)
(if (= (rem n 2) 0)
-1
1
)
)
)
)
收藏了，下来学习