add RTree classes
[anneal.git] / src / com / infomatiq / jsi / rtree / RTree.java
1 //   RTree.java
2 //   Java Spatial Index Library
3 //   Copyright (C) 2002 Infomatiq Limited
4 //  
5 //  This library is free software; you can redistribute it and/or
6 //  modify it under the terms of the GNU Lesser General Public
7 //  License as published by the Free Software Foundation; either
8 //  version 2.1 of the License, or (at your option) any later version.
9 //  
10 //  This library is distributed in the hope that it will be useful,
11 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13 //  Lesser General Public License for more details.
14 //  
15 //  You should have received a copy of the GNU Lesser General Public
16 //  License along with this library; if not, write to the Free Software
17 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18
19 package com.infomatiq.jsi.rtree;
20
21 import gnu.trove.TIntArrayList;
22 import gnu.trove.TIntObjectHashMap;
23 import gnu.trove.TIntProcedure;
24 import gnu.trove.TIntStack;
25 import java.util.Properties;
26
27 import com.infomatiq.jsi.IntProcedure;
28 import com.infomatiq.jsi.Point;
29 import com.infomatiq.jsi.Rectangle;
30 import com.infomatiq.jsi.SpatialIndex;
31
32 /**
33  * <p>This is a lightweight RTree implementation, specifically designed 
34  * for the following features (in order of importance): 
35  * <ul>
36  * <li>Fast intersection query performance. To achieve this, the RTree 
37  * uses only main memory to store entries. Obviously this will only improve
38  * performance if there is enough physical memory to avoid paging.</li>
39  * <li>Low memory requirements.</li>
40  * <li>Fast add performance.</li>
41  * </ul></p> 
42  * 
43  * <p>The main reason for the high speed of this RTree implementation is the 
44  * avoidance of the creation of unnecessary objects, mainly achieved by using
45  * primitive collections from the trove4j library.</p>
46  * 
47  * @author aled.morris@infomatiq.co.uk
48  * @version 1.0b2
49  */
50 public class RTree implements SpatialIndex {
51   
52   private static final String version = "1.0b2";
53   
54   // parameters of the tree
55   private final static int DEFAULT_MAX_NODE_ENTRIES = 10;
56   int maxNodeEntries;
57   int minNodeEntries;
58   
59   // map of nodeId -> node object
60   // [x] TODO eliminate this map - it should not be needed. Nodes
61   // can be found by traversing the tree.
62   private TIntObjectHashMap nodeMap = new TIntObjectHashMap();
63   
64   // internal consistency checking - set to true if debugging tree corruption
65   private final static boolean INTERNAL_CONSISTENCY_CHECKING = false;
66   
67   // used to mark the status of entries during a node split
68   private final static int ENTRY_STATUS_ASSIGNED = 0;
69   private final static int ENTRY_STATUS_UNASSIGNED = 1; 
70   private byte[] entryStatus = null;
71   private byte[] initialEntryStatus = null;
72   
73   // stacks used to store nodeId and entry index of each node 
74   // from the root down to the leaf. Enables fast lookup
75   // of nodes when a split is propagated up the tree.
76   private TIntStack parents = new TIntStack();
77   private TIntStack parentsEntry = new TIntStack();
78   
79   // initialisation
80   private int treeHeight = 1; // leaves are always level 1
81   private int rootNodeId = 0;
82   private int size = 0;
83   
84   // Enables creation of new nodes
85   private int highestUsedNodeId = rootNodeId; 
86   
87   // Deleted node objects are retained in the nodeMap, 
88   // so that they can be reused. Store the IDs of nodes
89   // which can be reused.
90   private TIntStack deletedNodeIds = new TIntStack();
91   
92   // List of nearest rectangles. Use a member variable to
93   // avoid recreating the object each time nearest() is called.
94   private TIntArrayList nearestIds = new TIntArrayList();
95   
96   // Inner class used as a bridge between the trove4j TIntProcedure
97   // and the SpatialIndex IntProcedure. This is used because 
98   // the nearest rectangles must be stored as they are found, in
99   // case a closer one is found later. 
100   // 
101   // A single instance of this class is used to avoid creating a new 
102   // one every time nearest() is called.
103   private class TIntProcedureVisit implements TIntProcedure {
104         public IntProcedure m_intProcedure = null;
105         
106         public void setProcedure(IntProcedure ip) {
107           m_intProcedure = ip;
108         }
109         public boolean execute(int i) {
110           m_intProcedure.execute(i);
111           return true;
112         }
113   }; 
114   private TIntProcedureVisit visitProc = new TIntProcedureVisit();
115   
116   /**
117    * Constructor. Use init() method to initialize parameters of the RTree.
118    */
119   public RTree() {  
120     return; // NOP    
121   }
122   
123   //-------------------------------------------------------------------------
124   // public implementation of SpatialIndex interface:
125   //  init(Properties)
126   //  add(Rectangle, int)
127   //  delete(Rectangle, int)
128   //  nearest(Point, IntProcedure, float)
129   //  intersects(Rectangle, IntProcedure)
130   //  contains(Rectangle, IntProcedure)
131   //  size()
132   //-------------------------------------------------------------------------
133   /**
134    * <p>Initialize implementation dependent properties of the RTree.
135    * Currently implemented properties are:
136    * <ul>
137    * <li>MaxNodeEntries</li> This specifies the maximum number of entries
138    * in a node. The default value is 10, which is used if the property is
139    * not specified, or is less than 2.
140    * <li>MinNodeEntries</li> This specifies the minimum number of entries
141    * in a node. The default value is half of the MaxNodeEntries value (rounded
142    * down), which is used if the property is not specified or is less than 1.
143    * </ul></p>
144    * 
145    * @see com.infomatiq.jsi.SpatialIndex#init(Properties)
146    */
147   public void init(Properties props) {
148     maxNodeEntries = Integer.parseInt(props.getProperty("MaxNodeEntries", "0"));
149     minNodeEntries = Integer.parseInt(props.getProperty("MinNodeEntries", "0"));
150     
151     // Obviously a node with less than 2 entries cannot be split.
152     // The node splitting algorithm will work with only 2 entries
153     // per node, but will be inefficient.
154     if (maxNodeEntries < 2) { 
155       System.err.println("Invalid MaxNodeEntries = " + maxNodeEntries + " Resetting to default value of " + DEFAULT_MAX_NODE_ENTRIES);
156       maxNodeEntries = DEFAULT_MAX_NODE_ENTRIES;
157     }
158     
159     // The MinNodeEntries must be less than or equal to (int) (MaxNodeEntries / 2)
160     if (minNodeEntries < 1 || minNodeEntries > maxNodeEntries / 2) {
161       System.err.println("MinNodeEntries must be between 1 and MaxNodeEntries / 2");
162       minNodeEntries = maxNodeEntries / 2;
163     }
164     
165     entryStatus = new byte[maxNodeEntries];  
166     initialEntryStatus = new byte[maxNodeEntries];
167     
168     for (int i = 0; i < maxNodeEntries; i++) {
169       initialEntryStatus[i] = ENTRY_STATUS_UNASSIGNED;
170     }
171     
172     Node root = new Node(rootNodeId, 1, maxNodeEntries);
173     nodeMap.put(rootNodeId, root);
174     
175   }
176   
177   /**
178    * @see com.infomatiq.jsi.SpatialIndex#add(Rectangle, int)
179    */
180   public void add(Rectangle r, int id) {
181       /*
182     if (log.isDebugEnabled()) {
183     //////log.debug("Adding rectangle " + r + ", id " + id);
184     }
185       */
186     add(r.copy(), id, 1); 
187     
188     size++; 
189   }
190   
191   /**
192    * Adds a new entry at a specified level in the tree
193    */
194   private void add(Rectangle r, int id, int level) {
195     // I1 [Find position for new record] Invoke ChooseLeaf to select a 
196     // leaf node L in which to place r
197     Node n = chooseNode(r, level);
198     Node newLeaf = null;
199     
200     // I2 [Add record to leaf node] If L has room for another entry, 
201     // install E. Otherwise invoke SplitNode to obtain L and LL containing
202     // E and all the old entries of L
203     if (n.entryCount < maxNodeEntries) {
204       n.addEntryNoCopy(r, id);
205     } else {
206       newLeaf = splitNode(n, r, id);  
207     }
208     
209     // I3 [Propagate changes upwards] Invoke AdjustTree on L, also passing LL
210     // if a split was performed
211     Node newNode = adjustTree(n, newLeaf); 
212
213     // I4 [Grow tree taller] If node split propagation caused the root to 
214     // split, create a new root whose children are the two resulting nodes.
215     if (newNode != null) {
216       int oldRootNodeId = rootNodeId;
217       Node oldRoot = getNode(oldRootNodeId);
218       
219       rootNodeId = getNextNodeId();
220       treeHeight++;
221       Node root = new Node(rootNodeId, treeHeight, maxNodeEntries);
222       root.addEntry(newNode.mbr, newNode.nodeId);
223       root.addEntry(oldRoot.mbr, oldRoot.nodeId);
224       nodeMap.put(rootNodeId, root);
225     }
226     
227     if (INTERNAL_CONSISTENCY_CHECKING) {
228       checkConsistency(rootNodeId, treeHeight, null);
229     }
230   } 
231   
232   /**
233    * @see com.infomatiq.jsi.SpatialIndex#delete(Rectangle, int)
234    */
235   public boolean delete(Rectangle r, int id) {
236    // FindLeaf algorithm inlined here. Note the "official" algorithm 
237    // searches all overlapping entries. This seems inefficient to me, 
238    // as an entry is only worth searching if it contains (NOT overlaps)
239    // the rectangle we are searching for.
240    //
241    // Also the algorithm has been changed so that it is not recursive.
242     
243     // FL1 [Search subtrees] If root is not a leaf, check each entry 
244     // to determine if it contains r. For each entry found, invoke
245     // findLeaf on the node pointed to by the entry, until r is found or
246     // all entries have been checked.
247         parents.clear();
248         parents.push(rootNodeId);
249         
250         parentsEntry.clear();
251         parentsEntry.push(-1);
252         Node n = null;
253         int foundIndex = -1;  // index of entry to be deleted in leaf
254         
255         while (foundIndex == -1 && parents.size() > 0) {
256           n = getNode(parents.peek());
257           int startIndex = parentsEntry.peek() + 1;
258       
259       if (!n.isLeaf()) {
260           /*
261             //deleteLog.debug("searching node " + n.nodeId + ", from index " + startIndex);
262           */
263                   boolean contains = false;
264         for (int i = startIndex; i < n.entryCount; i++) {
265                     if (n.entries[i].contains(r)) {
266                       parents.push(n.ids[i]);
267                       parentsEntry.pop();
268                       parentsEntry.push(i); // this becomes the start index when the child has been searched
269                       parentsEntry.push(-1);
270                       contains = true;
271             break; // ie go to next iteration of while()
272                     }
273                   }
274         if (contains) {
275           continue;
276         }
277       } else {
278         foundIndex = n.findEntry(r, id);        
279       }
280       
281       parents.pop();
282       parentsEntry.pop();
283         } // while not found
284         
285         if (foundIndex != -1) {
286           n.deleteEntry(foundIndex, minNodeEntries);
287       condenseTree(n);
288       size--;
289         }
290     
291     return (foundIndex != -1);
292   }
293   
294   /**
295    * @see com.infomatiq.jsi.SpatialIndex#nearest(Point, IntProcedure, float)
296    */
297   public void nearest(Point p, IntProcedure v, float furthestDistance) {
298     Node rootNode = getNode(rootNodeId);
299    
300     nearest(p, rootNode, furthestDistance);
301    
302     visitProc.setProcedure(v);
303     nearestIds.forEach(visitProc);
304     nearestIds.clear();
305   }
306    
307   /**
308    * @see com.infomatiq.jsi.SpatialIndex#intersects(Rectangle, IntProcedure)
309    */
310   public void intersects(Rectangle r, IntProcedure v) {
311     Node rootNode = getNode(rootNodeId);
312     intersects(r, v, rootNode);
313   }
314
315   /**
316    * @see com.infomatiq.jsi.SpatialIndex#contains(Rectangle, IntProcedure)
317    */
318   public void contains(Rectangle r, IntProcedure v) {
319     // find all rectangles in the tree that are contained by the passed rectangle
320     // written to be non-recursive (should model other searches on this?)
321         
322     parents.clear();
323     parents.push(rootNodeId);
324     
325     parentsEntry.clear();
326     parentsEntry.push(-1);
327     
328     // TODO: possible shortcut here - could test for intersection with the 
329     // MBR of the root node. If no intersection, return immediately.
330     
331     while (parents.size() > 0) {
332       Node n = getNode(parents.peek());
333       int startIndex = parentsEntry.peek() + 1;
334       
335       if (!n.isLeaf()) {
336         // go through every entry in the index node to check
337         // if it intersects the passed rectangle. If so, it 
338         // could contain entries that are contained.
339         boolean intersects = false;
340         for (int i = startIndex; i < n.entryCount; i++) {
341           if (r.intersects(n.entries[i])) {
342             parents.push(n.ids[i]);
343             parentsEntry.pop();
344             parentsEntry.push(i); // this becomes the start index when the child has been searched
345             parentsEntry.push(-1);
346             intersects = true;
347             break; // ie go to next iteration of while()
348           }
349         }
350         if (intersects) {
351           continue;
352         }
353       } else {
354         // go through every entry in the leaf to check if 
355         // it is contained by the passed rectangle
356         for (int i = 0; i < n.entryCount; i++) {
357           if (r.contains(n.entries[i])) {
358             v.execute(n.ids[i]);
359           } 
360         }                       
361       }
362       parents.pop();
363       parentsEntry.pop();  
364     }
365   }
366
367   /**
368    * @see com.infomatiq.jsi.SpatialIndex#size()
369    */
370   public int size() {
371     return size;
372   }
373
374   /**
375    * @see com.infomatiq.jsi.SpatialIndex#getBounds()
376    */
377   public Rectangle getBounds() {
378     Rectangle bounds = null;
379     
380     Node n = getNode(getRootNodeId());
381     if (n != null && n.getMBR() != null) {
382       bounds = n.getMBR().copy();
383     }
384     return bounds;
385   }
386     
387   /**
388    * @see com.infomatiq.jsi.SpatialIndex#getVersion()
389    */
390   public String getVersion() {
391     return "RTree-" + version;
392   }
393   //-------------------------------------------------------------------------
394   // end of SpatialIndex methods
395   //-------------------------------------------------------------------------
396   
397   /**
398    * Get the next available node ID. Reuse deleted node IDs if
399    * possible
400    */
401   private int getNextNodeId() {
402     int nextNodeId = 0;
403     if (deletedNodeIds.size() > 0) {
404       nextNodeId = deletedNodeIds.pop();
405     } else {
406       nextNodeId = 1 + highestUsedNodeId++;
407     }
408     return nextNodeId;
409   }
410
411   /**
412    * Get a node object, given the ID of the node.
413    */
414   public Node getNode(int index) {
415     return (Node) nodeMap.get(index);
416   }
417
418   /**
419    * Get the highest used node ID
420    */  
421   public int getHighestUsedNodeId() {
422     return highestUsedNodeId;
423   }
424
425   /**
426    * Get the root node ID
427    */
428   public int getRootNodeId() {
429     return rootNodeId; 
430   }
431       
432   /**
433    * Split a node. Algorithm is taken pretty much verbatim from
434    * Guttman's original paper.
435    * 
436    * @return new node object.
437    */
438   private Node splitNode(Node n, Rectangle newRect, int newId) {
439     // [Pick first entry for each group] Apply algorithm pickSeeds to 
440     // choose two entries to be the first elements of the groups. Assign
441     // each to a group.
442     
443     // debug code
444     float initialArea = 0;
445     /*
446     if (log.isDebugEnabled()) {
447       Rectangle union = n.mbr.union(newRect);
448       initialArea = union.area();
449     }
450     */ 
451     System.arraycopy(initialEntryStatus, 0, entryStatus, 0, maxNodeEntries);
452     
453     Node newNode = null;
454     newNode = new Node(getNextNodeId(), n.level, maxNodeEntries);
455     nodeMap.put(newNode.nodeId, newNode);
456     
457     pickSeeds(n, newRect, newId, newNode); // this also sets the entryCount to 1
458     
459     // [Check if done] If all entries have been assigned, stop. If one
460     // group has so few entries that all the rest must be assigned to it in 
461     // order for it to have the minimum number m, assign them and stop. 
462     while (n.entryCount + newNode.entryCount < maxNodeEntries + 1) {
463       if (maxNodeEntries + 1 - newNode.entryCount == minNodeEntries) {
464         // assign all remaining entries to original node
465         for (int i = 0; i < maxNodeEntries; i++) {
466           if (entryStatus[i] == ENTRY_STATUS_UNASSIGNED) {
467             entryStatus[i] = ENTRY_STATUS_ASSIGNED;
468             n.mbr.add(n.entries[i]);
469             n.entryCount++;
470           }
471         }
472         break;
473       }   
474       if (maxNodeEntries + 1 - n.entryCount == minNodeEntries) {
475         // assign all remaining entries to new node
476         for (int i = 0; i < maxNodeEntries; i++) {
477           if (entryStatus[i] == ENTRY_STATUS_UNASSIGNED) {
478             entryStatus[i] = ENTRY_STATUS_ASSIGNED;
479             newNode.addEntryNoCopy(n.entries[i], n.ids[i]);
480             n.entries[i] = null;
481           }
482         }
483         break;
484       }
485       
486       // [Select entry to assign] Invoke algorithm pickNext to choose the
487       // next entry to assign. Add it to the group whose covering rectangle 
488       // will have to be enlarged least to accommodate it. Resolve ties
489       // by adding the entry to the group with smaller area, then to the 
490       // the one with fewer entries, then to either. Repeat from S2
491       pickNext(n, newNode);   
492     }
493       
494     n.reorganize(this);
495     
496     // check that the MBR stored for each node is correct.
497     if (INTERNAL_CONSISTENCY_CHECKING) {
498       if (!n.mbr.equals(calculateMBR(n))) {
499           throw new Error("Error: splitNode old node MBR wrong");
500       }
501       
502       if (!newNode.mbr.equals(calculateMBR(newNode))) {
503         throw new Error("Error: splitNode new node MBR wrong");
504       }
505     }
506     
507     // debug code
508     /*
509     if (log.isDebugEnabled()) {
510       float newArea = n.mbr.area() + newNode.mbr.area();
511       float percentageIncrease = (100 * (newArea - initialArea)) / initialArea;
512       //log.debug("Node " + n.nodeId + " split. New area increased by " + percentageIncrease + "%");   
513     }
514     */
515     return newNode;
516   }
517   
518   /**
519    * Pick the seeds used to split a node.
520    * Select two entries to be the first elements of the groups
521    */
522   private void pickSeeds(Node n, Rectangle newRect, int newId, Node newNode) {
523     // Find extreme rectangles along all dimension. Along each dimension,
524     // find the entry whose rectangle has the highest low side, and the one 
525     // with the lowest high side. Record the separation.
526     float maxNormalizedSeparation = 0;
527     int highestLowIndex = 0;
528     int lowestHighIndex = 0;
529     
530     // for the purposes of picking seeds, take the MBR of the node to include
531     // the new rectangle aswell.
532     n.mbr.add(newRect);
533     
534     /*
535     if (log.isDebugEnabled()) {
536         //log.debug("pickSeeds(): NodeId = " + n.nodeId + ", newRect = " + newRect);
537     }
538     */
539     
540     for (int d = 0; d < Rectangle.DIMENSIONS; d++) {
541       float tempHighestLow = newRect.min[d];
542       int tempHighestLowIndex = -1; // -1 indicates the new rectangle is the seed
543       
544       float tempLowestHigh = newRect.max[d];
545       int tempLowestHighIndex = -1;
546       
547       for (int i = 0; i < n.entryCount; i++) {
548         float tempLow = n.entries[i].min[d];
549         if (tempLow >= tempHighestLow) {
550            tempHighestLow = tempLow;
551            tempHighestLowIndex = i;
552         } else {  // ensure that the same index cannot be both lowestHigh and highestLow
553           float tempHigh = n.entries[i].max[d];
554           if (tempHigh <= tempLowestHigh) {
555             tempLowestHigh = tempHigh;
556             tempLowestHighIndex = i;
557           }
558         }
559       
560         // PS2 [Adjust for shape of the rectangle cluster] Normalize the separations
561         // by dividing by the widths of the entire set along the corresponding
562         // dimension
563         float normalizedSeparation = (tempHighestLow - tempLowestHigh) / (n.mbr.max[d] - n.mbr.min[d]);
564         
565         if (normalizedSeparation > 1 || normalizedSeparation < -1) {
566           throw new Error("Invalid normalized separation");
567         }
568       
569         /*
570         if (log.isDebugEnabled()) {
571             log.debug("Entry " + i + ", dimension " + d + ": HighestLow = " + tempHighestLow + 
572                     " (index " + tempHighestLowIndex + ")" + ", LowestHigh = " +
573                     tempLowestHigh + " (index " + tempLowestHighIndex + ", NormalizedSeparation = " + normalizedSeparation);
574         }
575         */
576         
577         // PS3 [Select the most extreme pair] Choose the pair with the greatest
578         // normalized separation along any dimension.
579         if (normalizedSeparation > maxNormalizedSeparation) {
580           maxNormalizedSeparation = normalizedSeparation;
581           highestLowIndex = tempHighestLowIndex;
582           lowestHighIndex = tempLowestHighIndex;
583         }
584       }
585     }
586       
587     // highestLowIndex is the seed for the new node.
588     if (highestLowIndex == -1) {
589       newNode.addEntry(newRect, newId);
590     } else {
591       newNode.addEntryNoCopy(n.entries[highestLowIndex], n.ids[highestLowIndex]);
592       n.entries[highestLowIndex] = null;
593       
594       // move the new rectangle into the space vacated by the seed for the new node
595       n.entries[highestLowIndex] = newRect;
596       n.ids[highestLowIndex] = newId;
597     }
598     
599     // lowestHighIndex is the seed for the original node. 
600     if (lowestHighIndex == -1) {
601       lowestHighIndex = highestLowIndex;
602     }
603     
604     entryStatus[lowestHighIndex] = ENTRY_STATUS_ASSIGNED;
605     n.entryCount = 1;
606     n.mbr.set(n.entries[lowestHighIndex].min, n.entries[lowestHighIndex].max);
607   }
608
609   /** 
610    * Pick the next entry to be assigned to a group during a node split.
611    * 
612    * [Determine cost of putting each entry in each group] For each 
613    * entry not yet in a group, calculate the area increase required
614    * in the covering rectangles of each group  
615    */
616   private int pickNext(Node n, Node newNode) {
617     float maxDifference = Float.NEGATIVE_INFINITY;
618     int next = 0;
619     int nextGroup = 0;
620     
621     maxDifference = Float.NEGATIVE_INFINITY;
622    
623     /*
624     if (log.isDebugEnabled()) {
625         //log.debug("pickNext()");
626     }
627     */
628    
629     for (int i = 0; i < maxNodeEntries; i++) {
630       if (entryStatus[i] == ENTRY_STATUS_UNASSIGNED) {
631         
632         if (n.entries[i] == null) {
633           throw new Error("Error: Node " + n.nodeId + ", entry " + i + " is null");
634         }
635         
636         float nIncrease = n.mbr.enlargement(n.entries[i]);
637         float newNodeIncrease = newNode.mbr.enlargement(n.entries[i]);
638         float difference = Math.abs(nIncrease - newNodeIncrease);
639          
640         if (difference > maxDifference) {
641           next = i;
642           
643           if (nIncrease < newNodeIncrease) {
644             nextGroup = 0; 
645           } else if (newNodeIncrease < nIncrease) {
646             nextGroup = 1;
647           } else if (n.mbr.area() < newNode.mbr.area()) {
648             nextGroup = 0;
649           } else if (newNode.mbr.area() < n.mbr.area()) {
650             nextGroup = 1;
651           } else if (newNode.entryCount < maxNodeEntries / 2) {
652             nextGroup = 0;
653           } else {
654             nextGroup = 1;
655           }
656           maxDifference = difference; 
657         }
658         /*
659         if (log.isDebugEnabled()) {
660             log.debug("Entry " + i + " group0 increase = " + nIncrease + ", group1 increase = " + newNodeIncrease +
661                     ", diff = " + difference + ", MaxDiff = " + maxDifference + " (entry " + next + ")");
662         }
663         */
664       }
665     }
666     
667     entryStatus[next] = ENTRY_STATUS_ASSIGNED;
668       
669     if (nextGroup == 0) {
670       n.mbr.add(n.entries[next]);
671       n.entryCount++;
672     } else {
673       // move to new node.
674       newNode.addEntryNoCopy(n.entries[next], n.ids[next]);
675       n.entries[next] = null;
676     }
677     
678     return next; 
679   }
680
681   /**
682    * Recursively searches the tree for the nearest entry. Other queries
683    * call execute() on an IntProcedure when a matching entry is found; 
684    * however nearest() must store the entry Ids as it searches the tree,
685    * in case a nearer entry is found.
686    * Uses the member variable nearestIds to store the nearest
687    * entry IDs.
688    * 
689    * [x] TODO rewrite this to be non-recursive?
690    */
691   private float nearest(Point p, Node n, float nearestDistance) {
692     for (int i = 0; i < n.entryCount; i++) {
693       float tempDistance = n.entries[i].distance(p);  
694       if (n.isLeaf()) { // for leaves, the distance is an actual nearest distance 
695         if (tempDistance < nearestDistance) {
696           nearestDistance = tempDistance;
697           nearestIds.clear();
698         }
699         if (tempDistance <= nearestDistance) {
700           nearestIds.add(n.ids[i]);
701         }     
702       } else { // for index nodes, only go into them if they potentially could have
703                // a rectangle nearer than actualNearest
704          if (tempDistance <= nearestDistance) {
705            // search the child node
706            nearestDistance = nearest(p, getNode(n.ids[i]), nearestDistance);
707          }
708       }
709     }
710     return nearestDistance;
711   }
712   
713   /** 
714    * Recursively searches the tree for all intersecting entries.
715    * Immediately calls execute() on the passed IntProcedure when 
716    * a matching entry is found.
717    * 
718    * [x] TODO rewrite this to be non-recursive? Make sure it
719    * doesn't slow it down.
720    */
721   private void intersects(Rectangle r, IntProcedure v, Node n) {
722     for (int i = 0; i < n.entryCount; i++) {
723       if (r.intersects(n.entries[i])) {
724         if (n.isLeaf()) {
725           v.execute(n.ids[i]);
726         } else {
727           Node childNode = getNode(n.ids[i]);
728           intersects(r, v, childNode);
729         }
730       }
731     }
732   }
733
734   /**
735    * Used by delete(). Ensures that all nodes from the passed node
736    * up to the root have the minimum number of entries.
737    * 
738    * Note that the parent and parentEntry stacks are expected to
739    * contain the nodeIds of all parents up to the root.
740    */
741     private Rectangle oldRectangle = new Rectangle(0, 0, 0, 0, 0, 0);
742   private void condenseTree(Node l) {
743     // CT1 [Initialize] Set n=l. Set the list of eliminated
744     // nodes to be empty.
745     Node n = l;
746     Node parent = null;
747     int parentEntry = 0;
748     
749     TIntStack eliminatedNodeIds = new TIntStack();
750   
751     // CT2 [Find parent entry] If N is the root, go to CT6. Otherwise 
752     // let P be the parent of N, and let En be N's entry in P  
753     while (n.level != treeHeight) {
754       parent = getNode(parents.pop());
755       parentEntry = parentsEntry.pop();
756       
757       // CT3 [Eliminiate under-full node] If N has too few entries,
758       // delete En from P and add N to the list of eliminated nodes
759       if (n.entryCount < minNodeEntries) {
760         parent.deleteEntry(parentEntry, minNodeEntries);
761         eliminatedNodeIds.push(n.nodeId);
762       } else {
763         // CT4 [Adjust covering rectangle] If N has not been eliminated,
764         // adjust EnI to tightly contain all entries in N
765         if (!n.mbr.equals(parent.entries[parentEntry])) {
766           oldRectangle.set(parent.entries[parentEntry].min, parent.entries[parentEntry].max);
767           parent.entries[parentEntry].set(n.mbr.min, n.mbr.max);
768           parent.recalculateMBR(oldRectangle);
769         }
770       }
771       // CT5 [Move up one level in tree] Set N=P and repeat from CT2
772       n = parent;
773     }
774     
775     // CT6 [Reinsert orphaned entries] Reinsert all entries of nodes in set Q.
776     // Entries from eliminated leaf nodes are reinserted in tree leaves as in 
777     // Insert(), but entries from higher level nodes must be placed higher in 
778     // the tree, so that leaves of their dependent subtrees will be on the same
779     // level as leaves of the main tree
780     while (eliminatedNodeIds.size() > 0) {
781       Node e = getNode(eliminatedNodeIds.pop());
782       for (int j = 0; j < e.entryCount; j++) {
783         add(e.entries[j], e.ids[j], e.level); 
784         e.entries[j] = null;
785       }
786       e.entryCount = 0;
787       deletedNodeIds.push(e.nodeId);
788     }
789   }
790
791   /**
792    *  Used by add(). Chooses a leaf to add the rectangle to.
793    */
794   private Node chooseNode(Rectangle r, int level) {
795     // CL1 [Initialize] Set N to be the root node
796     Node n = getNode(rootNodeId);
797     parents.clear();
798     parentsEntry.clear();
799      
800     // CL2 [Leaf check] If N is a leaf, return N
801     while (true) {
802       if (n == null) {
803         throw new Error("Could not get root node (" + rootNodeId + ")");  
804       }
805    
806       if (n.level == level) {
807         return n;
808       }
809       
810       // CL3 [Choose subtree] If N is not at the desired level, let F be the entry in N 
811       // whose rectangle FI needs least enlargement to include EI. Resolve
812       // ties by choosing the entry with the rectangle of smaller area.
813       float leastEnlargement = n.getEntry(0).enlargement(r);
814       int index = 0; // index of rectangle in subtree
815       for (int i = 1; i < n.entryCount; i++) {
816         Rectangle tempRectangle = n.getEntry(i);
817         float tempEnlargement = tempRectangle.enlargement(r);
818         if ((tempEnlargement < leastEnlargement) ||
819             ((tempEnlargement == leastEnlargement) && 
820              (tempRectangle.area() < n.getEntry(index).area()))) {
821           index = i;
822           leastEnlargement = tempEnlargement;
823         }
824       }
825       
826       parents.push(n.nodeId);
827       parentsEntry.push(index);
828     
829       // CL4 [Descend until a leaf is reached] Set N to be the child node 
830       // pointed to by Fp and repeat from CL2
831       n = getNode(n.ids[index]);
832     }
833   }
834   
835   /**
836    * Ascend from a leaf node L to the root, adjusting covering rectangles and
837    * propagating node splits as necessary.
838    */
839   private Node adjustTree(Node n, Node nn) {
840     // AT1 [Initialize] Set N=L. If L was split previously, set NN to be 
841     // the resulting second node.
842     
843     // AT2 [Check if done] If N is the root, stop
844     while (n.level != treeHeight) {
845     
846       // AT3 [Adjust covering rectangle in parent entry] Let P be the parent 
847       // node of N, and let En be N's entry in P. Adjust EnI so that it tightly
848       // encloses all entry rectangles in N.
849       Node parent = getNode(parents.pop());
850       int entry = parentsEntry.pop(); 
851       
852       if (parent.ids[entry] != n.nodeId) {
853         throw new Error("Error: entry " + entry + " in node " + 
854              parent.nodeId + " should point to node " + 
855              n.nodeId + "; actually points to node " + parent.ids[entry]);
856       }
857       
858       if (!parent.entries[entry].equals(n.mbr)) {
859         parent.entries[entry].set(n.mbr.min, n.mbr.max);
860         parent.mbr.set(parent.entries[0].min, parent.entries[0].max);
861         for (int i = 1; i < parent.entryCount; i++) {
862           parent.mbr.add(parent.entries[i]);
863         }
864       }
865       
866       // AT4 [Propagate node split upward] If N has a partner NN resulting from 
867       // an earlier split, create a new entry Enn with Ennp pointing to NN and 
868       // Enni enclosing all rectangles in NN. Add Enn to P if there is room. 
869       // Otherwise, invoke splitNode to produce P and PP containing Enn and
870       // all P's old entries.
871       Node newNode = null;
872       if (nn != null) {
873         if (parent.entryCount < maxNodeEntries) {
874           parent.addEntry(nn.mbr, nn.nodeId);
875         } else {
876           newNode = splitNode(parent, nn.mbr.copy(), nn.nodeId);
877         }
878       }
879       
880       // AT5 [Move up to next level] Set N = P and set NN = PP if a split 
881       // occurred. Repeat from AT2
882       n = parent;
883       nn = newNode;
884       
885       parent = null;
886       newNode = null;
887     }
888     
889     return nn;
890   }
891   
892   /**
893    * Check the consistency of the tree.
894    */
895   private void checkConsistency(int nodeId, int expectedLevel, Rectangle expectedMBR) {
896     // go through the tree, and check that the internal data structures of 
897     // the tree are not corrupted.    
898     Node n = getNode(nodeId);
899     
900     if (n == null) {
901       throw new Error("Error: Could not read node " + nodeId);
902     }
903     
904     if (n.level != expectedLevel) {
905       throw new Error("Error: Node " + nodeId + ", expected level " + expectedLevel + ", actual level " + n.level);
906     }
907     
908     Rectangle calculatedMBR = calculateMBR(n);
909     
910     if (!n.mbr.equals(calculatedMBR)) {
911       throw new Error("Error: Node " + nodeId + ", calculated MBR does not equal stored MBR");
912     }
913     
914     if (expectedMBR != null && !n.mbr.equals(expectedMBR)) {
915       throw new Error("Error: Node " + nodeId + ", expected MBR (from parent) does not equal stored MBR");
916     }
917     
918     // Check for corruption where a parent entry is the same object as the child MBR
919     if (expectedMBR != null && n.mbr.sameObject(expectedMBR)) {
920       throw new Error("Error: Node " + nodeId + " MBR using same rectangle object as parent's entry");
921     }
922     
923     for (int i = 0; i < n.entryCount; i++) {
924       if (n.entries[i] == null) {
925         throw new Error("Error: Node " + nodeId + ", Entry " + i + " is null");
926       }     
927       
928       if (n.level > 1) { // if not a leaf
929         checkConsistency(n.ids[i], n.level - 1, n.entries[i]); 
930       }   
931     } 
932   }
933   
934   /**
935    * Given a node object, calculate the node MBR from it's entries.
936    * Used in consistency checking
937    */
938   private Rectangle calculateMBR(Node n) {
939     Rectangle mbr = new Rectangle(n.entries[0].min, n.entries[0].max);
940     
941     for (int i = 1; i < n.entryCount; i++) {
942       mbr.add(n.entries[i]);
943     }
944     return mbr; 
945   }
946 }