[原创]：最小包围圆的最佳算法

highflybird · 发表于 2006-11-11 23:28:27

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？立即注册

×

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

;;;************************************
;;;求最小包围圆的lisp程序--------------
;;;其算法为参见了有关文献--------------
;;;这种算法在退化很严重的情况结果也正确
;;;其中程序主段是核心算法，其他的附加程
;;;序为取点，画点，画圆和其半径，可舍弃
;;;************************************
(defun C:mincir (/ p1 p2 p3 ptmax cen_r cen radius oldmode oce x)
;;定义错误函数和预处理---------------
(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)
;;定义取点函数-----------------------
(defun getp (/ m)
(setq m (fix (getint "\n请输入点数目:")))
(if (< m 1) nil
(progn
(setq m (1- m))
;;定义0-N的有序整数集函数------
(defun listn (a / n x)
(setq n 1)
(setq x (list 0))
(while (<= n a)
(setq x (append (list n) x))
(setq n (1+ n))
)
(reverse x)
)
(setq x (listn m))
;;取得点集---------------------
(setq x (mapcar '(lambda (x) (getpoint "\n输入点: ")) x))
)
)
)
;;定义中点函数，本来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))
)
)
;;取点，画点-------------------------
(command ".ucs" "W")
(setq x (vl-remove nil (getp)))
(mapcar
'(lambda (pt)
(entmake
(append
'((0 . "Point") (100 . "AcDbEntity") (100 . "AcDbPoint"))
(list (cons 10 pt))
)
)
)
x
)
;;判断有效点个数---------------------
(cond
((< (length x) 1)
(progn
(alert "点的有效数目太小，请重新输入!")
(command ".ucs" "p")
(setvar "osmode" oldmode)
(setvar "cmdecho" oce)
nil
)
)
((= (length x) 1)
(progn
(alert "点集合为一点,最小圆半径为0")
(command ".ucs" "p")
(setvar "osmode" oldmode)
(setvar "cmdecho" oce)
(list (car x) 0)
)
)
((= (length x) 2)
(progn
(alert "点集合为两点，最小圆直径为其两点距离，\n圆心为其连线中点")
(setq cen (midp (car x) (cadr x)) radius (/ (distance (car x) (cadr x)) 2))
(entmake
(append
'((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine"))
(list (cons 10 cen))(list (cons 11 (car x)))(list (cons 62 1))
)
)
(entmake
(append
'((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle"))
(list (cons 10 cen))(list (cons 40 radius))(list (cons 62 1))
)
)
(command ".ucs" "p")
(setvar "osmode" oldmode)
(setvar "cmdecho" oce)
(list cen radius)
)
)
(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)))
(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 x) p2 (cadr x) p3 (caddr x))
(setq cen_r (3pc p1 p2 p3))
(setq ptmax (maxd-cir x (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 x (car cen_r)))
)
;;画圆，画半径，并列出最小圆的圆心、半径值
(entmake
(append
'((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle"))
(list (cons 10 (car cen_r)))(list (cons 40 (cadr cen_r)))(list (cons 62 1))
)
)
(entmake
(append
'((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine"))
(list (cons 10 (car cen_r)))(list (cons 11 ptmax))(list (cons 62 1))
)
)
(command ".ucs" "p")
(setvar "osmode" oldmode)
(setvar "cmdecho" oce)
(list (car cen_r) (cadr cen_r))
);;for progn
);; for t
);; for cond
);; for defun

