点集的直径.lsp

计算几何中的几种lisp语言的算法

;;;************************************************************************
;;;一个求点集合的凸包的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"))))                  ;选择点集
  (setq sel (ssget))
  (if (setq ptlist (getpt sel))                             ;构造点集
    (progn
      (setq t1 (getvar "CDATE"))                            ;计时开始
      ;;(setq t0 (getvar "TDUSRTIMER"))                       ;开始计时
      (setq p1 (hull ptlist))                               ;求凸包
      (setq t2 (getvar "CDATE"))                            ;计时结束
      (setq pp (Max-distance p1))
      (setq t3 (getvar "CDATE")) 
      (princ "\n求点集的凸包用时:")
      (princ (* (- t2 t1) 1e6))
      (princ "秒")
      (princ "\n凸包共有")
      (princ (length p1))
      (princ "个顶点")
      (princ "\n求凸包的直径用时:")
      (princ (* (- t3 t2) 1e6))
      (princ "秒")
      (princ "\n总用时=最大距离点对用时:")
      (princ (* (- t3 t1) 1e6))
      (princ "秒")      
      ;;(princ (* (- (getvar "TDUSRTIMER") t0) 86400))        ;结束计时
    )
  )
  (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 0))                                  ;闭合的
	(list (cons 62 1))                                  ;红色的
      )
    )
  )
  (princ)
)
;;;==========================
;;;程序主段，可以单独成为函数
;;;==========================

;;;右半部的凸包
(defun hull1 (pts MaxPt MinPt / nextPt hullPt)
  (if (< (length pts) 3)
    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 Max-distance (ptlist / maxD halfPi HullP1 l HullP2 midPts
		              i j Pi+1 Qi+1 Ai+1 D-i PairPt)
  (setq MaxD nil)
  (setq 2Pi  (* 2 pi))
  (setq halfPi (/ Pi 2))
  (setq HullP1 (Hull ptlist))
  (setq l      (1- (length HullP1)))
  (setq HullP2 (append (cdr HullP1) (list (car HullP1))))
  (setq midPts (mapcar 'mid-pt HullP1 HullP2))
  (setq i 1)
  (foreach pt (reverse (cdr (reverse HullP1)))
    (setq j i)
    (setq Pi+1 (nth i HullP1))
    (setq Qi+1 (nth i midPts))
    (setq Ai+1 (ang Qi+1 pt Pi+1))
    (while (and (< Ai+1 halfPi) (< j l))
      (setq j (1+ j))
      (setq Ai+1 (ang (nth j midPts) Pt (nth j HullP1)))  
    )
    (setq D-i  (distance pt (nth j HullP1)))
    (setq MAXD (cons (list D-i (1- i) j) MAXD)) 	  
    (setq i (1+ i))  
  )
  (setq PairPt (assoc (apply 'Max (mapcar 'car MaxD)) MaxD))
  (list
    (nth (cadr  PairPt) HullP1)
    (nth (caddr PairPt) HullP1)
  )
)
(defun mid-pt (p1 p2)
  (list
    (* (+ (car  p1) (car  p2)) 0.5)
    (* (+ (cadr p1) (cadr p2)) 0.5)
  )
)
;;;============
;;;程序主段结束
;;;============

;;;依据晓东网站的代码改写而成的取点函数
(defun getpt1 (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
)
;;定义三点的夹角函数
(defun ang (p1 p2 p3 / x)
  (setq x (abs (- (angle p1 p3) (angle p1 p2))))
  (if (< x Pi) x (- 2Pi x)) 
)
(defun C:tt (/ p1 p2 p3)
  (initget 1)
  (setq p1 (getpoint "\n输入第一点:"))
  (setq p2 (getpoint "\n输入第二点:"))
  (setq p3 (getpoint "\n输入第三点:"))
  (ang p1 p2 p3)
)  
;;;取点函数2
(defun getpt (ss / i listpp a b c d)
  (setq	i 0)
  (if ss
    (repeat (sslength ss)
      (setq a (ssname ss i))
      (setq b (entget a))
      (setq ename (cdr (assoc 0 b)))
      (cond
	( (= ename "LWPOLYLINE")
	  (setq c (get-LWpolyline-vertexs b))
	  (setq listpp (append c listpp))
	)
	( (= ename "LINE")
	  (setq c (cdr (assoc 10 b)))
	  (setq d (cdr (assoc 11 b)))
	  (setq c (list (car c) (cadr c)))
	  (setq d (list (car d) (cadr d)))
	  (setq listpp (cons c listpp))
	  (setq listpp (cons d listpp))
	)
	( (= ename "POINT")
	  (setq c (cdr (assoc 10 b)))
	  (setq c (list (car c) (cadr c)))
	  (setq listpp (cons c listpp))
	)
      )
      (setq i (1+ i))
    )
  )
  listpp
)
;;取得多边形顶点
(defun get-LWpolyline-vertexs (entlst / n lst)
  (foreach n entlst
    (if (= (car n) 10)
      (setq lst (cons (cdr n) lst))
    )
  )
  lst
)