2011-07-14 29 views
25

我在Mathematica中实现了一个quadtree。我是用Mathematica等函数式编程语言进行编码的新手,我想知道是否可以通过更好地使用模式来改进它或使其更加紧凑。在Mathematica中实现四叉树

(据我所知,我也许可以通过修剪未使用的节点优化树,并有可能会像对空间分解KD树更好的数据结构。)

而且,我仍然不舒服复制的想法每次添加新点时,整个树/表达式。但我的理解是,对整个表达式进行操作而不修改这些部分是功能性编程方式。我很感谢在这方面的任何澄清。

MV

守则

ClearAll[qtMakeNode, qtInsert, insideBox, qtDraw, splitBox, isLeaf, qtbb, qtpt]; 

(* create a quadtree node *) 
qtMakeNode[{{xmin_,ymin_}, {xmax_, ymax_}}] := 
{{}, {}, {}, {}, qtbb[{xmin, ymin}, {xmax, ymax}], {}} 

(* is pt inside box? *) 
insideBox[pt_, bb_] := If[(pt[[1]] <= bb[[2, 1]]) && (pt[[1]] >= bb[[1, 1]]) && 
    (pt[[2]] <= bb[[2, 2]]) && (pt[[2]] >= bb[[1, 2]]), 
    True, False] 

(* split bounding box into 4 children *) 
splitBox[{{xmin_,ymin_}, {xmax_, ymax_}}] := { 
{{xmin, (ymin+ymax)/2}, {(xmin+xmax)/2, ymax}}, 
{{xmin, ymin},{(xmin+xmax)/2,(ymin+ymax)/2}}, 
{{(xmin+xmax)/2, ymin},{xmax, (ymin+ymax)/2}}, 
{{(xmin+xmax)/2, (ymin+ymax)/2},{xmax, ymax}} 
} 

