确定包含平面点集的最小直径圆

snoopychen

最近有个朋友说有这样一道题
平面上有许多点，如何确定包含这些点集的最小圆
用了两个小时写下如下一些代码，见笑了
主要用了menzi的凸包函数还有一个自己想的还不成熟的想法。
后来查了一下书，原来是计算几何方面的内容，于是去借了一本周培德的《计算几何算法》第二版和一本M de berg的计算几何，发现原来计算几何这个坑这么深
eachy版主，xiaolongxin还有之前的杭州都在很久前就在晓东讨论了，要好好学习一下了～
写的仓促，验证不足，请批评
找外接圆的时候用command，有点懒
用的算法是简单的穷举，希望各位高手提出高效率的算法，谢谢。
[pcode=lisp,true]
;;; ========================================================================
;;; Thanks to Menzi's great code for find the convex hull of points        ;
;;; ========================================================================
;;; Some of the following codes are writen by QJCHEN                       ;
;;; Civil engineering Department, South China University of Technology     ;
;;;                                       ;
;;; Purpose: To find the minimum radius circle that contain a set of points;
;;;          in z=0 plane                                                  ;
;;; Algorithm: First, use Menzi's great code to find the convex hull of the;
;;;           points. Then, in the new point list, repeat each three points;
;;;           get the minimum radius circle that contain each three points ;
;;;           then the maximum radius cirlce of this set is the minimum    ;
;;;           radius circle that contain the points.                       ;
;;; The second step algorithm is my own thought, and havent been proved    ;
;;;           strictly, but it seems right, I will find a more strict way  ;
;;; The efficiency of the second step is slowly but after all it can run   ;
;;; Version: 0.1                                                           ;
;;; Platform: R2000 and after                                              ;
;;; Original post :www.Theswamp.org                                        ;
;;; 2006.06.14                                                             ;
;;; Thanks to Menzi's great code about the convex hull which I get from    ;
;;; Google group                                                           ;
;;; His website is http://www.menziengineering.ch/                         ;
;;; ========================================================================

