【越飞越高讲堂12】最小包围盒和最大距离点对
<p>最小包围盒和最大距离点对在实际生活中有很多运用。</p><p>譬如，最小包围盒指的是能包围一个多边形或者一些多边形的最小面积（可能周长）的这样的一个矩形，这个帖子我只是讨论二维的，所以说是矩形，</p>
<p>求得最小包围和就能合理利用材料，选择合适的截面。</p>
<p>在这个算法中，采用了“游标卡尺”的算法，所以运行速度很快，在求出一个形状的凸包之后，就能迅速得到最小矩形。O(n *log(n))的时间复杂度，</p>
<p>大量运算成为可能。</p>
<p>对于CAD本身来说，用getboundingBox得出的结果仅仅是平行与WCS的坐标系统的矩形，而且有时候很不准确，特别对于spline来说，有时候相差很大。</p>
<p>然而，WCS下的boundingBox一般来说不是最小矩形，如果要得到最小的，必须要旋转这个图形很多次才能得到结果，但这时一个费事费力的事情。</p>
<p>所以，我们有必要考虑一种算法。</p>
<p>这些天通过几个比较晚的晚上，终于得出了算法。为了增加讨论效果，我在这里先贴出arx文件，适用于2004-2006平台。至于LISP，以后在贴出。</p>
<p>速度之快，几乎是0秒。</p>
<p> </p>
<p>用法： 先输入命令 Hull,</p>
<p>然后输入<font style="background-color: rgb(255, 255, 255);" face="Verdana">分弧精度---指的是对样条曲线或者弧形的分段数目,取值在100-2000比较合理,太大容易引起问题且也不能有效提高精度。</font></p>
<p>选择图形中的多边形或者样条曲线，然后你就可以看到结果了。最远距离用红线标出。</p>
<p></p>
<p></p> 在上篇帖子中我已经说明了这个算法的大致思想。下面我贴出源码，具体算法我要强调的一点是：
我构造了一个循环的凸包点集，这是很重要的，因为能极大程度地提高lisp的运算速度，特别是在while,repeat之类的语句中时候，需要大量循环
的时候，如果对点集的读取频繁的时候，就显得极为重要，它只是用到了car,cdr,cons这样的基本函数，而没有用append,nth之类的函数。
我是基于这样的考虑：Lisp是链表结构，不像 C的数组，随即存取速度极其快，读取C的某个元素，经常要从表头开始，读到位置所在为止。
其他的算法说明在注释中解释得已经详尽，如果各位有兴趣，不妨提供建议。
另外我附加了求最远距离的程序，和一个用Graham 扫描法求凸包的程序，这两个程序我写的自认为很简洁。仅仅短短的几行，就能达到很高的效率。
;;;=======================================================
;;;Function : Find the minimum area of encasing rectangle.
;;;Arguments : A CCW HULL
;;;Return: The Four points of Rectangle and its Area
;;;=======================================================
(defun MinAreaRectangle (ptlist / AAAI BB D1
D2 EDGEI I1XI1Y I2X I2Y
IL INF IX IYJ1 J2 MINA
MINHMINWNORHNORMPI1 PI2 PTI0
PTI1PTI2PTJ1PTK1PTM1PTS1PTS2
PTS3PTS4REC1REC2REC3REC4RECT
VECHVECLVJ12VM12
)
(setq INF 1e309)
(setq minA INF) ;Initiating the Minimum area is infinite
(setq pti0 (car ptlist)) ;the first point of Hull.
(setq pts1 (append ptlist (list pti0))) ;add the first point at back of Hull
(setq pts2 (cdr (append ptlist ptlist (list pti0)))) ;Construct a loop for the hull
(setq i 0)
;;Find area of encasing rectangle anchored on each edge.
(repeat (length ptlist)
(setq pi1 (car pts1)
pi2 (cadrpts1)
i1x (car pi1)
i1y (cadrpi1)
i2x (car pi2)
i2y (cadrpi2)
ix(- i2x i1x)
iy(- i2y i1y)
il(distance (list ix iy) '(0.0 0.0))
)
;;寻找最左点
;;Find a vertex on on first perpendicular line of support
(while (> (DOTPR ix iy pts2) 0.0)
(setq pts2 (cdr pts2))
)
;;寻找最上点
;;Find a vertex on second perpendicular line of suppoer
(if (= i 0)
(setq pts3 pts2)
)
(while (> (CROSSPR ix iy pts3) 0.0)
(setq pts3 (cdr pts3))
)
;;寻找最右点
;;Find a vertex on second perpendicular line of suppoer
(if (= i 0)
(setq pts4 pts3)
)
(while (< (DOTPR ix iy pts4) 0.0)
(setq pts4 (cdr pts4))
)
;;得出了每边的矩形
;;Find distances between parallel and perpendicular lines of support
(cond
((equal i1x i2x 1e-4) ;如果边两点的X值相同
(setq d1 (- (caar pts3) i1x) ;那么矩形的高就是最上点与边的X的差值
d2 (- (cadar pts4) (cadar pts2)) ;矩形的宽就是最左和最右的Y的差值
)
)
((equal i1y i2y 1e-4) ;如果边两点的Y值相同
(setq d1 (- (cadar pts3) i1y) ;那么矩形的高就是最上点与边的Y的差值
d2 (- (caar pts4) (caar pts2)) ;矩形的宽就是最左和最右的X的差值
)
)
(T
(setq aa (det pi1 pi2 (car pts3))) ;否则计算边和最上点构成的面积的二倍(det)
(setq d1 (/ aa il)) ;高就是det值除以边长
(setq j1 (car pts2)) ;最右边点
(setq j2 (list (- (car j1) iy) (+ (cadr j1) ix)));通过最右边点的垂直边的点
(setq bb (det j1 j2 (car pts4))) ;最右边点，上面的点和最左边的点
(setq d2 (/ bb il)) ;这三点的det除以边长就是宽
)
)
;;计算矩形的面积，必要时更新最小面积
;;Compute area of encasing rectangle anchored on current edge.
;;if the area is smaller than the old Minimum area,then update,and record the width,height and five points.
(setq Ai (abs (* d1 d2))) ;面积就是高和宽的积
(if (< Ai MinA) ;如果面积小于先前的最小面积，则记录：
(setq MinA Ai ;更新最小面积
MinH d1 ;最小面积的高
MinW d2 ;最小面积的宽
pti1 pi1 ;最小面积的边的第一个端点
pti2 pi2 ;最小面积的边的第二个端点
ptj1 (car pts2) ;最右边的点
ptk1 (car pts3) ;最上面的点
ptm1 (car pts4) ;最左边的点
)
)
(setq pts1 (cdr pts1)) ;检测下一条边
(setq i (1+ i)) ;计数器加一
)
;;according to the result ,draw the Minimum Area Rectangle
(setq edge (mapcar '- pti2 pti1)) ;最小面积的边对应的向量
(setq VecL (distance edge '(0.0 0.0))) ;最小面积的边的长度
(setq NorH (abs (/ MinH VecL))) ;这边的法线
(setq Norm (list (- (cadr edge)) (car edge))) ;这边的垂直向量
(setq vj12 (mapcar '+ ptj1 Norm)) ;通过最右点的垂直向量
(setq vm12 (mapcar '+ ptm1 Norm)) ;通过最左点的垂直向量
(setq vecH (mapcar '* (list NorH NorH) Norm))
(setq rec1 (inters pti1 pti2 ptj1 vj12 nil)) ;矩形的第一点
(setq rec4 (inters pti1 pti2 ptm1 vm12 nil)) ;矩形的第四点
(setq rec2 (mapcar '+ rec1 vecH)) ;矩形的第二点
(setq rec3 (mapcar '+ rec4 vecH)) ;矩形的第三点
(setq rect (list Rec1 rec2 rec3 rec4)) ;矩形的点表
(cons rect MinA) ;返回这个矩形的点表和最大距离
)
另外，我已经编写好了arx代码，基于这里只是讨论lisp，故不贴出源文件。
需要附加说明的是:
如你需要转载，请你说明帖子的来源和原创者。这是尊重我的劳动的起码要求，也是这个论坛的共享“人人为我，我为人人”的体现，而不是一味的索求源码。强烈鄙视拿这些代码卖钱的行为！
本帖最后由 作者 于 2010-6-2 22:38:13 编辑
首先我谈谈最远距离的思路：如果求出了凸壳，那么我们假设这凸壳是逆时针的，从凸壳的每一条边开始，寻找它的最大距离点。
表面上看起来是一个O（n^2)的算法，但实际不是；我引用网上的一段说明：
旋转卡壳可以用于求凸包的直径、宽度，两个不相交凸包间的最大距离和最小距离等。虽然算法的思想不难理解，但是实现起来真的很容易让人“卡壳”。
拿凸包直径(也就是凸包上最远的两点的距离）为例，原始的算法是这样子：
1，计算凸壳的最左点ymin 和最右点ymax.
2，构造两条平行线通过 ymin和ymax. 既然这已经是一条对踵点，计算这两点的距离,设为初始的最大值.
3，旋转这两条平行线直到与凸壳的一边重合.
4，步进新的对锺点， 计算新的距离，并跟前面的最大值比较，,如大于则更新最大值.
5，重复步骤3和 4，直到对踵点又达到 (ymin,ymax) .
6，输出最大距离和对踵点。
http://cgm.cs.mcgill.ca/~orm/rotcal.frame.html
直接按照这个描述可以实现旋转卡壳算法，但是代码肯定相当冗长。逆向思考，如果qa，qb是凸包上最远两点，必然可以分别过qa，qb画出一对平行线。通过旋转这对平行线，我们可以让它和凸包上的一条边重合,如图中蓝色直线，可以注意到，qa是凸包上离p和qb所在直线最远的点。于是我们的思路就是枚举凸包上的所有边，对每一条边找出凸包上离该边最远的顶点，计算这个顶点到该边两个端点的距离，并记录最大的值。直观上这是一个O(n2)的算法，和直接枚举任意两个顶点一样了。但是注意到当我们逆时针枚举边的时候，最远点的变化也是逆时针的，这样就可以不用从头计算最远点，而可以紧接着上一次的最远点继续计算(详细的证明可以参见上面链接中的论文)。于是我们得到了O(n)的算法。
实际实现它的代码真的很简洁。下面是程序的主要代码,需要说明的是我的算法是对上述算法的改进，来源于周培德的一本书，并在他的算法的基础上继续进行了改进，并使之适合于LISP语言。
;;;========================================
;;;求凸壳的直径的程序
;;;参数：逆时针的凸壳 H-------注意逆时针!!!
;;;返回值: 直径的两个端点和直径 Pair . MaxD
;;;========================================
(defun Max-distance (H / D M MAXD P PAIR Q U V W)
(setq Q (cdr (append H H (list (car H))))) ;构造一个首尾循环的凸集,且起始点为凸壳的第二点
(setq MaxD 0.0) ;初始化最小距离为0
(foreach U H ;依次检查凸壳的边
(setq V (car Q)) ;循环集的第一点
(setq W (cadr Q)) ;循环集的第二点
(setq M (mid-pt V W)) ;这两点的中点
(while (> (dot M U V) 0.0) ;如果夹角小于90度（即点积大于0）
(setq Q (cdr Q)) ;循环集推进
(setq V (car Q)) ;取下一点
(setq W (cadr Q)) ;下下一点
(setq M (mid-pt V W)) ;这两点的中点
)
(setq D (distance U V)) ;计算这时的最大距离
(if (> D MaxD) ;如果大于前面的最大距离
(setq MaxD D ;就替换前面的最大距离
Pair (list U V) ;并记录这对点
)
)
)
(cons Pair MaxD) ;返回这对点和最大距离
)
而对于最小面积的包围矩形和最小周长的包围矩形也采取了类似的游标卡尺的算法，它们的算法概念如下：
1，从凸壳的第一边，寻找对于这边的最左点，最右点，最上点，得到初始的最小面积
2，从第二条边奥开始，对于寻找最左边点，只要从上条边计算的点开始计算就可以了，对于最右点，和最上点也是如此。
3，计算这条边的包围矩形面积，如果小于前面的最小面积，则更新最小面积，并记录这条边和最左点，最右点，最上点。
4，直到凸壳的最后一条边。
5，根据最小面积时候记录的五个点，得出凸壳的最小面积包围矩形。
至于主要的源代码，我下次回帖则贴出。
非常好,谢谢分享 已经有了这样的程序。。。。。 <p>不错，很实用的工具</p>