【飞鸟集】最小包围圆的最佳算法

highflybir · 发表于 2006-11-11 23:43

这类话题本想在lisp开发版块中发布的，但因为与几何算法密切相关，故贴到此处。

最小包围圆的算法在实际和理论中都有值得探讨的必要。
在国内网站，对于此类算法鲜有介绍，今天完成了它的一个lisp程序，甚高兴。
《计算几何-算法与应用》中介绍的方法为随机增量式算法,可在O(n)的期望时间中算出来，
而这个算法有别于上种算法，其时间为O(|lg(d/R)|*n)，也就是说很可能比上种算法时间更少。
请大家指教指教检查。附件为其lisp程序，加载，运行test然后选取点对象即可。

;;;************************************
;;;求最小包围圆的lisp程序--------------
;;;其算法为参见了有关文献--------------
;;;这种算法在退化很严重的情况结果也正确
;;;其中程序主段是核心算法，其他的附加程
;;;序为取点，画点，画圆和半径，用来测试
;;;************************************
(defun C:test (/ olderr en errmsg oce
oldmodessp ptlist x cen radius ptmax)
;;定义错误函数和预处理---------------
(setvar "errno" 0)
(setq olderr *error*)
(defun *error* (msg / en errmsg)
    (setq en (getvar "errno"))
    (setq errmsg (strcat "errno=" (itoa en) "\nError:" msg))
    (alert errmsg)
    (setq *error* olderr)
)
(graphscr)
(setq oldmode (getvar "osmode"))
(setq oce (getvar "cmdecho"))
(setvar "cmdecho" 0)
(command ".ucs" "W")
;;取点，画点，并对函数用时计算-------
(setq ssp (ssget '((0 . "POINT"))))
(setq ptlist (ssgetpoint ssp))
(setq t1 (getvar "CDATE"))
(setq x (mincir ptlist))
(setq t2 (getvar "CDATE"))
(setq usetime (* (- t2 t1) 1e6))
(princ (strcat "\n用时=" (rtos usetime 2 6) "秒"))
(if (= nil x)
    (progn
      (alert "点的有效数目太小，请重新输入!")
      (command ".ucs" "p")
      (setvar "osmode" oldmode)
      (setvar "cmdecho" oce)
      (princ "\n")
      (princ)
    )
    (progn
      (setq cen (car x) radius (cadr x) ptmax (caddr x))
      ;;;画圆及半径，列出圆的圆心半径值
      (entmake
        (append
          '((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle"))
          (list (cons 10 cen))(list (cons 40 radius))(list (cons 62 1))
        )
      )
      (entmake
        (append
          '((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine"))
          (list (cons 10 cen))(list (cons 11 ptmax))(list (cons 62 1))
        )
      )
      (command ".ucs" "p")
      (setvar "osmode" oldmode)
      (setvar "cmdecho" oce)
      (princ "\n")
      (list cen radius)
    )
)
)
;;;************************************
;;;求最小包围圆的函数，空集返回空集，否
;;;则返回最小圆的圆心，半径和圆上的一点
;;;这是程序的主段----------------------
;;;************************************
(defun mincir (ptlist / p1 p2 p3 ptmax cen_r cen radius)
;;定义中点函数，本来R2004版中无须定义
;;但不知道为什么到R2006版没有定义了。
(defun mid (p1 p2)
    (if (or nil (= (length p1) 2) (= (length p2) 2))
      (list (/ (+ (car p1) (car p2)) 2.0) (/ (+ (cadr p1) (cadr p2)) 2.0) 0.0)
      (list (/ (+ (car p1) (car p2)) 2.0) (/ (+ (cadr p1) (cadr p2)) 2.0) (/ (+ (caddr p1) (caddr p2)) 2.0))
    )
)
;;判断有效点个数---------------------
(cond
    ((= (length ptlist) 0)
     nil
    )
    ((= (length ptlist) 1)
     (progn
       (alert "点集合为一点,最小圆半径为0")
       (list (car ptlist) 0 (car ptlist))
     )
    )
    ((= (length ptlist) 2)
     (progn
       (alert "点集合为两点，最小圆直径为其两点距离，\n圆心为其连线中点")
       (setq cen (mid (car ptlist) (cadr ptlist)) radius (/ (distance (car ptlist) (cadr ptlist)) 2))
       (list cen radius (car ptlist))
     )
    )
    (t
     (progn
       ;;上面啰嗦的一大段在实际情况中一般不会出现
       ;;判断点是否在圆内------------------------
       (defun in1 (pt cen r)
(if (> (- r (distance pt cen)) (- 1e-8))
    t
    nil
)
       )
       ;;判断点集是否在圆内----------------------
       (defun in2 (ptl cen r)
(if (apply 'and (mapcar '(lambda (x) (in1 x cen r)) ptl))
    t
    nil
)
       )
       ;;定义三点最小圆圆心及其半径，若是锐角三角
       ;;形，则是其三点圆，否则是其最大边的直径圆
       (defun 3pc (pa pb pc / a b c l p ja jb jc ppa ppb ppc cen radius)
(setq a (list (distance pb pc) pa))
(setq b (list (distance pc pa) pb))
(setq c (list (distance pa pb) pc))
(setq l (list a b c))
(setq p (/ (+ (car a) (car b) (car c)) 2))
(setq a (nth (car (vl-sort-i (mapcar 'car l) '>)) l))
(setq b (nth (cadr (vl-sort-i (mapcar 'car l) '>)) l))
(setq c (nth (caddr (vl-sort-i (mapcar 'car l) '>)) l))
(setq l (+ (* (car b) (car b)) (* (car c) (car c)) (* (car a) (car a) -1.0)))
;;上面l利用了余弦定理作为判断-----------
(if (< l 1e-8)
    (list (mid (cadr b) (cadr c))(/ (car a) 2)(list (cadr b) (cadr c) (cadr a)))
    (progn
      (setq ja (angle pb pc))
      (setq jb (angle pc pa))
      (setq jc (angle pa pb))
      (setq ppc (polar (mid pa pb) (+ (/ pi 2) jc) p))
      (setq ppa (polar (mid pb pc) (+ (/ pi 2) ja) p))
      (setq ppb (polar (mid pc pa) (+ (/ pi 2) jb) p))
      (setq cen (inters ppc (mid pa pb) ppa (mid pb pc) nil))
      (setq radius (distance cen pa))
      (list cen radius (list pa pb pc))
    )
)
       )
       ;;定义四点的最小圆圆心半径，并返回三点坐标
       (defun 4pc (p1 p2 p3 ptmax / pts 3pt)
(setq pts (list (3pc p1 p2 p3) (3pc p1 p2 ptmax) (3pc p1 p3 ptmax) (3pc p2 p3 ptmax)))
(setq 3pt (vl-sort-i (mapcar 'cadr pts) '<))
(setq pts (list (nth (car 3pt) pts) (nth (cadr 3pt) pts)
           (nth (caddr 3pt) pts)(nth (cadddr 3pt) pts)))
(nth (vl-position t (mapcar '(lambda (x) (in2 (list p1 p2 p3 ptmax) (car x) (cadr x))) pts)) pts)
       )
       ;;定义求点集中离圆心最远的点的函数--------
       (defun maxd-cir (ptl cen / distl)
(setq distl (mapcar '(lambda (x) (distance x cen)) ptl))
(nth (car (vl-sort-i distl '>)) ptl)
       )
       ;;开始递归运算----------------------------
       (setq p1 (car ptlist) p2 (cadr ptlist) p3 (caddr ptlist))
       (setq cen_r (3pc p1 p2 p3))
       (setq ptmax (maxd-cir ptlist (car cen_r)))
       (while (= nil (in1 ptmax (car cen_r) (cadr cen_r)))
         (setq cen_r (4pc p1 p2 p3 ptmax))
         (setq p1 (car (caddr cen_r)) p2 (cadr (caddr cen_r)) p3 (caddr (caddr cen_r)))
(setq ptmax (maxd-cir ptlist (car cen_r)))
       )
       (list (car cen_r) (cadr cen_r) ptmax)
      );;for progn
    );; for t
);;    for cond
);;      for defun
;;以下代码来自晓东
;;定义取点函数----
(defun ssgetpoint (ss / i listpp a b c)
(setq i 0 listpp nil )
(if ss
    (repeat (sslength ss)
      (setq a (ssname ss i))
      (setq b (entget a))
      (setq c (cdr (assoc 10 b)))
      (setq listpp (cons c listpp))
      (setq i (1+ i))
    )
)
listpp
)

highflybird · 发表于 2010-7-13 11:20

改进了一下，速度提高不少

;;;************************************
;;;求最小包围圆的lisp程序--------------
;;;其算法为参见了有关文献--------------
;;;这种算法在退化很严重的情况结果也正确
;;;其中程序主段是核心算法，其他的附加程
;;;序为取点，画点，画圆和半径，用来测试
;;;************************************
(defun C:test (/ CEN PTLIST PTMAX RADIUS SL SS T0 X)
;;取点，画点，并对函数用时计算-------
(setq sl '((0 . "POINT,LINE,POLYLINE,LWPOLYLINE")))
(setq ss (ssget sl))
(setq ptlist (ssgetpoint ss))
(setq t0 (getvar "TDUSRTIMER"))
(setq x (mincir ptlist))
(princ "\n用时")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400)) ;结束计时
(princ "秒")
(if (null x)
(alert "点的有效数目太小，请重新输入!")
(progn
(setq cen (car x)
radius (cadr x)
ptmax (caddr x)
)
;;画圆及半径，列出圆的圆心半径值
(entmake
(append
'((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle"))
(list (cons 10 cen))
(list (cons 40 radius))
(list (cons 62 1))
)
)
(entmake
(append
'((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine"))
(list (cons 10 cen))
(list (cons 11 ptmax))
(list (cons 62 1))
)
)
(list cen radius)
)
)
)
;;;************************************
;;;求最小包围圆的函数，空集返回空集，否
;;;则返回最小圆的圆心，半径和圆上的一点
;;;这是程序的主段----------------------
;;;************************************
(defun mincir (ptlist / CEN CEN_R P1 P2 P3 PTMAX R RADIUS X i)
;;判断有效点个数---------------------
(cond
((= (length ptlist) 0)
nil
)
((= (length ptlist) 1)
(alert "点集为一点,最小圆半径为0")
(list (car ptlist) 0 (car ptlist))
)
((= (length ptlist) 2)
(alert "点集为两点，最小圆为过两点的圆")
(setq cen (mid (car ptlist) (cadr ptlist))
radius (/ (distance (car ptlist) (cadr ptlist)) 2)
)
(list cen radius (car ptlist))
)
(t
;;开始递归运算----------------------------
(setq p1 (car ptlist)
p2 (cadr ptlist)
p3 (caddr ptlist)
)
(setq cen_r (3pc p1 p2 p3))
(setq ptmax (maxd-cir ptlist (car cen_r)))
(setq i 0)
(while (null (in1 ptmax (car cen_r) (cadr cen_r)))
(setq cen_r (4pc p1 p2 p3 ptmax))
(setq p1 (car (caddr cen_r))
p2 (cadr (caddr cen_r))
p3 (caddr (caddr cen_r))
)
(setq ptmax (maxd-cir ptlist (car cen_r)))
(setq i (1+ i))
)
(list (car cen_r) (cadr cen_r) ptmax)
)
)
)
(defun make-line (p q)
(entmake
(list
'(0 . "LINE")
(cons 10 p)
(cons 11 q)
)
)
)
;;以下代码来自晓东
;;定义取点函数----
(defun ssgetpoint (ss / i l a b c)
(setq i 0)
(if ss
(repeat (sslength ss)
(setq a (ssname ss i))
(setq i (1+ i))
(setq b (entget a))
(setq c (cdr (assoc 10 b)))
(setq l (cons c l))
)
)
(reverse l)
)
(defun mid (p1 p2)
(list
(* (+ (car p1) (car p2)) 0.5)
(* (+ (cadr p1) (cadr p2)) 0.5)
(* (+ (caddr p1) (caddr p2)) 0.5)
)
)
;;判断点是否在圆内------------------------
(defun in1 (pt cen r)
(< (- (distance pt cen) r) 1e-8)
)
;;判断点集是否在圆内----------------------
(defun in2 (ptl cen r / pts pt)
(setq pts ptl)
(while (and (setq pt (car pts))
(in1 pt cen r)
)
(setq pts (cdr pts))
)
(null pt)
)
;;定义三点最小圆圆心及其半径，若是锐角三角
;;形，则是其三点圆，否则是其最大边的直径圆
(defun 3pc (pa pb pc / D MIDPT)
(cond
((in1 pc (setq midpt (mid pa pb)) (setq d (/ (distance pa pb) 2)))
(list midpt d (list pa pb pc))
)
((in1 pa (setq midpt (mid pb pc)) (setq d (/ (distance pb pc) 2)))
(list midpt d (list pb pc pa))
)
((in1 pb (setq midpt (mid pc pa)) (setq d (/ (distance pc pa) 2)))
(list midpt d (list pc pa pb))
)
(t
(3pcircle pa pb pc)
)
)
)
;;; 三点圆函数
(defun 3PCirCle (P0 P1 P2 / X0 Y0 X1 Y1 X2 Y2 DX1 DY1 DX2 DY2 D 2D C1 C2 CE)
(setq X0 (car P0)
Y0 (cadr P0)
X1 (car P1)
Y1 (cadr P1)
X2 (car P2)
Y2 (cadr P2)
DX1 (- X1 X0)
DY1 (- Y1 Y0)
DX2 (- X2 X0)
DY2 (- Y2 Y0)
)
(setq D (- (* DX1 DY2) (* DX2 DY1)))
(if (/= D 0.0)
(progn
(setq 2D (+ D D)
C1 (+ (* DX1 (+ X0 X1)) (* DY1 (+ Y0 Y1)))
C2 (+ (* DX2 (+ X0 X2)) (* DY2 (+ Y0 Y2)))
CE (List (/ (- (* C1 DY2) (* C2 DY1)) 2D)
(/ (- (* C2 DX1) (* C1 DX2)) 2D)
)
)
(list CE (distance CE P0) (list p0 p1 p2))
)
)
)
;;定义四点的最小圆圆心半径，并返回三点坐标
(defun 4pc (p1 p2 p3 ptmax / pts mind minr r 4ps)
(setq pts (list (3pc p1 p2 ptmax)
(3pc p1 p3 ptmax)
(3pc p2 p3 ptmax)
)
)
(setq 4ps (list p1 p2 p3 ptmax))
(setq minr 1e308)
(foreach n pts
(setq r (cadr n))
(if (and (< r minr)
(in2 4ps (car n) r)
)
(setq mind n)
)
)
mind
)
;;定义求点集中离圆心最远的点的函数--------
(defun maxd-cir (ptl cen / pmax dmax d)
(setq dmax 0.0)
(foreach pt ptl
(if (> (setq d (distance pt cen)) dmax)
(setq dmax d
pmax pt
)
)
)
pmax
)

;;;************************************
;;;求最小包围圆的lisp程序--------------
;;;其算法为参见了有关文献--------------
;;;这种算法在退化很严重的情况结果也正确
;;;其中程序主段是核心算法，其他的附加程
;;;序为取点，画点，画圆和半径，用来测试
;;;************************************
(defun C:test (/ CEN PTLIST PTMAX RADIUS SL SS T0 X)
  ;;取点，画点，并对函数用时计算-------
  (setq sl '((0 . "POINT,LINE,POLYLINE,LWPOLYLINE")))
  (setq ss (ssget sl))
  (setq ptlist (ssgetpoint ss))
  (setq t0 (getvar "TDUSRTIMER"))
  (setq x (mincir ptlist))
  (princ "\n用时")
  (princ (* (- (getvar "TDUSRTIMER") t0) 86400)) ;结束计时
  (princ "秒")
  (if (null x)
    (alert "点的有效数目太小，请重新输入!")
    (progn
      (setq cen    (car x)
     radius (cadr x)
     ptmax  (caddr x)
      )
      ;;画圆及半径，列出圆的圆心半径值
      (entmake
 (append
   '((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle"))
   (list (cons 10 cen))
   (list (cons 40 radius))
   (list (cons 62 1))
 )
      )
      (entmake
 (append
   '((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine"))
   (list (cons 10 cen))
   (list (cons 11 ptmax))
   (list (cons 62 1))
 )
      )
      (list cen radius)
    )
  )
)
;;;************************************
;;;求最小包围圆的函数，空集返回空集，否
;;;则返回最小圆的圆心，半径和圆上的一点
;;;这是程序的主段----------------------
;;;************************************
(defun mincir (ptlist / CEN CEN_R P1 P2 P3 PTMAX R RADIUS X i)
  ;;判断有效点个数---------------------
  (cond
    ((= (length ptlist) 0)
     nil
    )
    ((= (length ptlist) 1)
     (alert "点集为一点,最小圆半径为0")
     (list (car ptlist) 0 (car ptlist))
    )
    ((= (length ptlist) 2)
     (alert "点集为两点，最小圆为过两点的圆")
     (setq cen   (mid (car ptlist) (cadr ptlist))
    radius (/ (distance (car ptlist) (cadr ptlist)) 2)
     )
     (list cen radius (car ptlist))
    )
    (t
     ;;开始递归运算----------------------------
     (setq p1 (car ptlist)
    p2 (cadr ptlist)
    p3 (caddr ptlist)
     )
     (setq cen_r (3pc p1 p2 p3))
     (setq ptmax (maxd-cir ptlist (car cen_r)))
     (setq i 0)
     (while (null (in1 ptmax (car cen_r) (cadr cen_r)))
       (setq cen_r (4pc p1 p2 p3 ptmax))
       (setq p1 (car (caddr cen_r))
      p2 (cadr (caddr cen_r))
      p3 (caddr (caddr cen_r))
       )
       (setq ptmax (maxd-cir ptlist (car cen_r)))
       (setq i (1+ i))
     )
     (list (car cen_r) (cadr cen_r) ptmax)
    )
  )
)
(defun make-line (p q)
  (entmake
    (list
      '(0 . "LINE")
      (cons 10 p)
      (cons 11 q)
    )
  )
)
;;以下代码来自晓东
;;定义取点函数----
(defun ssgetpoint (ss / i l a b c)
  (setq i 0)
  (if ss
    (repeat (sslength ss)
      (setq a (ssname ss i))
      (setq i (1+ i))
      (setq b (entget a))
      (setq c (cdr (assoc 10 b)))
      (setq l (cons c l))
    )
  )
  (reverse l)
)
(defun mid (p1 p2)
  (list
    (* (+ (car p1) (car p2)) 0.5)
    (* (+ (cadr p1) (cadr p2)) 0.5)
    (* (+ (caddr p1) (caddr p2)) 0.5)
  )
)
;;判断点是否在圆内------------------------
(defun in1 (pt cen r)
  (< (- (distance pt cen) r) 1e-8)
)
;;判断点集是否在圆内----------------------
(defun in2 (ptl cen r / pts pt)
  (setq pts ptl)
  (while (and (setq pt (car pts))
       (in1 pt cen r)
  )
    (setq pts (cdr pts))
  )
  (null pt)
)
;;定义三点最小圆圆心及其半径，若是锐角三角
;;形，则是其三点圆，否则是其最大边的直径圆
(defun 3pc (pa pb pc / D MIDPT)
  (cond  
    ((in1 pc (setq midpt (mid pa pb)) (setq d (/ (distance pa pb) 2)))
     (list midpt d (list pa pb pc))
    )
    ((in1 pa (setq midpt (mid pb pc)) (setq d (/ (distance pb pc) 2)))
     (list midpt d (list pb pc pa))
    )
    ((in1 pb (setq midpt (mid pc pa)) (setq d (/ (distance pc pa) 2)))
     (list midpt d (list pc pa pb))
    )
    (t
      (3pcircle pa pb pc)
    )
  )
)
;;; 三点圆函数
(defun 3PCirCle (P0 P1 P2 / X0 Y0 X1 Y1 X2 Y2 DX1 DY1 DX2 DY2 D 2D C1 C2 CE)
  (setq X0  (car  P0)
 Y0  (cadr P0)
 X1  (car  P1)
 Y1  (cadr P1)
 X2  (car  P2)
 Y2  (cadr P2)
 DX1 (- X1 X0)
 DY1 (- Y1 Y0)
 DX2 (- X2 X0)
 DY2 (- Y2 Y0)
  )
  (setq D (- (* DX1 DY2) (* DX2 DY1)))
  (if (/= D 0.0)
    (progn
      (setq 2D (+ D D)
     C1 (+ (* DX1 (+ X0 X1)) (* DY1 (+ Y0 Y1)))
     C2 (+ (* DX2 (+ X0 X2)) (* DY2 (+ Y0 Y2)))
     CE (List (/ (- (* C1 DY2) (* C2 DY1)) 2D)
       (/ (- (* C2 DX1) (* C1 DX2)) 2D)
        )
      )
      (list CE (distance CE P0) (list p0 p1 p2))
    )
  )
)
;;定义四点的最小圆圆心半径，并返回三点坐标
(defun 4pc (p1 p2 p3 ptmax / pts mind minr r 4ps)
  (setq pts (list (3pc p1 p2 ptmax)
    (3pc p1 p3 ptmax)
    (3pc p2 p3 ptmax)
     )
  )
  (setq 4ps (list p1 p2 p3 ptmax))
  (setq minr 1e308)
  (foreach n pts
    (setq r (cadr n))
    (if (and (< r minr)
      (in2 4ps (car n) r)
 )
      (setq mind n)
    )
  )
  mind
)
;;定义求点集中离圆心最远的点的函数--------
(defun maxd-cir (ptl cen / pmax dmax d)
  (setq dmax 0.0)
  (foreach pt ptl
    (if (> (setq d (distance pt cen)) dmax)
      (setq dmax d
     pmax pt
      )
    )
  )
  pmax
)

clivx · 发表于 2010-7-19 19:21

曲高，难和啊

jh1005 · 发表于 2009-2-26 22:12

谢谢版主，我用MEASURE命令把多义线等分成N个点，再配合你的程序，终于可以画出任意多义线的最小外接圆了，误差很少，以后就不用天天画这个圆了，嘿嘿！不过用MEASURE描成N个点后有点慢。

highflybir · 发表于 2009-2-27 02:41

原以为这样的文章可能没多少实际用途，被冷藏了这么久。想不到还有点用处啊。

谢谢版主，我用MEASURE命令把多义线等分成N个点，再配合你的程序，终于可以画出任意多义线的最小外接圆了，误差很少，以后就不用天天画这个圆了，嘿嘿！不过用MEASURE描成N个点后有点慢。

下面的程序可以给你解决这个烦恼。

只需要输入test命令，然后你选择一个或者多个多段线，就给你画出来了。不用你去measure了。

另外因为没有仔细考虑优化的问题，所以精度最好不要输入过大（不输入也可以），测试1000和30000的结果都相差不大。

qjchen · 发表于 2010-7-16 17:01

:) 谢谢highlfybird版主的精彩代码

顺便学习了版主的C版代码

http://www.objectarx.net/forum.php?mod=viewthread&tid=4870

谢谢

winabcd · 发表于 2010-7-22 15:42

好东西，回去好好研究一下下。

hgf876 · 发表于 2010-7-28 11:05

呵呵，感觉不错，

cxjzxh · 发表于 2010-9-2 11:22

思路不错，值得学习

cumtjh · 发表于 2010-9-4 00:04

思路不错，值得学习
<script type="text/javascript">var reload=1;</script>

账号		自动登录	找回密码
密码			注册

[越飞越高] 【飞鸟集】最小包围圆的最佳算法

本帖子中包含更多资源

评分

本帖被以下淘专辑推荐:

相关帖子

评分

本帖子中包含更多资源