(* is node a leaf? *) 
isLeaf[qt_] := If[ And @@((# == {})& /@ Join[qt[[1;;4]], {List @@ qt[[6]]}]),True, False] 

(*--- insert methods ---*) 

(* qtInsert #1 - return input if pt is out of bounds *) 
qtInsert[qtree_, pt_] /; !insideBox[pt, List @@ qtree[[5]]]:= qtree 

(* qtInsert #2 - if leaf, just add pt to node *) 
qtInsert[qtree_, pt_] /; isLeaf[qtree] := 
{qtree[[1]],qtree[[2]],qtree[[3]],qtree[[4]],qtree[[5]], qtpt @@ pt} 

(* qtInsert #3 - recursively insert pt *) 
qtInsert[qtree_, pt_] := 
    Module[{cNodes, currPt}, 
    cNodes = qtree[[1;;4]]; 
    (* child nodes not created? *) 
    If[And @@ ((# == {})& /@ cNodes), 
    (* compute child node bounds *) 
    (* create child nodes with above bounds*) 
    cNodes = qtMakeNode[#]& /@ splitBox[List @@ qtree[[5]]]; 
    ]; 
    (* move curr node pt (if not empty) into child *) 
    currPt = List @@ qtree[[6]]; 
    If[currPt != {}, 
    cNodes = qtInsert[#, currPt]& /@ cNodes; 
    ]; 
(* insert new pt into child *) 
cNodes = qtInsert[#, pt]& /@ cNodes; 
(* return new quadtree *) 
{cNodes[[1]],cNodes[[2]], cNodes[[3]], cNodes[[4]], qtree[[5]], {}} 
] 

(* draw quadtree *) 
qtDraw[qt_] := Module[{pts, bboxes}, 
    pts = Cases[qt, _qtpt, Infinity] /. qtpt :> List; 
    bboxes = Cases[qt, _qtbb, Infinity] /. qtbb :> List; 
    Graphics[{ 
    EdgeForm[Black],Hue[0.2], Map[Disk[#, 0.01]&, pts], 
    Hue[0.7],EdgeForm[Red], FaceForm[],(Rectangle @@ #) & /@ bboxes 
    }, 
    Frame->True 
] 
] 

使用

Clear[qt]; 
len = 50; 
pts = RandomReal[{0, 2}, {len, 2}]; 
qt = qtMakeNode[{{0.0, 0.0}, {2.0, 2.0}}]; 
Do[qt = qtInsert[qt, pts[[i]]], {i, 1, len}] 
qtDraw[qt] 

输出

enter image description here

+1

我没有时间与你的代码的发挥,但在应对,” ......我仍然不满意每次添加新点时复制整个树/表达式的想法...'---您是否正在寻找一种方法为您的函数添加内存?如果是这样,搜索'记住他们的价值的功能'。如果不是你要找的东西,无视,今天晚些时候我会玩。看起来不错。 :) – telefunkenvf14

+0

不,我的意思是,当我添加一个点时,而不是将节点插入到退出的树中,实际上,我们正在生成一个新的树,丢弃旧的树。我只是想知道这种内存管理的含义,当这些表达式真的很大。很难摆脱C++的心态! ;-) 谢谢。 –

+1

如果我正确理解你,你想通过引用传递,因为你不喜欢'qtInsert'中的'Return'行来创建一个新的嵌套列表(可能)巨大的旧树和少量新节点。Mathematica中有一种方法可以用'Attributes [qtInsert] = HoldAll;'来做这样的事情,其缺点是'qtInsert'的所有参数必须是变量,而不是文字值。 – bbtrb

回答

11

这里是一个更紧凑的版本。它使用与原始版本相同的数据结构。功能splitBoxinsideBox基本上也是一样的(只是写法略有不同)。

而不是逐个添加点,初始框包含开始处的所有点,因此不需要qtInsert例程。在每个递归步骤中,包含多个点的框被分割并且点分布在子框上。这意味着所有具有多个点的节点都是叶子,因此不需要检查。

qtMakeNode[bb_, pts_] := {{}, {}, {}, {}, qtbb @@ bb, pts} 

splitBox[bx_] := splitBox[{min_, max_}] := {min + #, max + #}/2 & /@ 
    Tuples[Transpose[{min, max}]] 


insideBox[pt_, bb_] := bb[[1, 1]] <= pt[[1]] <= bb[[2, 1]] && 
    bb[[1, 2]] <= pt[[2]] <= bb[[2, 2]] 

distribute[qtree_] := Which[ 
    Length[qtree[[6]]] == 1, 
    (* no points in node -> return node unchanged *) 
    qtree, 

    Length[qtree[[6]]] == 1, 
    (* one point in node -> replace head of point with qtpt and return node *) 
    ReplacePart[qtree, 6 -> qtpt @@ qtree[[6, 1]]], 

    Length[qtree[[6]]] > 1, 
    (* multiple points in node -> create sub-nodes and distribute points *) 
    (* apply distribute to sub-nodes *) 
    Module[{spl = splitBox[qtree[[5]]], div, newtreelist}, 
    div = Cases[qtree[[6]], a_ /; insideBox[a, #], 1] & /@ spl; 
    ReplacePart[qtree, 
    Join[Table[i -> distribute[qtMakeNode[spl[[i]], div[[i]]]], {i, 4}], 
     {6 -> {}}]]]] 

例子(使用的qtDraw原始版本):

len = 50; 
pts = RandomReal[{0, 2}, {len, 2}]; 
qt = makeTree[qtMakeNode[{{0.0, 0.0}, {2.0, 2.0}}, pts]]; 
qtDraw[qt] 

结果:

Quadtree example

+0

感谢您的代码 - 很多有趣的用法。但我认为能够逐个添加点是很重要的。例如,考虑粒子逐个进入的2D模拟。 –

+0

@喜欢你的图像是空白的,或者我的机器有一些困难。 –

+0

@belisarius:奇怪,在我的屏幕上它是可见的。当你关注这个链接时,你能看到它吗? http://i.stack.imgur.com/kNyrg.jpg – Heike

3

这可能不是你想要做什么,但最近[ ]可以创建一个内建四叉树结构的NearestFunction []。

35

我认为你的代码不像你期望的那样饥饿。它确实打破和改革名单,但它往往保持大多数的子列表完好无损。

正如其他人所指出的,使用Hold包装器和/或HoldXXX属性仍然可以做得更好,以便模拟逐个引用。

对于硬核的方法来一些相关的数据结构的实现,见

http://library.wolfram.com/infocenter/MathSource/7619/

相关的代码是在笔记本Hemmecke-final.nb(如此命名是因为它实现了由于散光Groebner基算法R. Hemmecke和合着者)。

我用Hold ...属性刺穿了重新实现,但是我并不擅长这一点,并且在代码刺穿了我的时候放弃了它(错过了,但是杀死了我的Mathematica会话)。相反,我有一个使用未记录的“原始”Mathematica数据类型的实现,该类型是惰性的,因此适用于逐引用行为。

由于通用的Mathematica数据结构是“expr”,所以这个结构被称为“expr bag”。它就像一个List,但是(1)它可以在一端增长(虽然不缩小),(2)与其他原始表达类型(例如版本8中的图)一样,它具有可以通过提供的功能访问和/或更改的组件(一个API,可以这么说)。它的基本“元素”是惰性的,因为它们可以引用ANY表达式(包括包本身),并且可以用我将在下面指出的方式进行操作。

上面的第一项提供了实施Sow/Reap的底层技术。这是第二个将在下面的代码感兴趣。最后,我会沿着解释数据结构的方向加入一些评论,因为没有正式的文档。

我把代码或多或少地保留在与原文相同的样式中,特别是它仍然是一个在线版本(也就是说,元素不需要全部进入,但可以单独添加)。改了几个名字。由类似于

节点的基本结构如果存在子节点,则该值字段为空(边界框,值,零个或四个子节点)

。 box和value字段用通常的Mathematica List表达式表示,尽管使用专用的头部并且使其更类似于C结构风格也许是有意义的。我的确在做一些类似的事情来命名各种访问/设置功能。

需要注意的是,这种原始数据类型消耗的内存开销比例如一个列表。所以我下面的变体会使用比原来发布的代码更多的内存。不是渐近地更多,只是一个恒定的因素。在访问或设置元素值方面,它也需要一个不变的因子,比如一个可比较的C结构。所以这不是一个神奇的子弹,只是一种行为不应该给予渐近式惊喜的数据类型。


AppendTo[$ContextPath, "Internal`"]; 

makeQuadTreeNode[bounds_] := Bag[{bounds, {}, {}}] 

(*is pt inside box?*) 

insideBox[pt_, box_] := 
And @@ Thread[box[[1]] <= (List @@ pt) <= box[[2]]] 

(*split bounding box into 4 children*) 

splitBox[{{xmin_, ymin_}, {xmax_, ymax_}}] := 
Map[makeQuadTreeNode, {{{xmin, (ymin + ymax)/2}, {(xmin + xmax)/2, 
    ymax}}, {{xmin, 
    ymin}, {(xmin + xmax)/2, (ymin + ymax)/2}}, {{(xmin + xmax)/2, 
    ymin}, {xmax, (ymin + ymax)/2}}, {{(xmin + xmax)/ 
     2, (ymin + ymax)/2}, {xmax, ymax}}}] 

bounds[qt_] := BagPart[qt, 1] 
value[qt_] := BagPart[qt, 2] 
children[qt_] := BagPart[qt, 3] 

isLeaf[qt_] := value[qt] =!= {} 
isSplit[qt_] := children[qt] =!= {} 
emptyNode[qt_] := ! isLeaf[qt] && ! isSplit[qt] 

(*qtInsert #1-return input if pt is out of bounds*) 

qtInsert[qtree_, pt_] /; ! insideBox[pt, bounds[qtree]] := qtree 

(*qtInsert #2-empty node (no value,no children)*) 

qtInsert[qtree_, pt_] /; emptyNode[qtree] := value[qtree] = pt 

(*qtInsert #2-currently a leaf (has a value and no children)*) 

qtInsert[qtree_, pt_] /; isLeaf[qtree] := Module[ 
    {kids = splitBox[bounds[qtree]], currval = value[qtree]}, 
    value[qtree] = {}; 
    children[qtree] = kids; 
    Map[(qtInsert[#, currval]; qtInsert[#, pt]) &, kids]; 
    ] 

(*qtInsert #4-not a leaf and has children*) 

qtInsert[qtree_, pt_] := Map[qtInsert[#, pt] &, children[qtree]]; 

getBoxes[ee_Bag] := 
Join[{bounds[ee]}, Flatten[Map[getBoxes, children[ee]], 1]] 
getPoints[ee_Bag] := 
Join[{value[ee]}, Flatten[Map[getPoints, children[ee]], 1]] 

qtDraw[qt_] := Module[ 
    {pts, bboxes}, 
    pts = getPoints[qt] /. {} :> Sequence[]; 
    bboxes = getBoxes[qt]; 
    Graphics[{EdgeForm[Black], Hue[0.2], Map[Disk[#, 0.01] &, pts], 
    Hue[0.7], EdgeForm[Red], 
    FaceForm[], (Rectangle @@ #) & /@ bboxes}, Frame -> True]] 

下面是一个例子。我会注意到缩放是合理的。也许O(n log(n))左右。肯定比O(n^2)好。

len = 4000; 
pts = RandomReal[{0, 2}, {len, 2}]; 
qt = makeQuadTreeNode[{{0.0, 0.0}, {2.0, 2.0}}]; 
Timing[Do[qtInsert[qt, pts[[i]]], {i, 1, len}]] 

{1.6, Null} 

一般EXPR袋笔记。这些都是旧的,所以我并不认为这一切仍然按照指示进行。

这些函数存在于内部环境中。

包 创建一个expr包,可选带有预设元素。

BagPart 获取expr包的部分,类似于普通零件 exprs。也可用于lhs,例如重置一个值。

StuffBag 将元素追加到包的结尾。

我们也有BagLength。用于遍历一个包。

这些功能是有两个原因extrememly有用。

首先,这是为了让一个可扩展的表 数学的好方法。

其次,袋的内容进行评估,但然后置于 原料EXPR,因此被屏蔽。因此,人们可以利用这些作为 “指针”(在C意义上的),而不是作为对象,这 不需要保留等。下面是一些例子:

a = {1,2,a} (* gives infinite recursion *) 

如果我们改用袋我们得到了一个自我 - 参照结构。

In[1]:= AppendTo[$ContextPath, "Internal`"]; 

In[2]:= a = Bag[{1,2,a}] 
Out[2]= Bag[<3>] 

In[3]:= expr1 = BagPart[a, All] 
Out[3]= {1, 2, Bag[<3>]} 

In[4]:= expr2 = BagPart[BagPart[a, 3], All] 
Out[4]= {1, 2, Bag[<3>]} 

In[5]:= expr1 === expr2 
Out[5]= True 

这很难在Mathematica中以任何其他方式效仿。 人们需要使用稀疏表(哈希)在一些 不是非常透明的方式。

这是一个相关的例子,没有完全调试。我们基本上 实现链表,其中一个可以修改破坏性尾巴 ,更换子列表等

tail[ll_] := BagPart[ll,2] 
settail[ll_, ll2_] := BagPart[ll,2] = ll2 
contents[ll_] := BagPart[ll,1] 
setcontents[ll_, elem_] := BagPart[ll,1] = elem 

createlinkedlist[elems__] := Module[ 
    {result, elist={elems}, prev, el}, 
    result = Bag[{elist[[1]],Bag[]}]; 
    prev = result; 
    Do [el = Bag[{elist[[j]],Bag[]}]; 
     settail[prev, el]; 
     prev = el, 
     {j,2,Length[elist]}]; 
    result 
    ] 

In[18]:= tt = createlinkedlist[vv,ww,xx] 
Out[18]= Bag[<2>] 

In[20]:= BagPart[tt,All] 
Out[20]= {vv, Bag[<2>]} 

所以TT是一个链表,第一个元素是VV,接下来就是 本身就是一个链表等。我没有使用Lisp 术语(car/cdr等),因为我无法回想起Lisp的列表操作是否具有破坏性。但是 你得到了一般想法。

按照类似的思路,我已经使用EXPR袋来实现二进制 树。这是有用的,因为我们可以在 固定的时间做破坏性的变化(假设我们已经有了关于插入/缺失的点 一个“把手”),而且expr的“原始”自然 袋意味着我们完全避免了数学的无限评价 语义。

另一个应用程序,也许。

Pointer = Internal`Bag 
Contents[aa_Pointer, j_Integer] /;0<j<=Internal`BagLength[aa] := 
    Internal`BagPart[aa,j] 
SetContents[aa_Pointer, j_Integer, e_] /; 0<j<=Internal`BagLength[aa] := 
    Internal`BagPart[aa,j] = e 
SetContents[aa_Pointer, j_Integer, e_] /; j>BagLength[aa] := 
    (Do[Internal`StuffBag[aa,Null], {k,Internal`BagLength[aa]+1,j-1}]; 
    Internal`StuffBag[aa,e]) 

尝试用

a = Bag[{1,2,a,6,t,y,99,Bag[{a,q,3,r,a,5,t}]}] 
expr1 = BagPart[a, All] 
expr2 = BagPart[BagPart[a, 3], All] 

Contents[a, 4] 
SetContents[a, 7, Contents[a,7]+5] 
SetContents[a,11,33] 

丹尼尔Lichtblau 沃尔夫勒姆研究

+0

优秀的信息!似乎是建立列表的非常快速的方式 – acl

+0

+1,我刚刚通读这第一次。好的信息,谢谢。 – rcollyer