;;;************************************************************************
;;;一個求點集合的凸包的lisp程序--------------------------------------------
;;;------采用的算法為禮品包扎法--------------------------------------------
;;;方法為最右端的點開始處理，將該點作為凸包邊界的第一個點P1，從最初的垂直線
;;;方向繞P1順時針旋轉，直到碰到另一個P2這就是凸包邊界的第二個點P2，依此類推
;;;p2求得p3......直到又重新回到p1，已經考慮了各種退化情況和浮點運算，其算法
;;;時間不超過O(n.h),其中h是凸包的復雜度，時間還是很快的。大家不妨驗證。    
;;;參考文獻<<計算幾何-算法及其應用>>(第二版),以及參考了其他網站的一些源代碼
;;;------------------------------------------------------------------------
;;;其中程序主段是核心算法，其他的附加程序為取得點集，畫凸包邊界線，測試大量
;;;點集函數處理所花費的時間。----------------------------------------------
;;;用法: 加載lisp運行test選取點，直線段，或多義線(全是直線段組成)即可。----
;;;************************************************************************
(defun C:test (/ sel t0 ptlist pp 2Pi)
  (setq 2Pi (* 2 pi))
  (setq	sel (ssget (list '(0 . "POINT"))))                  ;選擇點集
  (if (setq ptlist (getpt sel))                             ;構造點集
    (progn
      (setq t0 (getvar "TDUSRTIMER"))                       ;開始計時
      (setq pp (hull ptlist))                               ;求凸包
      (princ "\n用時")  
      (princ (* (- (getvar "TDUSRTIMER") t0) 86400))        ;結束計時
      (princ "秒")
    )
  )
  (if (null pp)
    (alert "點的有效數目太小，請重新輸入!")
    (entmake                                                ;畫凸包
      (append
       '( (0 . "LWPOLYLINE")
	  (100 . "AcDbEntity")
	  (100 . "AcDbPolyline")
	)
	(list (cons 90 (length pp)))                        ;頂點個數
	(mapcar '(lambda (x) (cons 10 x)) pp)               ;多段線頂點
	(list (cons 70 1))                                  ;閉合的
	(list (cons 62 1))                                  ;紅色的
      )
    )
  )
  (princ)
)
;;;==========================
;;;程序主段，可以單獨成為函數
;;;==========================

;;;右半部的凸包
(defun hull1 (pts MaxPt MinPt / nextPt hullPt)
  (if pts 
    (progn
      (setq nextPt (Max-angle1 pts MaxPt))                  ;從最上面的點開始
      (setq hullPt (cons nextPt (cons MaxPt hullPt)))       ;順時針求得第一點
      (while (not (equal nextPt MinPt 1e-8))                ;到最下面的點為止
        (setq nextPt (Max-angle pts nextPt))                ;循環求凸包每一點
	(setq hullPt (cons nextPt hullPt))                  ;把每點加入凸包集
      )
    )
  )
)
;;;左半部的凸包
(defun hull (ptlist / revPts 2ndPts maxYp1 minYp1 maxYp2 minYp2
	              ptlst1 ptlst2 +ptlst -ptlst hullp1 hullp2)
  (setq revPts (mapcar 'reverse ptlist))                    ;點表的X和Y交換
  (setq 2ndPts (mapcar 'cadr ptlist))                       ;點表的Y值的表
  (setq maxYp1 (reverse (assoc (apply 'max 2ndPts) revPts)));最上面的點
  (setq minYp1 (reverse (assoc (apply 'min 2ndPts) revPts)));最下面的點
  (setq maxYp2 (list (- (car maxYp1)) (cadr maxYp1)))       ;鏡像后最上面的點
  (setq minYp2 (list (- (car minYp1)) (cadr minYp1)))       ;鏡像后最下面的點
  (foreach n ptlist                                         ;把點表分成兩部分
    (if (> (det minYp1 n maxYp1) 0)                         ;如果左轉
      (setq ptlst1 (cons n ptlst1))                         ;加入右半部分
      (setq ptlst2 (cons n ptlst2))                         ;否則左半部分
    )
  )
  (setq +ptlst (cons minYp1 (cons maxYp1 ptlst1)))          ;右半部分
  (setq -ptlst (Mirror-list-X ptlst2))                      ;左半部分以Y軸鏡像
  (setq hullp1 (hull1 +ptlst maxYp1 minYp1))                ;右半部分的凸包
  (setq hullp2 (hull1 -ptlst maxYp2 minYp2))                ;左半部分鏡像的凸包
  (setq hullp2 (cdr (reverse (cdr hullp2))))
  (setq hullp2 (Mirror-list-X hullp2))                      ;左半部分的凸包
  (append hullp1 hullp2)                                    ;把凸包左右相加
)
;;;鏡像左半部分
(defun Mirror-list-X (ptlist)
  (mapcar (function (lambda (x)(list (- (car x))(cadr x)))) ptlist)
)
;;;求點集中夾角的最大值的點
(defun Max-angle (ptlist pt / An)
  (setq An (mapcar (function (lambda (x) (angle pt x))) ptlist))
  (nth (- (length An) (length (member (apply 'max An) An))) ptlist)
)
(defun Max-angle1 (ptlist pt / An)
  (setq An
    (mapcar
      (function
        (lambda (x)
	  (if
	    (and
	      (equal (cadr x) (cadr pt) 1e-8)
	      (> (car x) (car pt))
	    )
	    (+ 2Pi (- (car x) (car pt)))
	    (angle pt x)
	  )	   
	)
      )
      ptlist
    )
  )
  (nth (- (length An) (length (member (apply 'max An) An))) ptlist)
)
;;定義三點的行列式,即三點之倍面積
(defun det (p1 p2 p3 / dx1 dy1 dx2 dy2)
  (setq	dx1 (- (car  p2) (car  p1))
	dy1 (- (cadr p2) (cadr p1))
	dx2 (- (car  p3) (car  p1))
	dy2 (- (cadr p3) (cadr p1))
  )
  (- (* dx1 dy2) (* dx2 dy1))
)				     
;;;============
;;;程序主段結束
;;;============

;;;依據曉東網站的代碼改寫而成的取點函數
(defun getpt (ss / a b c d i p)
  (setq	i 0)
  (if ss
    (repeat (sslength ss)
      (setq a (ssname ss i))
      (setq b (entget a))
      (setq c (cdr (assoc 10 b)))
      (setq c (list (car c) (cadr c)))
      (setq p (cons c p))
      (setq i (1+ i))
    )
  )
  p
)