;;;************************************ 
;;;求最小包围圆的lisp程序-------------- 
;;;其算法为参见了有关文献-------------- 
;;;这种算法在退化很严重的情况结果也正确 
;;;其中程序主段是核心算法，其他的附加程 
;;;序为取点，画点，画圆和其半径，可舍弃 
;;;************************************ 
(defun C:mincir        (/ p1 p2 p3 ptmax cen_r cen radius oldmode oce x) 
  ;;定义错误函数和预处理--------------- 
  (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) 
  ;;定义取点函数----------------------- 
  (defun getp (/ m) 
    (setq m (fix (getint "\n请输入点数目:"))) 
    (if        (< m 1) nil 
      (progn 
        (setq m (1- m)) 
        ;;定义0-N的有序整数集函数------ 
        (defun listn (a / n x) 
          (setq n 1) 
          (setq x (list 0)) 
          (while (<= n a) 
            (setq x (append (list n) x)) 
            (setq n (1+ n)) 
          ) 
          (reverse x) 
        ) 
        (setq x (listn m)) 
        ;;取得点集--------------------- 
        (setq x (mapcar '(lambda (x) (getpoint "\n输入点: ")) x)) 
      ) 
    ) 
  ) 
  ;;定义中点函数，本来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)) 
    ) 
  ) 
  ;;取点，画点------------------------- 
  (command ".ucs" "W") 
  (setq x (vl-remove nil (getp))) 
  (mapcar 
    '(lambda (pt) 
       (entmake 
         (append 
           '((0 . "Point") (100 . "AcDbEntity") (100 . "AcDbPoint")) 
           (list (cons 10 pt)) 
         ) 
       ) 
     ) 
    x 
  ) 
  ;;判断有效点个数--------------------- 
  (cond 
    ((< (length x) 1) 
     (progn 
       (alert "点的有效数目太小，请重新输入!") 
       (command ".ucs" "p") 
       (setvar "osmode" oldmode) 
       (setvar "cmdecho" oce) 
       nil 
     ) 
    ) 
    ((= (length x) 1) 
     (progn 
       (alert "点集合为一点,最小圆半径为0") 
       (command ".ucs" "p") 
       (setvar "osmode" oldmode) 
       (setvar "cmdecho" oce) 
       (list (car x) 0) 
     ) 
    ) 
    ((= (length x) 2) 
     (progn 
       (alert "点集合为两点，最小圆直径为其两点距离，\n圆心为其连线中点") 
       (setq cen (midp (car x) (cadr x)) radius (/ (distance (car x) (cadr x)) 2)) 
       (entmake 
         (append 
           '((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine")) 
           (list (cons 10 cen))(list (cons 11 (car x)))(list (cons 62 1)) 
         ) 
       ) 
       (entmake 
         (append 
           '((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle")) 
           (list (cons 10 cen))(list (cons 40 radius))(list (cons 62 1)) 
         ) 
       ) 
       (command ".ucs" "p") 
       (setvar "osmode" oldmode) 
       (setvar "cmdecho" oce) 
       (list cen radius) 
     ) 
    ) 
    (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))) 
         (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 x) p2 (cadr x) p3 (caddr x)) 
       (setq cen_r (3pc p1 p2 p3)) 
       (setq ptmax (maxd-cir x (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 x (car cen_r))) 
       ) 
       ;;画圆，画半径，并列出最小圆的圆心、半径值 
       (entmake 
         (append 
           '((0 . "circle") (100 . "AcDbEntity") (100 . "AcDbCircle")) 
           (list (cons 10 (car cen_r)))(list (cons 40 (cadr cen_r)))(list (cons 62 1)) 
         ) 
       ) 
       (entmake 
         (append 
           '((0 . "line") (100 . "AcDbEntity") (100 . "AcDbLine")) 
           (list (cons 10 (car cen_r)))(list (cons 11 ptmax))(list (cons 62 1)) 
         ) 
       ) 
       (command ".ucs" "p") 
       (setvar "osmode" oldmode) 
       (setvar "cmdecho" oce) 
       (list (car cen_r) (cadr cen_r)) 
      );;for progn 
    );;  for t    
  );;    for cond 
);;      for defun

附件为其lisp程序，加载，运行mincir 即可。

snoopychen · 发表于 2006-11-12 09:48:41

:)

前段时间也匆匆忙忙写过一个
http://www.xdcad.net/forum/showthread.php?s=&threadid=571449

但是是赶时间,用比较暴力的方法,比不上兄台的优质算法.
《计算几何-算法与应用》确实是本好书, 计算几何确实挺有趣.

希望能再多共享一些这样的程序,谢谢:)

fsxm · 发表于 2006-11-12 11:38:29

highflybird ^)^写的不错！你与snoopychen 老大一起给我上了很多非常好的几何课！
在这之前我要将几何忘记的差不多了~多谢两位！

提点个人的小意见一家之言只供参考！
在你的子函数4pc中算法可以优化~
3pc调用次数太多了可能是这个程序主要费时函数！
而事实上在求一点pt是不是在过a b c 的圆的内，外，上可以不用求出cen,r的啦！
比如说用同弧所对圆周角等于弧角度一半，大于则在圆内，小于则在外~
p在圆上角apc=角abc
p在圆外角apc<角abc
p在圆内角apc>角abc
角的取值在什么范围也是要处理一下！
我没测验理论上可以速率比求cen&R的要快点
事实上快不快就只有测试下知道啦！只供参考！

hhc_1 · 发表于 2006-11-12 15:56:53

如果结果只有两个点在圆上,这两个点一定是距离最远的点,位于直径上

highflybird · 发表于 2006-11-12 17:28:35

