[原创]：递归、分治和分类以及最小距离
给定平面上的一个数量为n的点集，如何能有效地找出距离最近的点对呢？（在实际中有着应用，而且给出的是点对的集合，也就是在误差范围内的所有点对都找出来）。这个问题很容易理解，似乎也不难解决。我们只要将每一点与其他n-1个点的距离算出，找出达到最小距离的两个点即可。然而，这样做效率太低，需要O(n^2)的计算时间。当数量规模较小时，尚能解决，但一旦规模达到万级以上，其速度之慢，时间之长，无法令人忍受。下面我给出这个问题的一个θ(nlogn)算法。
在这个算法中，我利用了递归、分治和分类思想。
递归算法是自身调用自身函数的一种算法，例如求阶乘：
(defun jc(n) (if (= n 0) 1 (* n (jc (1- n)))))
分治算法是对于一个规模为n的问题，把它分解成为K个规模较小的子问题，这些子问题互相独立，且结构与原来问题的结构相同。在解这些子问题时，又对于每一个子问题进行进一步的分解，直到某个阈值为止。递归地解决这些子问题，再把各个子问题的解合并起来，就得到原来问题的解。因此递归一般和分治联系在一起。
对于求最小距离的解，我参考了一些资料，他们给出的大多是C++ 程序，我看不懂，只好按照算法和思想来设计lisp程序。
闲话少说，先看源程序:加载程序，运行TE
;;程序主段-------
(defun f2 (ptlist / l p1 p2 p3 dd 3pmind 3plist ptlist1 ptlist2 ptlist3 ptlist4
n m midpt mind1 mind2 mindt a b c d Dismin Dnmin nplist mindi)
(setq l (length ptlist))
(cond
( (= l 2);;两点还用说
(list (distance (car ptlist) (cadr ptlist)) (list ptlist))
)
( (= l 3);;三点最小距离直接求解点对
(progn
(setq p1 (car ptlist) p2 (cadr ptlist) p3 (caddr ptlist))
(setq dd
(list (list (distance p1 p2) (list p1 p2))
(list (distance p1 p3) (list p1 p3))
(list (distance p2 p3) (list p2 p3))
)
)
(setq 3pmind (apply 'min (mapcar 'car dd)))
(setq 3plist nil)
(foreach 3name dd
(if (equal (car 3name) 3pmind 1e-8)
(setq 3plist (cons (cadr 3name) 3plist))
)
)
(list 3pmind 3plist)
)
)
( (> l 3)
(progn
(setq n (/ l 2) m (- l n));;分治
(setq ptlist1 (cut ptlist 0 (1- m)))
(setq ptlist2 (cut ptlist m l))
(setq midpt (last ptlist1))
(setq mind1 (f2 ptlist1));;递归左边
(setq mind2 (f2 ptlist2));;递归右边
(setq mindT
(cond
((equal (car mind1) (car mind2) 1e-8)(list (car mind1) (append (cadr mind1) (cadr mind2))))
((< (car mind1) (car mind2)) mind1)
(t mind2)
)
)
(setq mindi (car mindT))
(setq a (- (car midpt) mindi) b (car midpt))
(setq ptlist3 (searchX ptlist1 a b))
(if (/= ptlist3 nil)
(progn
(setq Dismin nil)
(foreach name ptlist3
(setq a (car midpt) b (+ (car midpt) mindi) c (- (cadr name) mindi) d (+ (cadr name) mindi))
(setq ptlist4 (searchXY ptlist2 a b c d))
(if (/= ptlist4 nil)
(setq Dismin (cons (6ptmin ptlist4 name) Dismin))
)
)
(if (= Dismin nil)
mindT
(progn
(setq Dnmin (apply 'min (mapcar 'car Dismin)) nplist nil)
(foreach npname Dismin
(if (equal (car npname) Dnmin 1e-8)
(setq nplist (append (cadr npname) nplist))
)
)
(cond
((equal (car mindT) Dnmin 1e-8) (list Dnmin (append nplist (cadr mindT))))
((< (car mindT) Dnmin) mindT)
(t (list Dnmin nplist))
);;for inest cond
);;for inest if-progn
);;for inest if
);;for if-progn
mindT
);;for if
);;for cond-last-progn
);;for cond-last
);;for cond
);;for defun
;;***************
;;以下为测试和附加程序：
(defun C:te (/ olderr en errmsg oldmode oce sl ss t0 ptlist pp pp1)
;;定义错误函数和预处理
(setvar "errno" 0)
(setq olderr *error*)
(defun *error* (msg)
(setq en (getvar "errno"))
(setq errmsg (strcat "errno=" (itoa en) "\nError:" msg))
(alert errmsg)
(setq *error* olderr)
)
(graphscr)
(setq oldmode (getvar "osmode"))
(setq oce (getvar "cmdecho"))
(setvar "cmdecho" 0)
(command ".ucs" "W")
;;也可以用其他方式取得点集
(setq sl '((0 . "POINT")))
(setq t0 (getvar "TDUSRTIMER"))
(setq ss (ssget sl))
(setq ptlist (getpt ss))
;;分类
(setq t0 (getvar "TDUSRTIMER"))
(setq ptlist (sortx ptlist))
(princ "\n函数排序用时")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ "秒")
;;函数用时估算，以了解函数性能
(setq t0 (getvar "TDUSRTIMER"))
(setq pp1 (f2 ptlist) pp (car (cdr pp1)))
(princ "\n函数查找用时")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ "秒")
(if (= nil pp)
(progn
(alert "不存在有最小距离的一对点!")
(command ".ucs" "p")
(setvar "osmode" oldmode)
(setvar "cmdecho" oce)
(princ)
)
(progn
;;画最短距离的点对集的连线，可能有多条
(setvar "osmode" 0)
(foreach nn pp
(entmake
(append
'((0 . "line")(100 . "AcDbEntity")(100 . "AcDbLine"))
(list (cons 10 (carnn)))
(list (cons 11 (cadr nn)))
(list (cons 62 1))
)
)
)
(command ".ucs" "P")
(setvar "osmode" oldmode)
(setvar "cmdecho" oce)
(princ)
)
)
)
;;取点函数，其中i为点的编号
(defun getpt (ss / i listpp a b c)
(setq i 0 listpp nil )
(if ss
(repeat (sslength ss)
(setq a (ssname ss i))
(setq b (entget a))
(setq c (cdr (assoc 10 b)))
(setq listpp (cons c listpp))
(setq i (1+ i))
)
)
(reverse listpp)
)
;;从J到K的表
(defun cut (ptlist j k / i ptlist1)
(setq i 0 ptlist1 nil)
(foreach n ptlist
(if (and (>= i j) (<= i k) )
(setq ptlist1 (cons n ptlist1))
)
(setq i (1+ i))
)
(reverse ptlist1)
)
;;对X排序
(defun sortX (ptlist)
(mapcar '(lambda (x) (nth x ptlist))
(vl-sort-i ptlist '(lambda (e1 e2) (< (car e1)(car e2))))
)
)
;;在带形区域查找
(defun searchX (ptlist1 x1 x2 / pp n)
(setq pp nil)
(foreach n ptlist1
(if (and (>= (car n) x1)(<= (car n) x2))
(setq pp (cons n pp))
)
)
(reverse pp)
)
;;在矩形区域查找
(defun searchXY (ptlist2 x1 x2 y1 y2 / pp n)
(setq pp nil)
(foreach n ptlist2
(if (and (>= (carn) x1)
(<= (carn) x2)
(>= (cadr n) y1)
(<= (cadr n) y2)
)
(setq pp (cons n pp))
)
)
(reverse pp)
)
;;最多6点最小距离
(defun 6ptmin (ptlist4 pt / 6pmin 6plist)
(setq 6pmin (mapcar '(lambda (x) (distance x pt)) ptlist4))
(setq 6pmin (apply 'min 6pmin) 6plist nil)
(foreach 6name ptlist4
(if (equal (distance 6name pt) 6pmin 1e-8)
(setq 6plist (cons (list pt 6name) 6plist))
)
)
(list 6pmin 6plist)
)
在这个论坛上，比较了一些别人写的代码，（大都是平方级别的，有的甚至更高），发现在时间上还是比平方级以上的要快很多。但是还没有做最优处理，希望大家多多提意见给我。
highflybird@sina.com
2006-11-25,kunming 这是一个很有意思的问题，分区、分治在 Lisp 中对表处理应该结合 Lisp 的特点来灵活运用，针对这个命题，可以用下面的思路还可以提高几十倍的效率
1 如果有重复点，不用说了，那些重复的距离是最小的 0
2 对点排序（由左到右，由下到上）
3 先任意求出一个最小距离MinD作参照
4 用这个MinD对排序后的点表依次向后求满足MinD的每个点的点集
5 逐一处理这个点集，查找是否有更小的 MinD 多谢斑竹指点，对于算法和lisp理解的还不透，还望多多指点。
对这个题的算法吸取斑竹的想法正在改进中，但有一点我觉得可能无需改进：
“1 如果有重复点，不用说了，那些重复的距离是最小的0”
因为我这样考虑的，找到重复点也要花一定时间。
而这个程序对于即使最小距离为0的重复点对都可以找出来。
实际上，如果这个程序稍微改一下，不考虑如果存在多对相同最小距离的点对这种情况的话，要快很多。 highflybird老兄真很厉害！
递归偶还真的不太会啊！以后还请多多指教啊！
推荐你如果有空请帮忙优化一下合并直线，弧线的程序啊！
看看最优化的表处理效果转化到实际工作中来
会不会带来令人惊喜的最新的效率纪录!
因为偶暂时发现对表处理最好的体现就是合并重线的程序了啦！ 不好意思,没看明白eachy提到的算法,提点意见和大家探讨一下,
虽说highflybird 的代码确实有值得改进的地方,不过应该很难提高几十倍的效率
查了一下资料:
highflybird 用的是Shamos 和 Hoey 在 1975 年提出了在 O(n lg n) 时间内寻找平面内最近点对的分治标准算法。该算法准确的说复杂度为 3n lg n
1998年上饶师范专科学校数学系的周玉林和复旦大学计算机科学系的熊鹏荣 发表的论文<求平面点集最近点对的一个改进算法>,将该算法改进为2nlgn
August 22, 2005 ,复旦的Qi Ge, Hai-Tao Wang, and Hong Zhu发表的<求平面点集最近点对的改进算法>将复杂度降为 (3n lg n)/2 您的查询字词都已标明如下：最近点对 算法(点击查询词，可以跳到它在文中首次出现的位置)
如果您想保存该页面，可以添加到搜藏
(百度和网页http://cse.seu.edu.cn/people/jmdeng/resource/algorithm/5-Randomized%20Algorithms/3-Randomized%20Algorithm%20of%20Finding%20the%20Closet%20Pair%20of%20Points.pdf的作者无关，不对其内容负责。百度快照谨为网络故障时之索引，不代表被搜索网站的即时页面。)
--------------------------------------------------------------------------------
求最近点对的随机算法
分治法求解需O(nlogn)时间,而用随机算法求解可达到O(n)时间.
对于给定的平面上的n个点(xi,yi)( i=1,2,…,n)和一个距离δ,
以δ为尺寸可以构造一个(逻辑上的)网格,
该网格以(min{xi},min{yj})为网格的最左下点,
之后以步长δ不断向右,向上伸展为长方形网格;
使得(max{xi},max{yj})也含在其中.
由于所有的点均在长方形网格中,最近点对也必在其中.
对每个尺寸为δ×δ的方格,分别向左向上,向右向上,向左向下,
向右向下各扩展一格,可得4个2δ×2δ的方格(i=1,2,3,4). i
δ2Γ
定理:设已知平面上的某两点的距离为δ.
若最近点对中的一个点落在某个δ×δ方格δΓ中,
则在上述的4个2δ× 2δ方格(i=1,2,3,4)中, i
δ2Γ
至少有一个同时含有最近点对的两个点.
推论:若已知平面上的某两点的距离为δ,
且算法对长方形网格中每一个2δ×2δ方格中的所有点对,
均一一计算其距离,则该算法一定可以找出最近点对.
算法:
1) 在点集S(|S|=n)中随机地取一个子集T,使得|T|= n .
(按随机取样算法的分析,n个中取m个(m≤n)时间为O(n))
2) 对T中的所有点对逐一计算距离,求出最近点对距离δ(T).
(点对数为 n ( n -1)/2,每个点对计算时间为O(1),
故总的计算时间为O( n 2)= O(n).)
3) 以δ(T)为尺寸构造一个(逻辑上的)长方形网格,
(设长方形的横向共有m格,纵向共有r格.)
4) 如果m*r≤cn(c可取2～10,取何值合适见后面的讨论),
则开设一个m×r的指针矩阵P,和一个m×r的计数矩阵Q,
然后逐一扫描S中的点sk:
如果sk落在方格中(i,j分别表示所在方格横向,纵向位置), ij
δΓ
则把该点放入P所指的链表内,Q增1,
直至n个点全部检查为止.
(每个点计算时间为O(1),故总的计算时间为O(n).)
然后扫描Q数组,找到含顶点最多的方格.(时间为O(n)) ij
δΓ
如果中点数≤ij
δΓn,则逐一计算中的点距离(时间为O(n)), ij
δΓ
求出其中的最近点对及新的δ'(如果δ'比δ小的话).
如果中点数>ij
δΓn,则把一分为4,成为4个ij
δΓ
2
δ×
2
δ的方格,
找出其中含点最多的方格,一直拆分到方格中的点数≤n,
求出其中的最近点对及新的δ'(时间复杂度期望值是O(n)).
若m*r>cn,则在素数表中找一个略小于cn的素数p,
然后逐一检查S中的点sk,找到其所在方格(时间为O(n)), ij
δΓ
用散列函数H(i,j,p)把sk散列至长为p的桶中:
桶中的每个槽只对应一个方格(若发生冲突,用冲突处理机制), ij
δΓ
每个槽要记录方格中点的个数,各个点用链表连接起来.
S的检查完成以后,找到含点数最多的方格(时间为O(n)), ij
δΓ
用前述方法计算出该方格中的最近点对和新δ'(如果δ'比δ小的话).
5) 以新δ'为尺寸构造一个(逻辑上的)新网格,
设网格中的小方格数为m'×r',则在这个网格中
恰有(m'-1)×(r'-1)个尺寸为2δ'×2δ'的方格.
由定理知,最近点对必定落在某个之中. δ2Γ
逐一检查S中的点,依次把这些点放入(m'-1)×(r'-1)个指针
所指的链表中(或长度为p的散列表中).(时间复杂度是O(n))
6) 逐一检查含有2个点以上的方格, δ2Γ
对方格中的所有点对进行计算,与当前最小的δ*比较,
若小于δ*,则更新δ*,更新当前的最近点对.
当所有含2个点以上的2δ×2δ方格检查完毕时,
点集中的最近点对也就找到了.根据概率分析,
此项工作的时间复杂度期望值是O(n).
关于c的取值(c可取2～10):
c取值大的好处是需要散列的可能性降低,可直接处理二维表,
且每一个格子中的顶点数相对要少,第6步中工作量减少;
坏处是循环次数增多.c取的小的好处,坏处则反之.
一般来说,若顶点较密集,则c要取的稍大. 很好的一个优化算法研究,收藏一下,后面再提建议. 今天来整理这些老帖。
关于这个地方的算法：
我今天重新更新了算法：
;;; The main routine for The Closest Pair of Points Problem.
;;; Highflybird .Made in China :-)
;;; 2006.12. 1st Edition2012.1. revised Edition.
(defun C:test (/ TOL sl to ss Points T0 pp)
(setq TOL 1e-4) ;input Tolerance
(setq ss (ssget '((0 . "POINT")))) ;Get the points
(if (setq Points (GetPoints ss)) ;Gather the points.
(progn
(setq Points1 Points)
(setq t0 (getvar "TDUSRTIMER"))
(setq Points (Sort-By-X Points)) ;sort by X value
(princ "\nIt takes: ")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ " Seconds for Sorting these points. ")
(setq t0 (getvar "TDUSRTIMER"))
(setq pp (CPP Points TOL)) ;to find the closest pair
(princ "\nIt takes: ")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ " Seconds for finding the closest pair of points. ")
(setq t0 (getvar "TDUSRTIMER"))
(setq pq (eea-cpp Points1)) ;to find the closest pair
(princ "\nIt takes: ")
(princ (* (- (getvar "TDUSRTIMER") t0) 86400))
(princ " Seconds for finding the closest pair of points. ")
(and pp (mapcar 'DrawLine (cadr pp))) ;draw the line(s) between the pair(s).
)
)
)
;;; Use recursion to solve this problem (divide-and-conquer)
(defun CPP (Points TOL / Count p1 p2 p3 S3
CPP_3P 3Pts Pts1 Pts2 Pts3 Pts4 n m
midpointmind1mind2 mindt a b c d
Dismin Dnminnplist CloserDistancepts5 Median_X Y
)
(setq Count (length Points))
(cond
( (> Count 3)
(setq n (/ Count 2))
(setq m (- Count n))
(setq Pts1 (Split Points n)) ; Divide
(setq Pts2 (cdr pts1))
(setq Pts1 (car pts1))
(setq midpoint (car Pts2))
(setq mind1 (CPP Pts1 TOL)) ; Recurse the left
(setq mind2 (CPP Pts2 TOL)) ; Recurse the right
(cond
( (equal (car mind1) (car mind2) TOL)
(setq MindT (list (car mind1) (append (cadr mind1) (cadr mind2))))
)
( (< (car mind1) (car mind2))
(setq MinDT mind1)
)
(t
(setq MinDT mind2)
)
)
(setq CloserDistance (car mindT))
(setq Median_X (car midpoint))
(setq Pts3 (searchX Pts1 (- Median_X CloserDistance TOL) Median_X))
(setq pts5 (SearchX pts2 Median_X (+ Median_X CloserDistance)))
(if (null Pts3)
mindT
(progn
(foreach Point Pts3
(setq Y (cadr Point))
(setq Pts4 (searchY Pts5 (- Y CloserDistance) (+ Y CloserDistance)))
(if Pts4
(setq Dismin (cons (6ptmin Pts4 Point TOL) Dismin))
)
)
(if (null Dismin)
mindT
(progn
(setq Dnmin (apply 'min (mapcar 'car Dismin)))
(foreach npname Dismin
(if (equal (car npname) Dnmin TOL)
(setq nplist (append (cadr npname) nplist))
)
)
(cond
( (equal (car mindT) Dnmin TOL)
(list CloserDistance (append nplist (cadr mindT)))
)
( (< (car mindT) Dnmin) mindT)
(t
(list Dnmin nplist)
)
);for inest cond
);for inest if-progn
);for inest if
);for if-progn
);for if
);for cond-last
( (= Count 3) ; Solve directly if three points
(setq p1 (car Points)
p2 (cadr Points)
p3 (caddr Points)
S3 (list (list (distance p1 p2) (list p1 p2))
(list (distance p2 p3) (list p2 p3))
(list (distance p3 p1) (list p3 p1))
)
)
(setq CPP_3P (apply 'min (mapcar 'car S3)))
(foreach e S3
(if (equal (car e) CPP_3P TOL)
(setq 3Pts (cons (cadr e) 3Pts))
)
)
(list (+ CPP_3P TOL) 3Pts)
)
( (= Count 2) ; List directly if just two points
(list (+ (apply 'distance Points) TOL)
(list Points)
)
)
);for cond
);for defun
;;; Gather the points.
(defun GetPoints (ss / i l a b c)
(if ss
(repeat (setq i (sslength ss))
(setq a (ssname ss (setq i (1- i))))
(setq b (entget a))
(setq c (cdr (assoc 10 b)))
(setq l (cons c l))
)
)
)
;;; Split the points into two in position k.
(defun Split (Points k / Right Left)
(setq Right Points)
(repeat k
(setq Left (cons (car Right) Left))
(setq Right (cdr Right))
)
(cons (reverse Left) Right)
)
;;; Sort points by X value.
(defun Sort-By-X (Points)
(vl-sort Points '(lambda (e1 e2) (< (car e1) (car e2))))
)
;;; Sort points by X value.
(defun Sort-By-Y (Points)
(vl-sort Points '(lambda (e1 e2) (< (cadr e1) (cadr e2))))
)
;;; Search in the Tolerance
(defun SearchX (Pts x1 x2 / pp l n)
(setq l pts)
(while (and l (< (caar l) x1))
(setq l (cdr l))
)
(while (and l (> x2 (caar l)))
(setq pp (cons (car l) pp))
(setq l (cdr l))
)
(reverse pp)
)
;;; Search in the Tolerance Zone
(defun SearchY (pts y1 y2 / pp)
(foreach pt pts
(if (<=y1 (cadr pt) y2)
(setq pp (cons pt pp))
)
)
(reverse pp)
)
;;; There are six points at the most in the Tolerance Zone.
(defun 6ptmin (Pts4 pt TOL / 6pmin 6plist)
(setq 6pmin (mapcar (function (lambda (x) (distance x pt))) Pts4))
(setq 6pmin (apply 'min 6pmin))
(foreach 6name Pts4
(if (equal (distance 6name pt) 6pmin TOL)
(setq 6plist (cons (list pt 6name) 6plist))
)
)
(list (+ 6pmin TOL) 6plist)
)
;;; Draw a line according the coordinates of pair
(defun DrawLine (pair)
(entmakeX
(list
'(0 . "LINE")
(cons 10 (car Pair))
(cons 11 (cadr Pair))
(cons 62 1)
)
)
)
完整的附件见下面：
附件里面还有个俄国的高手ElpanovEvgeniy 的一个算法
同样很精彩。
有点深奥，收下了 本帖最后由 dwg001 于 2013-5-6 23:40 编辑
难能可贵的坚持，赞一个！ 好东西，学习学习先 本帖最后由 炫翔 于 2016-11-20 17:24 编辑
Highflybird 发表于 2013-4-9 16:58
测试 CPP 误差设为4的时候距离为5的两个点也包含了
仔细看了代码 是我弄错了，飞鸟版主算法强悍
大家好，有个问题请教一下，在Highflybird编制的这个代码中，怎么忽略距离小于10mm的距离，或者重合的点。再次感谢。
页:
[1]