点集的直径和最小包围矩形

点集的最小包围矩形和点集的直径。最小包围盒指的是能包围一个多边形或者一些多边形的最小面积(可能周长)的这样的一个矩形。

下面是其LISP代码的实现:

  1. ;;;The procedure for Test
  2. (defun C:test (/ PP PTLIST SEL T0 n)
  3.   (princ "n命令是test(The command is Test).")
  4.   (princ "n请选择点、曲线类物体或者图块")
  5.   (setq sel (ssget '((0 . "POINT,LWPOLYLINE,LINE,SPLINE,ARC,CIRCLE,ELLIPSE,INSERT"))))  ;select curve or point
  6.   (initget 7)
  7.   (setq n (getint "n取样精度(100-2000):"))
  8.   (if sel                                                         
  9.     (progn
  10.       (setq ptlist (getpt sel n))                                       ;construct the set of points
  11.       (setq ptlist (Graham-scan ptlist))                                ;construct the CCW Hull of this set.
  12.       (if (<= (det (car ptlist) (cadr ptlist) (caddr ptlist)) 0.0)      ;ensure the hull is CCW.
  13.         (setq ptlist (reverse ptlist))                                  ;if it isn't CCW,then reverse it
  14.       ) 
  15.       (setq t0 (getvar "TDUSRTIMER"))                                   ;The start time of this algorithm
  16.       (setq pp (car (MinAreaRectangle ptlist)))                         ;start calculating
  17.       (princ "nIt takes :")                                            
  18.       (princ (* (- (getvar "TDUSRTIMER") t0) 86400))                    ;The End time
  19.       (princ "seconds")
  20.       (if pp
  21.         (make-poly pp)                                                  ;draw rectangle.
  22.       )
  23.     )
  24.   )
  25.   (princ)
  26. )
  27. ;;;=======================================================
  28. ;;;Function : Find the minimum area of encasing rectangle.
  29. ;;;Arguments : A CCW HULL                                 
  30. ;;;Return: The Four points of Rectangle and its Area      
  31. ;;;=======================================================
  32. (defun MinAreaRectangle (ptlist      /     AA    AI    BB    D1
  33.                          D2    EDGE  I     I1X   I1Y   I2X   I2Y
  34.                          IL    INF   IX    IY    J1    J2    MINA
  35.                          MINH  MINW  NORH  NORM  PI1   PI2   PTI0
  36.                          PTI1  PTI2  PTJ1  PTK1  PTM1  PTS1  PTS2
  37.                          PTS3  PTS4  REC1  REC2  REC3  REC4  RECT
  38.                          VECH  VECL  VJ12  VM12
  39.                         )
  40.   (setq INF 1e309)                                                      
  41.   (setq minA INF)                                                       ;Initiating the Minimum area is infinite
  42.   (setq pti0 (car ptlist))                                              ;the first point of Hull.
  43.   (setq pts1 (append ptlist (list pti0)))                               ;add the first point at back of Hull
  44.   (setq pts2 (cdr (append ptlist ptlist (list pti0))))                  ;Construct a loop for the hull
  45.   (setq i 0)                                                            
  46.  
  47.   ;;Find area of encasing rectangle anchored on each edge.
  48.   (repeat (length ptlist)
  49.     (setq pi1 (car   pts1)                                              
  50.           pi2 (cadr  pts1)
  51.           i1x (car   pi1)
  52.           i1y (cadr  pi1)
  53.           i2x (car   pi2)
  54.           i2y (cadr  pi2)
  55.           ix  (- i2x i1x)
  56.           iy  (- i2y i1y)
  57.           il  (distance (list ix iy) '(0.0 0.0))
  58.     )
  59.  
  60.     ;;寻找最左点
  61.     ;;Find a vertex on on first perpendicular line of support
  62.     (while (> (DOTPR ix iy pts2) 0.0)
  63.       (setq pts2 (cdr pts2))
  64.     )
  65.  
  66.     ;;寻找最上点
  67.     ;;Find a vertex on second perpendicular line of suppoer
  68.     (if (= i 0)
  69.       (setq pts3 pts2)
  70.     )
  71.     (while (> (CROSSPR ix iy pts3) 0.0)
  72.       (setq pts3 (cdr pts3))
  73.     )
  74.  
  75.     ;;寻找最右点
  76.     ;;Find a vertex on second perpendicular line of suppoer
  77.     (if (= i 0)
  78.       (setq pts4 pts3)
  79.     )
  80.     (while (< (DOTPR ix iy pts4) 0.0)
  81.       (setq pts4 (cdr pts4))
  82.     )
  83.  
  84.     ;;得出了每边的矩形
  85.     ;;Find distances between parallel and perpendicular lines of support
  86.     (cond
  87.       ((equal i1x i2x 1e-4)                                             ;如果边两点的X值相同
  88.        (setq d1 (- (caar pts3) i1x)                                     ;那么矩形的高就是最上点与边的X的差值
  89.              d2 (- (cadar pts4) (cadar pts2))                           ;矩形的宽就是最左和最右的Y的差值
  90.        )
  91.       )
  92.       ((equal i1y i2y 1e-4)                                             ;如果边两点的Y值相同
  93.        (setq d1 (- (cadar pts3) i1y)                                    ;那么矩形的高就是最上点与边的Y的差值
  94.              d2 (- (caar pts4) (caar pts2))                             ;矩形的宽就是最左和最右的X的差值
  95.        )
  96.       )
  97.  
  98.       (T
  99.        (setq aa (det pi1 pi2 (car pts3)))                               ;否则计算边和最上点构成的面积的二倍(det)
  100.        (setq d1 (/ aa il))                                              ;高就是det值除以边长
  101.        (setq j1 (car pts2))                                             ;最右边点
  102.        (setq j2 (list (- (car j1) iy) (+ (cadr j1) ix)))                ;通过最右边点的垂直边的点
  103.        (setq bb (det j1 j2 (car pts4)))                                 ;最右边点,上面的点和最左边的点
  104.        (setq d2 (/ bb il))                                              ;这三点的det除以边长就是宽
  105.       )
  106.     )
  107.  
  108.     ;;计算矩形的面积,必要时更新最小面积
  109.     ;;Compute area of encasing rectangle anchored on current edge.
  110.     ;;if the area is smaller than the old Minimum area,then update,and record the width,height and five points.
  111.     (setq Ai (abs (* d1 d2)))                                           ;面积就是高和宽的积
  112.     (if (< Ai MinA)                                                     ;如果面积小于先前的最小面积,则记录:
  113.       (setq MinA Ai                                                     ;更新最小面积
  114.             MinH d1                                                     ;最小面积的高
  115.             MinW d2                                                     ;最小面积的宽
  116.             pti1 pi1                                                    ;最小面积的边的第一个端点
  117.             pti2 pi2                                                    ;最小面积的边的第二个端点
  118.             ptj1 (car pts2)                                             ;最右边的点
  119.             ptk1 (car pts3)                                             ;最上面的点
  120.             ptm1 (car pts4)                                             ;最左边的点
  121.       )
  122.     )
  123.     (setq pts1 (cdr pts1))                                              ;检测下一条边
  124.     (setq i (1+ i))                                                     ;计数器加一
  125.   )
  126.  
  127.   ;;according to the result ,draw the Minimum Area Rectangle
  128.   (setq edge (mapcar '- pti2 pti1))                                     ;最小面积的边对应的向量
  129.   (setq VecL (distance edge '(0.0 0.0)))                                ;最小面积的边的长度
  130.   (setq NorH (abs (/ MinH VecL)))                                       ;这边的法线
  131.  
  132.   (setq Norm (list (- (cadr edge)) (car edge)))                         ;这边的垂直向量
  133.   (setq vj12 (mapcar '+ ptj1 Norm))                                     ;通过最右点的垂直向量
  134.   (setq vm12 (mapcar '+ ptm1 Norm))                                     ;通过最左点的垂直向量
  135.   (setq vecH (mapcar '* (list NorH NorH) Norm))                         
  136.  
  137.   (setq rec1 (inters pti1 pti2 ptj1 vj12 nil))                          ;矩形的第一点
  138.   (setq rec4 (inters pti1 pti2 ptm1 vm12 nil))                          ;矩形的第四点
  139.   (setq rec2 (mapcar '+ rec1 vecH))                                     ;矩形的第二点
  140.   (setq rec3 (mapcar '+ rec4 vecH))                                     ;矩形的第三点
  141.   (setq rect (list Rec1 rec2 rec3 rec4))                                ;矩形的点表
  142.   (cons rect MinA)                                                      ;返回这个矩形的点表和最大距离
  143. )
  144.  
  145. ;;;========================================
  146. ;;;求凸壳的直径的程序                      
  147. ;;;参数:逆时针的凸壳 H-------注意逆时针!!!
  148. ;;;返回值: 直径的两个端点和直径 Pair . MaxD
  149. ;;;========================================
  150. (defun Max-distance (H / D M MAXD P PAIR Q U V W)
  151.   (setq Q (cdr (append H H (list (car H)))))                            ;构造一个首尾循环的凸集,且起始点为凸壳的第二点
  152.   (setq MaxD 0.0)                                                       ;初始化最小距离为0
  153.   (foreach U H                                                          ;依次检查凸壳的边
  154.     (setq V (car Q))                                                    ;循环集的第一点
  155.     (setq W (cadr Q))                                                   ;循环集的第二点
  156.     (setq M (mid-pt V W))                                               ;这两点的中点
  157.     (while (> (dot M U V) 0.0)                                          ;如果夹角小于90度(即点积大于0)
  158.       (setq Q (cdr Q))                                                  ;循环集推进
  159.       (setq V (car Q))                                                  ;取下一点
  160.       (setq W (cadr Q))                                                 ;下下一点
  161.       (setq M (mid-pt V W))                                             ;这两点的中点
  162.     )
  163.     (setq D (distance U V))                                             ;计算这时的最大距离
  164.     (if (> D MaxD)                                                      ;如果大于前面的最大距离
  165.       (setq MaxD D                                                      ;就替换前面的最大距离
  166.             Pair (list U V)                                             ;并记录这对点
  167.       )
  168.     )
  169.   )
  170.   (cons Pair MaxD)                                                      ;返回这对点和最大距离
  171. )
  172.  
  173. ;;;==================
  174. ;;;Graham扫描法求凸包
  175. ;;;==================
  176. (defun Graham-scan (ptlist / hullpt maxXpt sortPt P Q)
  177.   (if (< (length ptlist) 3)                                             ;3点以下
  178.     ptlist                                                              ;是本集合
  179.     (progn
  180.       (setq maxXpt (assoc (apply 'max (mapcar 'car ptlist)) ptlist))    ;最右边的点
  181.       (setq sortPt (sort-by-angle-distance ptlist maxXpt))              ;分类点集
  182.       (setq hullPt (list (cadr sortPt) maxXpt))                         ;开始的两点      
  183.       (foreach n (cddr sortPt)                                          ;从第3点开始
  184.         (setq hullPt (cons n HullPt))                                   ;把Pi加入到凸集
  185.         (setq P (cadr hullPt))                                          ;Pi-1
  186.         (setq Q (caddr hullPt))                                         ;Pi-2
  187.         (while (and q (> (det n P Q) -1e-6))                            ;如果左转
  188.           (setq hullPt (cons n (cddr hullPt)))                          ;删除Pi-1点
  189.           (setq P (cadr hullPt))                                        ;得到新的Pi-1点
  190.           (setq Q (caddr hullPt))                                       ;得到新的Pi-2点
  191.         )
  192.       )
  193.       (reverse hullpt)                                                  ;返回凸集
  194.     )
  195.   )
  196. )
  197.  
  198. ;;;中点函数
  199. (defun mid-pt (p1 p2)
  200.   (list
  201.     (* (+ (car p1) (car p2)) 0.5)
  202.     (* (+ (cadr p1) (cadr p2)) 0.5)
  203.   )
  204. )
  205.  
  206. ;;;以某点为基点,按照角度和距离分类点集
  207. (defun sort-by-angle-distance (ptlist pt /)
  208.   (vl-sort
  209.     ptlist
  210.     (function
  211.       (lambda (e1 e2 / ang1 ang2)
  212.         (setq ang1 (angle pt e1))
  213.         (setq ang2 (angle pt e2))
  214.         (if (= ang1 ang2)
  215.           (< (distance pt e1) (distance pt e2))
  216.           (< ang1 ang2)
  217.         )
  218.       )
  219.     )
  220.   )
  221. )
  222.  
  223.  
  224. ;;;点积= x1*x2 + y1*y2
  225. (defun DOTPR (ix iy pts / pt1 pt2)
  226.   (setq pt1 (car pts))
  227.   (setq pt2 (cadr pts))
  228.   (+ (* ix (- (car pt2) (car pt1)))
  229.      (* iy (- (cadr pt2) (cadr pt1)))
  230.   )
  231. )
  232.  
  233. ;;;叉积= x1*y2 - x2*y1
  234. (defun CROSSPR (ix iy pts / pt1 pt2)
  235.   (setq pt1 (car pts))
  236.   (setq pt2 (cadr pts))
  237.   (- (* ix (- (cadr pt2) (cadr pt1)))
  238.      (* iy (- (car pt2) (car pt1)))
  239.   )
  240. )
  241.  
  242. ;;;定义三点的行列式,即三点之倍面积
  243. (defun det (p1 p2 p3 / x2 y2)
  244.   (setq x2 (car p2)
  245.         y2 (cadr p2)
  246.   )
  247.   (- (* (- x2 (car p3)) (- y2 (cadr p1)))
  248.      (* (- x2 (car p1)) (- y2 (cadr p3)))
  249.   )
  250. )
  251.  
  252. ;;;定义向量的点积函数
  253. (defun dot (p1 p2 p3 / x1 y1)
  254.   (setq x1 (car p1)
  255.         y1 (cadr p1)
  256.   )
  257.   (+ (* (- (car p2) x1) (- (car p3) x1))
  258.      (* (- (cadr p2) y1) (- (cadr p3) y1))
  259.   )
  260. )
  261. ;;;取点函数2
  262. (defun getpt (ss n / i s a b c d e)
  263.   (setq i 0)
  264.   (if ss
  265.     (repeat (sslength ss)
  266.       (setq a (ssname ss i))
  267.       (setq b (entget a))
  268.       (setq e (cdr (assoc 0 b)))
  269.       (cond
  270.         ( (= e "LWPOLYLINE")
  271.           (setq c (get-pline-vertexs a n))
  272.           (setq s (append c s))
  273.         )
  274.  
  275.         ( (= e "LINE")
  276.           (setq c (cdr (assoc 10 b)))
  277.           (setq d (cdr (assoc 11 b)))
  278.           (setq c (list (car c) (cadr c)))
  279.           (setq d (list (car d) (cadr d)))
  280.           (setq s (cons c s))
  281.           (setq s (cons d s))
  282.         )
  283.         ( (= e "POINT")
  284.           (setq c (cdr (assoc 10 b)))
  285.           (setq c (list (car c) (cadr c)))
  286.           (setq s (cons c s))
  287.         )
  288.         (t 
  289.           (setq c (get-spline-vertexs a n))
  290.           (setq s (append c s))
  291.         )
  292.       )
  293.       (setq i (1+ i))
  294.     )
  295.   )
  296.   s
  297. )
  298.  
  299. ;;取得多边形顶点
  300. (defun get-LWpolyline-vertexs (DXF / lst)
  301.   (foreach n DXF
  302.     (if (= (car n) 10)
  303.       (setq lst (cons (cdr n) lst))
  304.     )
  305.   )
  306.   (reverse lst)
  307. )
  308.  
  309. (defun get-3dpolyline-vertexs ( ent / p )
  310.   (if (and (setq ent (entnext ent)) (setq p (cdr (assoc 10 (entget ent)))))
  311.     (cons p (get-3dpolyline-vertexs ent))
  312.   )
  313.   p
  314. )
  315.  
  316. ;;;取得样条曲线的点
  317. (defun get-spline-vertexs (ent n / DIST ENDPAR I LEN OBJ PT PTS SEG)
  318.   (setq obj (vlax-ename->vla-object ent))
  319.   (setq endpar  (vlax-curve-getEndParam obj))
  320.   (setq len (vlax-curve-getDistAtParam obj endpar))
  321.   (setq seg (/ len n))
  322.   (setq dist 0)
  323.   (while (< dist len)   
  324.     (setq pt (vlax-curve-getPointAtDist obj dist))
  325.     (setq pts (cons pt pts))
  326.     (setq dist (+ seg dist))   
  327.   )
  328.   (if (= (vla-get-closed obj) :vlax-false)
  329.     (setq pt (vlax-curve-getEndPoint obj)
  330.           pts (cons pt pts)
  331.     )
  332.   )
  333.   (reverse pts)
  334. )
  335.  
  336. ;;;取得含有圆弧的多段线的点
  337. (defun get-pline-vertexs (ent n / BLG DIST ENDPAR I L1 L2 L3 LI OBJ PT PTS VEXNUM)
  338.   (setq obj (vlax-ename->vla-object ent))
  339.   (setq endpar (vlax-curve-getEndParam obj))
  340.   (setq vexNum (fix endPar))
  341.   (setq pts nil)
  342.   (setq i 0)
  343.   (repeat vexNum
  344.     (setq pt (vlax-curve-getPointAtParam obj i))
  345.     (setq pts (cons pt pts))
  346.     (setq blg (vla-getbulge obj i))
  347.     (if (/= blg 0.0)
  348.       (progn
  349.         (setq l1 (vlax-curve-getDistAtParam obj i))
  350.         (setq l2 (vlax-curve-getDistAtParam obj (1+ i)))
  351.         (setq l3 (- l2 l1))
  352.         (setq li (/ l3 n))
  353.         (setq dist l1)
  354.         (repeat (1- n)
  355.           (setq dist (+ dist li))
  356.           (setq pt (vlax-curve-getPointAtDist obj dist))
  357.           (setq pts (cons pt pts))
  358.         )
  359.       )
  360.     )
  361.     (setq i (1+ i))
  362.   )
  363.   (if (= (vla-get-closed obj) :vlax-false)
  364.     (setq pt (vlax-curve-getEndPoint obj)
  365.           pts (cons pt pts)
  366.     )
  367.   )
  368.   pts
  369. )
  370.  
  371. ;;;绘制多段线
  372. (defun Make-Poly (pp / x)
  373.   (entmake                                                              ;画凸包
  374.     (append
  375.       '((0 . "LWPOLYLINE")
  376.         (100 . "AcDbEntity")
  377.         (100 . "AcDbPolyline")
  378.        )
  379.       (list (cons 90 (length pp)))                                      ;顶点个数
  380.       (mapcar
  381.         (function (lambda (x) (cons 10 x)))
  382.         pp
  383.       )                                                                 ;多段线顶点
  384.       (list (cons 70 1))                                                ;闭合的
  385.       (list (cons 62 1))                                                ;红色的
  386.     )
  387.   )
  388. )

评论

点集的直径和最小包围矩形 — 一条评论

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注