多谢各位，我学lisp不久，向大家学习的地方很多。
fsxm 兄，很高兴能与我讨论。现我回答如下：
先说说程序思路：
1、找出三点ABC最小圆和圆心O。
2、再找出点集中离这个圆心O距离最远的那一点PTMAX。
3、如果这点在三点的最小圆内，这个圆为所求，算法结束。
4、否则三点和这个最远点PTMAX组成四点，求四点最小圆和圆心，并返回三点ABC（圆周上两点和另一点）
5、再求出离四点最小圆的圆心最远的点PTMAX，如果这点在四点最小圆内，则算法结束，否则以返回的三点和最远点PTMAX迭代进入第4 步。直到最远点PTMAX在跌代后的最小圆内。
因此程序中，3pc这个子函数是必不可少的，有了3pc，才有4pc，递推才能进行。
实际上迭代的次数并不是所取的点数n（可能几万个以上），跌代的次数很少，这种递归形成的级数收敛很快，开始前三点是随机选择的，如果选择得好的话，只要一两次迭代就足矣。
另外，我判断点在圆内圆外，并不是这个函数，而是根据：如果圆心到点距离大于半径则不在圆内，否则在圆内。是in1 和in2 这两个函数，in2 用来判断点集是否在圆内，但并不是判断点集合所有的元素，而是判断递归中形成的四点。
所以真正调用次数太多的并不是3pc，而是到圆心距离最远点的那个函数maxd-cir，因为这个函数每递归一次，就要判断点集合中的点到圆心得最远距离（实际上可以除去一些点，但那样的话又得把点集中的点判断一次），影响时间的是这个。我个人以为。
至于4pc这个函数，我也不知道如何优化，我是这样设计的：
对四点三点一组形成集合，每一组求其最小圆和半径(3pc) ，然后对半径按照从小到大排序，找到能包围四点的那一组。
希望大家能帮忙看看何处可以优化。

;;;************************************
;;;求最小包围圆的lisp程序--------------
;;;其算法为参见了有关文献--------------
;;;这种算法在退化很严重的情况结果也正确
;;;其中程序主段是核心算法，其他的附加程
;;;序为取点，画点，画圆和半径，用来测试
;;;************************************
(defun C:test (/ olderr en errmsg oldmode oce ssp 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 (append listpp (list c)))
(setq listpp (cons c listpp))
(setq i (1+ i))
)
)
listpp
)

;;;************************************ 
;;;求最小包围圆的lisp程序-------------- 
;;;其算法为参见了有关文献-------------- 
;;;这种算法在退化很严重的情况结果也正确 
;;;其中程序主段是核心算法，其他的附加程 
;;;序为取点，画点，画圆和半径，用来测试 
;;;************************************ 
(defun C:test (/ olderr en errmsg oldmode oce ssp 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 (append listpp (list c))) 
      (setq listpp (cons c listpp)) 
      (setq i (1+ i))   
    ) 
  ) 
  listpp 
)

我重新修改了程序，使得求最小圆成了一个函数。
然后在做了一个test 程序，我选择了50000个点，然后对其求值
发现，竟然取得点集合的那个时间占了绝大部分时间，而求最小圆的函数计算时间很短。
如果大家有其他算法的不妨作下比较（把局部变量变成全局的，然后从取得点集合后开始做比较）

eachy · 发表于 2006-11-12 18:44:05

这样改取点函数后试试
(setq listpp (append listpp (list c))) -> (setq listpp (cons c listpp)

highflybird · 发表于 2006-11-12 19:14:36

版主改得太妙了！（尽管我还不知道为什么）
我把snoopychen老大提供的暴力算法和我的算法作了下比较：（没别的意思）简直有天壤之别。
理论上用暴力算法的时间是O(n^4)，可见是很费时的。
看来除非万不得已最好莫用暴力算法。
向各位学习！

aeo · 发表于 2006-11-12 22:40:32

是不是一个最小圆，里面圈住所有的点?

highflybird · 发表于 2006-11-13 17:39:21

是这样的。

对于这个程序，只要再作稍微改动
1、把取得点集函数改成如下：

(defun ssgetpoint1 (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 d (cdr (assoc 11 b)));;直线端点
(setq listpp (cons c listpp))
(setq listpp (cons d listpp))
(setq i (1+ i))
)
)
listpp
)

在TEST中修改成选取直线对象：

(setq ssp (ssget '((0 . "LINE"))))
(setq ptlist (ssgetpoint1 ssp))

那么这个程序就可以成为直线集合的最小包围圆了。
我记得在明经上有人问过：如何求包含一个多边形的最小的圆?
那么这个程序基本可以解答了。

对于曲线（包括圆弧）集合，也可以通过对每个曲线用有限的点（不要太多）进行拟合，然后用这种方法来得到粗略解，然后细化，估计也行。但我想这样做肯定不是最佳解了。

hhc_1 · 发表于 2006-11-14 14:55:09

曾经发现:同样的功能如果不用mapcar这个函数,速度要快不少,不知是否都是这样

Highflybird · 发表于 2013-4-9 17:02:09

本帖最后由 Highflybird 于 2013-4-9 17:03 编辑

程序稍微做了一下改进。日期：2013-4-9.参见附件：

请点击此处下载

查看状态：需购买或无权限

您的用户组是：游客

文件名称:最小包围圆改进.rar
下载次数:112 文件大小:2.16 KB
下载权限: 不限以上 [免费赚D豆]

其实我以前还有个随机增量算法版本的。但是相比这个，在大数量时候还是慢了很多，故没贴出来。

adolfken · 发表于 2013-4-17 08:19:09

高手果然是高手啊

cuyongping · 发表于 2013-4-27 09:23:05

看不懂！还有待提高！

sicky111 · 发表于 2013-5-3 02:17:15

学习的好例子。

守仁格竹GM · 发表于 2013-5-22 11:25:22

还有一些问题

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

[原创]：最小包围圆的最佳算法

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

相关帖子