Implementing a Quadtree in Mathematica

Question

I have implemented a quadtree in Mathematica. I am new to coding in a functional programming language like Mathematica, and I was wondering if I could improve this or make it more compact by better use of patterns.

(I understand that I could perhaps optimize the tree by pruning unused nodes, and there might be better data structures like k-d trees for spatial decomposition.)

Also, I am still not comfortable with the idea of copying the entire tree/expression every time a new point is added. But my understanding is that operating on the expression as a whole and not modifying the parts is the functional programming way. I'd appreciate any clarification on this aspect.

MV

The Code

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
 ]
]

Usage

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]

Output

Selfreliant

I don't have time to play with your code, but in response to, '...I am still not comfortable with the idea of copying the entire tree/expression every time a new point is added...'---are you looking for a way to add memory to your function(s)? If so, search for 'functions that remember their values'. If not what you're looking for, disregard and I'll play with it later today. Looks cool. :) – Dicrotic

No, what I meant is, when I add a point, instead of nodes being inserted to an exiting tree, in effect, we are generating an entire new tree, discarding the old. I am just wondering about the memory management implications of this, when such expressions are really huge. Hard to get off the C++ mindset! ;-) Thanks. – Selfreliant

If I understand you correctly, you want to pass-by-reference because you don't like the Return line in qtInsert to create a new nested list out of the (possibly) huge old tree and the few new nodes. There is a way in Mathematica to do something like this with Attributes[qtInsert] = HoldAll; with the shortcoming, that all arguments of qtInsert then have to be variables, not literal values. – Kutuzov

Thanks @Kutuzov - I'll need to read up on this. – Selfreliant

Answer 1

Here is a more compact version. It uses the same data structure as the original version. The functions splitBox and insideBox are essentially the same as well (just written in a slightly different way).

Instead of adding points one-by-one, the initial box contains all the points at the beginning so there is no need for the qtInsert routines. In each recursion step, the boxes containing more than one point are split and the points are distributed over the sub-boxes. This means that all nodes with more than one point are leafs so there is no need to check for that either.

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 -> {}}]]]]

Example (using the original version of qtDraw):

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

Result:

Nostrum

Thanks for the code - lots of interesting usage. But I do think that the ability to add points one by one is important. For instance, consider a 2D simulation in which particles are coming in one by one. – Selfreliant

@Nostrum Either your image is blank or my machine has some difficulty with it. – Anemone

@belisarius: strange, on my screen it is visible. Can you see it when you following this link? i.stack.imgur.com/kNyrg.jpg – Nostrum

@Nostrum Don-t mind. When I went to the imgur.com home page I got Emergency maintenance - We're currently expirencing issues with uploading and viewing some images, and are working on correcting them now. Please check back soon for updates. So probably you are looking at a cached image on your machine. – Anemone

Answer 2

I think your code is not so memory hungry as you might expect. It does break and reform lists, but it tends to keep most sublists intact.

As others remarked, it might be possible to do better still using Hold wrappers and/or HoldXXX attributes, so as to emulate call-by-reference.

For a hard core approach to some related data structure implementations, see

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

The relevant code is in the notebook Hemmecke-final.nb (so named because it implements toric Groebner basis algorithm due to R. Hemmecke and coauthors).

I took a stab at reimplementing using Hold... attributes, but I'm not terribly good at that and gave it up when the code took a stab back at me (missed, but killed my Mathematica session). So instead I have an implementation that uses an undocumented "raw" Mathematica data type, one that is inert and thus amenable to call-by-reference behavior.

The structure in question is called an "expr bag" because the generic Mathematica data structure is the "expr". It is like a List but (1) It can grow at one end (though not shrink) and (2) like other raw expression types (e.g. graphs in version 8) it has components that can be accessed and/or changed via provided functions (an API, so to speak). Its underlying "elements" are inert in the sense that they can reference ANY expr (including the bag itself) and can be manipulated in ways I'll indicate below.

The first item above provides the underlying technology for the implementation of Sow/Reap. It is the second that will be of interest in the code below. At the end I'll include a few remarks along the lines of explaining the data structure, since there is no formal documentation for this.

I kept the code more or less in the same style as the original, and in particular it remains an on-line version (that is, elements need not all go in at the outset but can be added individually). Changed a few names. Made the basic structure akin to

node (bounding box, value, zero or four subnodes)

If there are subnodes then the value field is empty. The box and value fields are represented by the usual Mathematica List expression, though it might make sense to use dedicated heads and have it more akin to a C struct style. I did do something like that in naming the various field accessing/setting functions.

One caveat is that this raw data type consumes substantially more memory overhead than e.g. a List. So my variant below will use more memory than the originally posted code. Not asymptotically more, just by a constant factor. Also it requires a constant factor in overhead more than, say, a comparable C struct in terms of accessing or setting element value. So it's not a magic bullet, just a data type with behaviour that should not give asymptotic surprises.

---

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]]

Here is an example. I will note that the scaling is reasonable. Maybe O(n log(n)) or so. Definitely better than 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}

---

General expr bag notes. These are old so I do not claim that this all still works as indicated.

These functions live in Internal` context.

Bag Creates an expr bag, optionally with preset elements.

BagPart Obtains parts of an expr bag, similar to Part for ordinary exprs. Also can be used on a lhs, e.g. to reset a value.

StuffBag Appends elements to end of a bag.

We also have a BagLength. Useful for iterating over a bag.

These functions are extrememly useful for two reasons.

First, this is a good way to make an extensible table in Mathematica.

Second, the contents of bags are evaluated but then placed in a raw expr, hence are shielded. Thus one can use these as "pointers" (in the C sense) rather than as objects, and this requires no Hold etc. Here are some examples:

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

If we instead use bags we get a self-referential structure.

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

This is difficult to emulate in any other way in Mathematica. One would need to use sparse tables (hashing) in some not-very-transparent way.

Here is a related example, not fully debugged. We basically implement a linked list whereby one can destructively modify tails, replace sublists, etc.

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>]}

So tt is a linked list, the first element is vv, the next is itself a linked list, etc. I refrained from using Lisp terminology (car/cdr and the like) because I am unable to recall whether Lisp's list operations are destructive. But you get the general idea.

Along similar lines, I have used expr bags to implement binary trees. This is useful because we can do destructive changes in constant time (assuming we already have a "handle" on the point of insertion/deletion), and moreover the "raw" nature of expr bags means we completely avoid Mathematica's infinite evaluation semantics.

Another application, perhaps.

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])

Try with

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]

---

Daniel Lichtblau Wolfram Research

Amias

Excellent information! Seems to be a very fast way to build up lists – Dyeing

+1, I just read through this for the first time. Good info, thanks. – Acerate

Answer 3

This may not be what you're trying to do, but Nearest[] can create a NearestFunction[] that is a built in quadtree structure.

Champerty