【越飞越高讲堂12】最小包围盒和最大距离点对 - AutoLISP/Visual LISP 编程

admin

最小包围盒和最大距离点对在实际生活中有很多运用。
譬如，最小包围盒指的是能包围一个多边形或者一些多边形的最小面积（可能周长）的这样的一个矩形，这个帖子我只是讨论二维的，所以说是矩形，
求得最小包围和就能合理利用材料，选择合适的截面。
在这个算法中，采用了“游标卡尺”的算法，所以运行速度很快，在求出一个形状的凸包之后，就能迅速得到最小矩形。O(n *log(n))的时间复杂度，
大量运算成为可能。
对于CAD本身来说，用getboundingBox得出的结果仅仅是平行与WCS的坐标系统的矩形，而且有时候很不准确，特别对于spline来说，有时候相差很大。
然而，WCS下的boundingBox一般来说不是最小矩形，如果要得到最小的，必须要旋转这个图形很多次才能得到结果，但这时一个费事费力的事情。
所以，我们有必要考虑一种算法。
这些天通过几个比较晚的晚上，终于得出了算法。为了增加讨论效果，我在这里先贴出arx文件，适用于2004-2006平台。至于LISP，以后在贴出。
速度之快，几乎是0秒。

用法： 先输入命令 Hull,
然后输入分弧精度---指的是对样条曲线或者弧形的分段数目,取值在100-2000比较合理,太大容易引起问题且也不能有效提高精度。
选择图形中的多边形或者样条曲线，然后你就可以看到结果了。最远距离用红线标出。

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

复制代码

admin

1. ;;;========================================
2. ;;;求凸壳的直径的程序                     
3. ;;;参数：逆时针的凸壳 H-------注意逆时针!!!
4. ;;;返回值: 直径的两个端点和直径 Pair . MaxD
5. ;;;========================================
6. (defun Max-distance (H / D M MAXD P PAIR Q U V W)
7.   (setq Q (cdr (append H H (list (car H))))) ;构造一个首尾循环的凸集,且起始点为凸壳的第二点
8.   (setq MaxD 0.0)    ;初始化最小距离为0
9.   (foreach U H     ;依次检查凸壳的边
10.     (setq V (car Q))    ;循环集的第一点
11.     (setq W (cadr Q))    ;循环集的第二点
12.     (setq M (mid-pt V W))   ;这两点的中点
13.     (while (> (dot M U V) 0.0)   ;如果夹角小于90度（即点积大于0）
14.       (setq Q (cdr Q))    ;循环集推进
15.       (setq V (car Q))    ;取下一点
16.       (setq W (cadr Q))    ;下下一点
17.       (setq M (mid-pt V W))   ;这两点的中点
18.     )
19.     (setq D (distance U V))   ;计算这时的最大距离
20.     (if (> D MaxD)    ;如果大于前面的最大距离
21.       (setq MaxD D    ;就替换前面的最大距离
22.      Pair (list U V)   ;并记录这对点
23.       )
24.     )
25.   )
26.   (cons Pair MaxD)    ;返回这对点和最大距离
27. )

复制代码