(defun c:test (/ os ssp ssp1 ssp2 mmmax)
  (setq os (getvar "osmode"))
  (setvar "osmode" 0)
  (setvar "pdmode" 2)
  (setq ssp (ssget '((0 . "POINT"))))
  (setq ssp1 (ssgetpoint ssp))
  (setq ssp2 (MeGetConvexHull ssp1))
  (setq mmmax (getmin ssp2))
  (command "circle" (nth 1 mmmax) (nth 0 mmmax))
  (setvar "osmode" os)
)


;;;repeat the convex hull by each 3 points to find the circle
(defun getmin (lst / i n a j b k c mm lstmm mmmax)
  (setq i 0
mmmax (list nil nil)
  )
  (setq n (length lst))
  (repeat (- n 2)
    (setq a (nth i ssp2))
    (setq j (1+ i))
    (repeat (- (- n i) 2)
      (setq b (nth j ssp2))
      (setq k (1+ j))
      (repeat (- (- n j) 1)
(setq c (nth k ssp2))
(setq mm (3pcircle a b c))
(if (> (nth 0 mm) (nth 0 mmmax))
  (setq mmmax mm)
)        
(setq k (1+ k))
      )
      (setq j (1+ j))
    )
    (setq i (1+ i))
  )
  mmmax
)

;;; get the minimum circle of point a b c
(defun 3pcircle (a b c / ab bc ca index lst)
  (setq ab (distance a b))
  (setq bc (distance b c))
  (setq ca (distance a c))

  (if (and
(> (cosabc ab bc ca) 0)
(> (cosabc bc ab ca) 0)
(> (cosabc ca ab bc) 0)
      )
    (setq index 1)
    (setq index 0)
  )

  (if (= index 0)
    (setq lst (longestabc a b c ab bc ca))
  )
  (if (= index 1)
    (setq lst (longestabc1 a b c))
  )
  lst
)

;;; get the angle
(defun cosabc (a b c / d)
  (setq d (- (+ (* b b) (* c c)) (* a a)))
  d
)

;;; get the minimum circle of angle big than 90 degreed.
(defun longestabc (a b c ab bc ca / d r p)
  (setq d (max
    ab
    bc
    ca
  )
r (* 0.5 d)
  )
  (cond
    ((= d ab)
      (setq p (midpoint a b))
    )
    ((= d bc)
      (setq p (midpoint b c))
    )
    ((= d ca)
      (setq p (midpoint c a))
    )
  )
  (list r p)
)

;;; get the minimum circle of angle low than 90 degreed.
(defun longestabc1 (a b c)
  (command "circle" "3p" a b c)
  (setq d (entlast))
  (setq e (entget d))
  (setq p (cdr (assoc 10 e)))
  (setq r (cdr (assoc 40 e)))
  (entdel d)
  (list r p)
)
;;; midpoint function
(defun midpoint (p1 p2)
  (mapcar
    '(lambda (x)
       (/ x 2.)
     )
    (mapcar
      '+
      p1
      p2
    )
  )
)

(defun ssgetpoint (ss / i listpp a b c)
  (setq i 0
listpp nil
  )
  (if ss
    (progn
      (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 i (1+ i))
      )
    )
  )
  listpp
)
;;;
;;; == Function MeGetConvexHull
;;; Calculates the convex hull from a 'cloud' of points.
;;; Arguments [Type]:
;;;   Lst = Point list [LIST]
;;; Return [Type]:
;;;   > Hull points [LIST]
;;; Notes:
;;;   - Algorithm by Graham
;;;
(defun MeGetConvexHull (Lst / FstCnt LstLen MinCnt NxtCnt TmpLst)
  (setq LstLen (length Lst)
MinCnt 0
FstCnt 1
  )
  (while (< FstCnt LstLen)
    (if (< (cadr (nth FstCnt Lst)) (cadr (nth MinCnt Lst)))
      (setq MinCnt FstCnt)
    )
    (setq FstCnt (1+ FstCnt))
  )
  (setq FstCnt 0)
  (while (< FstCnt LstLen)
    (if (and
  (equal (cadr (nth FstCnt Lst)) (cadr (nth MinCnt Lst)) 1E-8)
  (> (car (nth FstCnt Lst)) (car (nth MinCnt Lst)))
)
      (setq MinCnt FstCnt)
    )
    (setq FstCnt (1+ FstCnt))
  )
  (setq TmpLst (MeSwapList Lst 0 MinCnt)
TmpLst (vl-sort TmpLst (function (lambda (e1 e2)
   (< (MeCalcTheta (nth 0 TmpLst) e1)
      (MeCalcTheta (nth 0 TmpLst) e2)
   )
)
       )
       )
TmpLst (cons (last TmpLst) TmpLst)
NxtCnt 3
FstCnt 4
  )
  (while (<= FstCnt LstLen)
    (while (>= (MeGetCcw (nth NxtCnt TmpLst) (nth (1- NxtCnt) TmpLst)
(nth FstCnt TmpLst)
       ) 0
   )
      (setq NxtCnt (1- NxtCnt))
    )
    (setq NxtCnt (1+ NxtCnt)
  TmpLst (MeSwapList TmpLst FstCnt NxtCnt)
  FstCnt (1+ FstCnt)
    )
  )
  (mapcar
    '(lambda (l)
       (nth l TmpLst)
     )
    (MeRepeatedList 0 (1- NxtCnt) 1)
  )
)
;;;
;;; == Function MeSwapList
;;; Swaps 2 atoms in a list.
;;; Arguments [Type]:
;;;   Lst = List to swap [LIST]
;;;   Fst = Source position [INT]
;;;   Nxt = Target position [INT]
;;; Return [Type]:
;;;   > Swapped list [LIST]
;;; Notes:
;;;   None
;;;
(defun MeSwapList (Lst Fst Nxt / FstVal NxtVal TmpLst)
  (setq FstVal (nth Fst Lst)
NxtVal (nth Nxt Lst)
TmpLst (MeEditList 0 Fst NxtVal Lst)
  )
  (MeEditList 0 Nxt FstVal TmpLst)
)
;;;
;;; == Function MeCalcTheta
;;; Calculates a ordinal number between 2 points.
;;; Arguments [Type]:
;;;   Pt1 = First point [LIST]
;;;   Pt2 = Second point [LIST]
;;; Return [Type]:
;;;   > A value between 0 and 360 [REAL]
;;; Notes:
;;;   None
;;;
(defun MeCalcTheta (Pt1 Pt2 / X__Abs Y__Abs X__Dif Y__Dif TheVal)
  (setq X__Dif (- (car Pt2) (car Pt1))
Y__Dif (- (cadr Pt2) (cadr Pt1))
X__Abs (abs X__Dif)
Y__Abs (abs Y__Dif)
TheVal (if (equal (+ X__Abs Y__Abs) 0 1E-5)
0
(/ Y__Dif (+ X__Abs Y__Abs))
       )
  )
  (if (< X__Dif 0)
    (setq TheVal (- 2.0 TheVal))
    (if (< Y__Dif 0)
      (setq TheVal (+ 4.0 TheVal))
    )
  )
  (* 90.0 TheVal)
)
;;;
;;; == Function MeGetCcw
;;; Determines the direction of 3 points.
;;; Arguments [Type]:
;;;   Pt1 = First point [LIST]
;;;   Pt1 = Second point [LIST]
;;;   Pt3 = Third point [LIST]
;;; Return [Type]:
;;;   >  1 = ccw [INT]
;;;   > -1 = cw [INT]
;;;   >  0 = Colinear [INT]
;;; Notes:
;;;   None
;;;
(defun MeGetCcw (Pt0 Pt1 Pt2 / X1_Dif X1_Sqr X2_Dif X2_Sqr Y1_Dif Y1_Sqr
     Y2_Dif Y2_Sqr
)
  (setq X1_Dif (- (car Pt1) (car Pt0))
Y1_Dif (- (cadr Pt1) (cadr Pt0))
X2_Dif (- (car Pt2) (car Pt0))
Y2_Dif (- (cadr Pt2) (cadr Pt0))
X1_Sqr (* X1_Dif X1_Dif)
Y1_Sqr (* Y1_Dif Y1_Dif)
X2_Sqr (* X2_Dif X2_Dif)
Y2_Sqr (* Y2_Dif Y2_Dif)
  )
  (cond
    ((> (* X1_Dif Y2_Dif) (* Y1_Dif X2_Dif))
      1
    )
    ((< (* X1_Dif Y2_Dif) (* Y1_Dif X2_Dif))
      -1
    )
    ((or
       (< (* X1_Dif X2_Dif) 0)
       (< (* Y1_Dif Y2_Dif) 0)
     )
      -1
    )
    ((< (+ X1_Sqr Y1_Sqr) (+ X2_Sqr Y2_Sqr))
      1
    )
    (0)
  )
)
;;;
;;; == Function MeRepeatedList
;;; Fills a list with numbers from Srt to End, using increment of Inc.
;;; Arguments [Type]:
;;;   Srt = Start number [INT] or [REAL]
;;;   End = End number [INT] or [REAL]
;;;   Inc = Increment to use [INT] or [REAL]
;;; Return [Type]:
;;;   > List containting all numbers between and including Srt and End
;;; Notes:
;;;   If End is not a repeated of End - Srt, End will not be included
;;;
(defun MeRepeatedList (Srt End Inc / TmpVal RetVal)
  (setq TmpVal (- Srt Inc))
  (while (< TmpVal End)
    (setq RetVal (append
   RetVal
   (list (setq TmpVal (+ TmpVal Inc)))
)
    )
  )
  RetVal
)
;;;
;;; == Function MeEditList
;;; Delete, replace or append a list.
;;; Arguments [Type]:
;;;   Mde = Mode  1 append  [INT]
;;;         Mode  0 replace [INT]
;;;         Mode -1 delete  [INT]
;;;   Pos = Position in list [INT]
;;;   Val = New value [INT/REAL/STR/LIST]
;;;   Lst = List [LIST]
;;; Return [Type]:
;;;   > Modified list [LIST]
;;; Notes:
;;;   - Delete and replace, require the position (Pos) in list.
;;;   - Append and replace, require the new value (Val).
;;;   - :vlax-null items are not allowed in the list argument
;;;
(defun MeEditList (Mde Pos Val Lst / Countr TmpLst TmpVal)
  (if (= Mde 1)
    (reverse (cons Val (reverse Lst)))
    (progn
      (setq Countr -1
    TmpVal (if (= Mde -1)
     :vlax-null
     Val
   )
    TmpLst (mapcar
     '(lambda (l)
(if (= Pos (setq Countr (1+ Countr)))
  TmpVal
  l
)
      )
     Lst
   )
      )
      (vl-remove :vlax-null TmpLst)
    )
  )
)

[/pcode]

xiao_longxin

1.先求点集的凸壳
2.在凸壳中求最远距离的两个点
3.以此两点中心为圆心，距离一半为半径作圆。
只是想像，未验证

雨箭风刀

楼上的方法，理论上是错误的

snoopychen

假如凸包是一个等边三角形的话，xiaolongxin兄的方法似乎是不对的

雨箭风刀

P0-P1 距离大于 p0－p3,
p3未能包含在p0 p1做成的园内