add RTree classes
authoradam <adam@megacz.com>
Fri, 14 Dec 2007 01:53:17 +0000 (17:53 -0800)
committeradam <adam@megacz.com>
Fri, 14 Dec 2007 01:53:17 +0000 (17:53 -0800)
darcs-hash:20071214015317-5007d-31f9de5bc2afdc9cade724105114fec68e940c5b.gz

14 files changed:
Makefile
src/com/infomatiq/jsi/IntProcedure.java [new file with mode: 0644]
src/com/infomatiq/jsi/Point.java [new file with mode: 0644]
src/com/infomatiq/jsi/Rectangle.java [new file with mode: 0644]
src/com/infomatiq/jsi/SpatialIndex.java [new file with mode: 0644]
src/com/infomatiq/jsi/rtree/Node.java [new file with mode: 0644]
src/com/infomatiq/jsi/rtree/RTree.java [new file with mode: 0644]
src/edu/berkeley/qfat/geom/HasBoundingBox.java [new file with mode: 0644]
src/edu/berkeley/qfat/geom/HasPoint.java
src/edu/berkeley/qfat/geom/IntervalTree.java [new file with mode: 0644]
src/edu/berkeley/qfat/geom/Point.java
src/edu/berkeley/qfat/geom/PointSet.java
src/edu/berkeley/qfat/geom/RTree.java [new file with mode: 0644]
src/edu/berkeley/qfat/geom/Triangle.java

index dec73e8..d42a2b7 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,5 @@
+jars = lib/kd.jar:lib/trove-0.1.8.jar:lib/sil-0.43b-am1.jar
 all:
        mkdir -p build
 all:
        mkdir -p build
-       javac -cp lib/kd.jar -d build `find src -name \*.java`
-       java -server -cp lib/kd.jar:build edu.berkeley.qfat.Main
+       javac -cp $(jars) -d build `find src -name \*.java`
+       java -server -cp $(jars):build edu.berkeley.qfat.Main
diff --git a/src/com/infomatiq/jsi/IntProcedure.java b/src/com/infomatiq/jsi/IntProcedure.java
new file mode 100644 (file)
index 0000000..6081136
--- /dev/null
@@ -0,0 +1,37 @@
+//   IntProcedure.java
+//   Java Spatial Index Library
+//   Copyright (C) 2002 Infomatiq Limited
+//  
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License, or (at your option) any later version.
+//  
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//  
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+package com.infomatiq.jsi;
+
+/**
+ * Interface that defines a procedure to be executed, that takes
+ * an int parameter
+ * 
+ * @author aled.morris@infomatiq.co.uk
+ * @version 1.0b2
+ */
+public interface IntProcedure {
+  /**
+   * @param id integer value
+   * 
+   * @return flag to indicate whether to continue executing the procedure. Return true to
+   *         continue executing, or false to prevent any more calls to this 
+   *         method.
+   */
+  public boolean execute(int id);
+}
diff --git a/src/com/infomatiq/jsi/Point.java b/src/com/infomatiq/jsi/Point.java
new file mode 100644 (file)
index 0000000..c608971
--- /dev/null
@@ -0,0 +1,51 @@
+//   Point.java
+//   Java Spatial Index Library
+//   Copyright (C) 2002 Infomatiq Limited
+//  
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License, or (at your option) any later version.
+//  
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//  
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+
+package com.infomatiq.jsi;
+
+/**
+ * Currently hardcoded to 2 dimensions, but could be extended.
+ * 
+ * @author  aled.morris@infomatiq.co.uk
+ * @version 1.0b2
+ */
+public class Point {
+  /**
+   * Number of dimensions in a point. In theory this
+   * could be exended to three or more dimensions.
+   */
+  private final static int DIMENSIONS = 3;
+  
+  /**
+   * The (x, y) coordinates of the point.
+   */
+  public float[] coordinates;
+  
+  /**
+   * Constructor.
+   * 
+   * @param x The x coordinate of the point
+   * @param y The y coordinate of the point
+   */
+  public Point(float x, float y, float z) {
+    coordinates = new float[DIMENSIONS];
+    coordinates[0] = x; 
+    coordinates[1] = y;
+    coordinates[2] = z;
+  }
+}
diff --git a/src/com/infomatiq/jsi/Rectangle.java b/src/com/infomatiq/jsi/Rectangle.java
new file mode 100644 (file)
index 0000000..737e6cf
--- /dev/null
@@ -0,0 +1,366 @@
+//   Rectangle.java
+//   Java Spatial Index Library
+//   Copyright (C) 2002 Infomatiq Limited
+//  
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License, or (at your option) any later version.
+//  
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//  
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+
+package com.infomatiq.jsi;
+
+import java.util.Arrays;
+
+/**
+ * Currently hardcoded to 2 dimensions, but could be extended.
+ * 
+ * @author  aled.morris@infomatiq.co.uk
+ * @version 1.0b2
+ */
+public class Rectangle {
+  /**
+   * Number of dimensions in a rectangle. In theory this
+   * could be exended to three or more dimensions.
+   */
+  public final static int DIMENSIONS = 3;
+  
+  /**
+   * array containing the minimum value for each dimension; ie { min(x), min(y) }
+   */
+  public float[] max;
+  
+  /**
+   * array containing the maximum value for each dimension; ie { max(x), max(y) }
+   */
+  public float[] min;
+
+  /**
+   * Constructor.
+   * 
+   * @param x1 coordinate of any corner of the rectangle
+   * @param y1 (see x1)
+   * @param x2 coordinate of the opposite corner
+   * @param y2 (see x2)
+   */
+  public Rectangle(float x1, float y1, float z1, float x2, float y2, float z2) {
+    min = new float[DIMENSIONS];
+    max = new float[DIMENSIONS]; 
+    set(x1, y1, z1, x2, y2, z2);
+  }
+
+  /**
+   * Constructor.
+   * 
+   * @param min array containing the minimum value for each dimension; ie { min(x), min(y) }
+   * @param max array containing the maximum value for each dimension; ie { max(x), max(y) }
+   */  
+  public Rectangle(float[] min, float[] max) {
+    if (min.length != DIMENSIONS || max.length != DIMENSIONS) {
+      throw new RuntimeException("Error in Rectangle constructor: " +
+                "min and max arrays must be of length " + DIMENSIONS); 
+    }
+   
+    this.min = new float[DIMENSIONS];
+    this.max = new float[DIMENSIONS];
+    
+    set(min, max);
+  }
+  
+ /**
+   * Sets the size of the rectangle.
+   * 
+   * @param x1 coordinate of any corner of the rectangle
+   * @param y1 (see x1)
+   * @param x2 coordinate of the opposite corner
+   * @param y2 (see x2)
+   */
+  public void set(float x1, float y1, float z1, float x2, float y2, float z2) {
+    min[0] = Math.min(x1, x2);
+    min[1] = Math.min(y1, y2);
+    min[2] = Math.min(z1, z2);
+    max[0] = Math.max(x1, x2);        
+    max[1] = Math.max(y1, y2); 
+    max[2] = Math.max(z1, z2);
+  }
+  
+  /**
+   * Sets the size of the rectangle.
+   * 
+   * @param min array containing the minimum value for each dimension; ie { min(x), min(y) }
+   * @param max array containing the maximum value for each dimension; ie { max(x), max(y) }
+   */
+  public void set(float[] min, float[] max) {
+    System.arraycopy(min, 0, this.min, 0, DIMENSIONS);
+    System.arraycopy(max, 0, this.max, 0, DIMENSIONS);
+  }
+  
+  /**
+   * Make a copy of this rectangle
+   * 
+   * @return copy of this rectangle
+   */
+  public Rectangle copy() {
+    return new Rectangle(min, max); 
+  }
+  
+  /**
+   * Determine whether an edge of this rectangle overlies the equivalent 
+   * edge of the passed rectangle
+   */
+  public boolean edgeOverlaps(Rectangle r) {
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (min[i] == r.min[i] || max[i] == r.max[i]) {
+        return true; 
+      } 
+    }  
+    return false;
+  }
+  
+  /**
+   * Determine whether this rectangle intersects the passed rectangle
+   * 
+   * @param r The rectangle that might intersect this rectangle
+   * 
+   * @return true if the rectangles intersect, false if they do not intersect
+   */
+  public boolean intersects(Rectangle r) {
+    // Every dimension must intersect. If any dimension
+    // does not intersect, return false immediately.
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (max[i] < r.min[i] || min[i] > r.max[i]) {
+        return false;
+      }
+    }
+    return true;
+  }
+  /**
+   * Determine whether this rectangle contains the passed rectangle
+   * 
+   * @param r The rectangle that might be contained by this rectangle
+   * 
+   * @return true if this rectangle contains the passed rectangle, false if
+   *         it does not
+   */
+  public boolean contains(Rectangle r) {
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (max[i] < r.max[i] || min[i] > r.min[i]) {
+        return false;
+      }
+    }
+    return true;     
+  }
+  /**
+   * Determine whether this rectangle is contained by the passed rectangle
+   * 
+   * @param r The rectangle that might contain this rectangle
+   * 
+   * @return true if the passed rectangle contains this rectangle, false if
+   *         it does not
+   */
+  public boolean containedBy(Rectangle r) {
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (max[i] > r.max[i] || min[i] < r.min[i]) {
+        return false;
+      }
+    }
+    return true;  
+  }
+  
+  /**
+   * Return the distance between this rectangle and the passed point.
+   * If the rectangle contains the point, the distance is zero.
+   * 
+   * @param p Point to find the distance to
+   * 
+   * @return distance beween this rectangle and the passed point.
+   */
+  public float distance(Point p) {
+    float distanceSquared = 0;
+    for (int i = 0; i < DIMENSIONS; i++) {
+      float greatestMin = Math.max(min[i], p.coordinates[i]);
+      float leastMax    = Math.min(max[i], p.coordinates[i]);
+      if (greatestMin > leastMax) {
+        distanceSquared += ((greatestMin - leastMax) * (greatestMin - leastMax)); 
+      }
+    }
+    return (float) Math.sqrt(distanceSquared);
+  }
+  
+  /**
+   * Return the distance between this rectangle and the passed rectangle.
+   * If the rectangles overlap, the distance is zero.
+   * 
+   * @param r Rectangle to find the distance to
+   * 
+   * @return distance between this rectangle and the passed rectangle
+   */
+
+  public float distance(Rectangle r) {
+    float distanceSquared = 0;
+    for (int i = 0; i < DIMENSIONS; i++) {
+      float greatestMin = Math.max(min[i], r.min[i]);
+      float leastMax    = Math.min(max[i], r.max[i]);
+      if (greatestMin > leastMax) {
+        distanceSquared += ((greatestMin - leastMax) * (greatestMin - leastMax)); 
+      }
+    }
+    return (float) Math.sqrt(distanceSquared);
+  }
+   
+  /**
+   * Return the squared distance from this rectangle to the passed point
+   */
+  private float distanceSquared(int dimension, float point) {
+    float distanceSquared = 0;
+    float tempDistance = point - max[dimension];
+    for (int i = 0; i < 2; i++) {
+      if (tempDistance > 0) {
+        distanceSquared = (tempDistance * tempDistance);
+        break;
+      } 
+      tempDistance = min[dimension] - point;
+    }
+    return distanceSquared;
+  }
+  
+  /**
+   * Return the furthst possible distance between this rectangle and
+   * the passed rectangle. 
+   * 
+   * Find the distance between this rectangle and each corner of the
+   * passed rectangle, and use the maximum.
+   *
+   */
+  public float furthestDistance(Rectangle r) {
+     float distanceSquared = 0;
+     
+     for (int i = 0; i < DIMENSIONS; i++) {
+       distanceSquared += Math.max(distanceSquared(i, r.min[i]), distanceSquared(i, r.max[i]));
+     }
+     
+     return (float) Math.sqrt(distanceSquared);
+  }
+  
+  /**
+   * Calculate the area by which this rectangle would be enlarged if
+   * added to the passed rectangle. Neither rectangle is altered.
+   * 
+   * @param r Rectangle to union with this rectangle, in order to 
+   *          compute the difference in area of the union and the
+   *          original rectangle
+   */
+  public float enlargement(Rectangle r) {
+    float enlargedArea = (Math.max(max[0], r.max[0]) - Math.min(min[0], r.min[0])) *
+                         (Math.max(max[1], r.max[1]) - Math.min(min[1], r.min[1]));
+                         
+    return enlargedArea - area();
+  }
+  
+  /**
+   * Compute the area of this rectangle.
+   * 
+   * @return The area of this rectangle
+   */
+  public float area() {
+    return (max[0] - min[0]) * (max[1] - min[1]);
+  }
+  
+  /**
+   * Computes the union of this rectangle and the passed rectangle, storing
+   * the result in this rectangle.
+   * 
+   * @param r Rectangle to add to this rectangle
+   */
+  public void add(Rectangle r) {
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (r.min[i] < min[i]) {
+        min[i] = r.min[i];
+      }
+      if (r.max[i] > max[i]) {
+        max[i] = r.max[i];
+      }
+    }
+  }
+  
+  /**
+   * Find the the union of this rectangle and the passed rectangle.
+   * Neither rectangle is altered
+   * 
+   * @param r The rectangle to union with this rectangle
+   */
+  public Rectangle union(Rectangle r) {
+    Rectangle union = this.copy();
+    union.add(r);
+    return union; 
+  }
+  
+  /**
+   * Determine whether this rectangle is equal to a given object.
+   * Equality is determined by the bounds of the rectangle.
+   * 
+   * @param The object to compare with this rectangle
+   */
+  public boolean equals(Object o) {
+    boolean equals = false;
+    if (o instanceof Rectangle) {
+      Rectangle r = (Rectangle) o;
+      if (Arrays.equals(r.min, min) && Arrays.equals(r.max, max)) {
+        equals = true;
+      }
+    } 
+    return equals;       
+  }
+
+  /** 
+   * Determine whether this rectangle is the same as another object
+   * 
+   * Note that two rectangles can be equal but not the same object, 
+   * if they both have the same bounds.
+   * 
+   * @param o The object to compare with this rectangle.
+   */  
+  public boolean sameObject(Object o) {
+    return super.equals(o); 
+  }
+  
+  /**
+   * Return a string representation of this rectangle, in the form: 
+   * (1.2, 3.4), (5.6, 7.8)
+   * 
+   * @return String String representation of this rectangle.
+   */
+  public String toString() {
+    StringBuffer sb = new StringBuffer();
+    
+    // min coordinates
+    sb.append('(');
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (i > 0) {
+        sb.append(", ");
+      }
+      sb.append(min[i]);
+    } 
+    sb.append("), (");
+    
+    // max coordinates
+    for (int i = 0; i < DIMENSIONS; i++) {
+      if (i > 0) {
+        sb.append(", ");
+      }
+      sb.append(max[i]);
+    } 
+    sb.append(')');
+    return sb.toString();
+  }
+}
diff --git a/src/com/infomatiq/jsi/SpatialIndex.java b/src/com/infomatiq/jsi/SpatialIndex.java
new file mode 100644 (file)
index 0000000..a385cf5
--- /dev/null
@@ -0,0 +1,126 @@
+//   SpatialIndex.java
+//   Java Spatial Index Library
+//   Copyright (C) 2002 Infomatiq Limited
+//  
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License, or (at your option) any later version.
+//  
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//  
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+package com.infomatiq.jsi;
+
+import java.util.Properties;
+
+/**
+ * Defines methods that must be implemented by all 
+ * spatial indexes. This includes the RTree and its variants.
+ * 
+ * @author  aled.morris@infomatiq.co.uk
+ * @version 1.0b2
+ */
+public interface SpatialIndex {
+  
+  /**
+   * Initializes any implementation dependent properties
+   * of the spatial index. For example, RTree implementations
+   * will have a NodeSize property.
+   * 
+   * @param props The set of properties used to initialize the spatial index.
+   */
+  public void init(Properties props);
+  
+  /**
+   * Adds a new rectangle to the spatial index
+   * 
+   * @param r  The rectangle to add to the spatial index.
+   * @param id The ID of the rectangle to add to the spatial index.
+   *           The result of adding more than one rectangle with
+   *           the same ID is undefined.
+   */ 
+  public void add(Rectangle r, int id);
+  
+  /**
+   * Deletes a rectangle from the spatial index
+   * 
+   * @param r  The rectangle to delete from the spatial index
+   * @param id The ID of the rectangle to delete from the spatial
+   *           index
+   * 
+   * @return true  if the rectangle was deleted
+   *         false if the rectangle was not found, or the 
+   *               rectangle was found but with a different ID
+   */
+  public boolean delete(Rectangle r, int id);
+   
+  /**
+   * Finds all rectangles that are nearest to the passed rectangle, and calls 
+   * execute() on the passed IntProcedure for each one. 
+   * 
+   * @param p The point for which this method finds the
+   * nearest neighbours.
+   * 
+   * @param ip The IntProcedure whose execute() method is is called
+   * for each nearest neighbour.
+   * 
+   * @param distance The furthest distance away from the rectangle
+   * to search. Rectangles further than this will not be found. 
+   * 
+   * This should be as small as possible to minimise
+   * the search time.
+   *                         
+   * Use Float.POSITIVE_INFINITY to guarantee that the nearest rectangle is found,
+   * no matter how far away, although this will slow down the algorithm.
+   */
+  public void nearest(Point p, IntProcedure v, float distance);
+  
+  /**
+   * Finds all rectangles that intersect the passed rectangle.
+   * 
+   * @param  r The rectangle for which this method finds
+   *           intersecting rectangles.
+   * 
+   * @param ip The IntProcedure whose execute() method is is called
+   *           for each intersecting rectangle.
+   */
+  public void intersects(Rectangle r, IntProcedure ip);  
+
+  /**
+   * Finds all rectangles contained by the passed rectangle.
+   * 
+   * @param r The rectangle for which this method finds
+   *           contained rectangles.
+   * 
+   * @param v The visitor whose visit() method is is called
+   *           for each contained rectangle.
+   */
+  public void contains(Rectangle r, IntProcedure ip); 
+  
+  /**
+   * Returns the number of entries in the spatial index
+   */
+  public int size();
+  
+  
+  /**
+   * Returns the bounds of all the entries in the spatial index,
+   * or null if there are no entries.
+   */
+  public Rectangle getBounds();
+  
+  /**
+   * Returns a string identifying the type of
+   * spatial index, and the version number, 
+   * eg "SimpleIndex-0.1"
+   */
+  public String getVersion();
+  
+}
diff --git a/src/com/infomatiq/jsi/rtree/Node.java b/src/com/infomatiq/jsi/rtree/Node.java
new file mode 100644 (file)
index 0000000..ac00dc0
--- /dev/null
@@ -0,0 +1,154 @@
+//   Node.java
+//   Java Spatial Index Library
+//   Copyright (C) 2002 Infomatiq Limited
+//  
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License, or (at your option) any later version.
+//  
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//  
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+package com.infomatiq.jsi.rtree;
+
+import com.infomatiq.jsi.Rectangle;
+
+/**
+ * <p>Used by RTree. There are no public methods in this class.</p>
+ * 
+ * @author aled.morris@infomatiq.co.uk
+ * @version 1.0b2
+ */
+public class Node {
+  int nodeId = 0;
+  Rectangle mbr = null;
+  Rectangle[] entries = null;
+  int[] ids = null;
+  int level;
+  int entryCount;
+
+  Node(int nodeId, int level, int maxNodeEntries) {
+    this.nodeId = nodeId;
+    this.level = level;
+    entries = new Rectangle[maxNodeEntries];
+    ids = new int[maxNodeEntries];
+  }
+   
+  void addEntry(Rectangle r, int id) {
+    ids[entryCount] = id;
+    entries[entryCount] = r.copy();
+    entryCount++;
+    if (mbr == null) {
+      mbr = r.copy();
+    } else {
+      mbr.add(r);
+    }
+  }
+  
+  void addEntryNoCopy(Rectangle r, int id) {
+    ids[entryCount] = id;
+    entries[entryCount] = r;
+    entryCount++;
+    if (mbr == null) {
+      mbr = r.copy();
+    } else {
+      mbr.add(r);
+    }
+  }
+  
+  // Return the index of the found entry, or -1 if not found
+  int findEntry(Rectangle r, int id) {
+    for (int i = 0; i < entryCount; i++) {
+       if (id == ids[i] && r.equals(entries[i])) {
+         return i;     
+       }
+    }
+    return -1;
+  }
+  
+  // delete entry. This is done by setting it to null and copying the last entry into its space.
+  void deleteEntry(int i, int minNodeEntries) {
+         int lastIndex = entryCount - 1;
+         Rectangle deletedRectangle = entries[i];
+    entries[i] = null;
+         if (i != lastIndex) {
+               entries[i] = entries[lastIndex];
+               ids[i] = ids[lastIndex];
+      entries[lastIndex] = null;
+         }
+    entryCount--;
+    
+    // if there are at least minNodeEntries, adjust the MBR.
+    // otherwise, don't bother, as the node will be 
+    // eliminated anyway.
+    if (entryCount >= minNodeEntries) {
+      recalculateMBR(deletedRectangle);
+    }
+  } 
+  
+  // oldRectangle is a rectangle that has just been deleted or made smaller.
+  // Thus, the MBR is only recalculated if the OldRectangle influenced the old MBR
+  void recalculateMBR(Rectangle deletedRectangle) {
+    if (mbr.edgeOverlaps(deletedRectangle)) { 
+      mbr.set(entries[0].min, entries[0].max);
+      
+      for (int i = 1; i < entryCount; i++) {
+        mbr.add(entries[i]);
+      }
+    }
+  }
+   
+  public int getEntryCount() {
+    return entryCount;
+  }
+  
+  public Rectangle getEntry(int index) {
+    if (index < entryCount) {
+      return entries[index];
+    }
+    return null;
+  }
+  
+  public int getId(int index) {
+    if (index < entryCount) {
+      return ids[index];
+    }
+    return -1;
+  }
+  
+  /**
+   * eliminate null entries, move all entries to the start of the source node
+   */
+  void reorganize(RTree rtree) {
+    int countdownIndex = rtree.maxNodeEntries - 1; 
+    for (int index = 0; index < entryCount; index++) {
+      if (entries[index] == null) {
+         while (entries[countdownIndex] == null && countdownIndex > index) {
+           countdownIndex--;
+         }
+         entries[index] = entries[countdownIndex];
+         ids[index] = ids[countdownIndex];    
+         entries[countdownIndex] = null;
+      }
+    }
+  }
+  
+  boolean isLeaf() {
+    return (level == 1);
+  }
+  
+  public int getLevel() {
+    return level; 
+  }
+  
+  public Rectangle getMBR() {
+    return mbr; 
+  }
+}
diff --git a/src/com/infomatiq/jsi/rtree/RTree.java b/src/com/infomatiq/jsi/rtree/RTree.java
new file mode 100644 (file)
index 0000000..ba9b130
--- /dev/null
@@ -0,0 +1,946 @@
+//   RTree.java
+//   Java Spatial Index Library
+//   Copyright (C) 2002 Infomatiq Limited
+//  
+//  This library is free software; you can redistribute it and/or
+//  modify it under the terms of the GNU Lesser General Public
+//  License as published by the Free Software Foundation; either
+//  version 2.1 of the License, or (at your option) any later version.
+//  
+//  This library is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+//  Lesser General Public License for more details.
+//  
+//  You should have received a copy of the GNU Lesser General Public
+//  License along with this library; if not, write to the Free Software
+//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+package com.infomatiq.jsi.rtree;
+
+import gnu.trove.TIntArrayList;
+import gnu.trove.TIntObjectHashMap;
+import gnu.trove.TIntProcedure;
+import gnu.trove.TIntStack;
+import java.util.Properties;
+
+import com.infomatiq.jsi.IntProcedure;
+import com.infomatiq.jsi.Point;
+import com.infomatiq.jsi.Rectangle;
+import com.infomatiq.jsi.SpatialIndex;
+
+/**
+ * <p>This is a lightweight RTree implementation, specifically designed 
+ * for the following features (in order of importance): 
+ * <ul>
+ * <li>Fast intersection query performance. To achieve this, the RTree 
+ * uses only main memory to store entries. Obviously this will only improve
+ * performance if there is enough physical memory to avoid paging.</li>
+ * <li>Low memory requirements.</li>
+ * <li>Fast add performance.</li>
+ * </ul></p> 
+ * 
+ * <p>The main reason for the high speed of this RTree implementation is the 
+ * avoidance of the creation of unnecessary objects, mainly achieved by using
+ * primitive collections from the trove4j library.</p>
+ * 
+ * @author aled.morris@infomatiq.co.uk
+ * @version 1.0b2
+ */
+public class RTree implements SpatialIndex {
+  
+  private static final String version = "1.0b2";
+  
+  // parameters of the tree
+  private final static int DEFAULT_MAX_NODE_ENTRIES = 10;
+  int maxNodeEntries;
+  int minNodeEntries;
+  
+  // map of nodeId -> node object
+  // [x] TODO eliminate this map - it should not be needed. Nodes
+  // can be found by traversing the tree.
+  private TIntObjectHashMap nodeMap = new TIntObjectHashMap();
+  
+  // internal consistency checking - set to true if debugging tree corruption
+  private final static boolean INTERNAL_CONSISTENCY_CHECKING = false;
+  
+  // used to mark the status of entries during a node split
+  private final static int ENTRY_STATUS_ASSIGNED = 0;
+  private final static int ENTRY_STATUS_UNASSIGNED = 1; 
+  private byte[] entryStatus = null;
+  private byte[] initialEntryStatus = null;
+  
+  // stacks used to store nodeId and entry index of each node 
+  // from the root down to the leaf. Enables fast lookup
+  // of nodes when a split is propagated up the tree.
+  private TIntStack parents = new TIntStack();
+  private TIntStack parentsEntry = new TIntStack();
+  
+  // initialisation
+  private int treeHeight = 1; // leaves are always level 1
+  private int rootNodeId = 0;
+  private int size = 0;
+  
+  // Enables creation of new nodes
+  private int highestUsedNodeId = rootNodeId; 
+  
+  // Deleted node objects are retained in the nodeMap, 
+  // so that they can be reused. Store the IDs of nodes
+  // which can be reused.
+  private TIntStack deletedNodeIds = new TIntStack();
+  
+  // List of nearest rectangles. Use a member variable to
+  // avoid recreating the object each time nearest() is called.
+  private TIntArrayList nearestIds = new TIntArrayList();
+  
+  // Inner class used as a bridge between the trove4j TIntProcedure
+  // and the SpatialIndex IntProcedure. This is used because 
+  // the nearest rectangles must be stored as they are found, in
+  // case a closer one is found later. 
+  // 
+  // A single instance of this class is used to avoid creating a new 
+  // one every time nearest() is called.
+  private class TIntProcedureVisit implements TIntProcedure {
+       public IntProcedure m_intProcedure = null;
+       
+       public void setProcedure(IntProcedure ip) {
+         m_intProcedure = ip;
+       }
+       public boolean execute(int i) {
+         m_intProcedure.execute(i);
+         return true;
+       }
+  }; 
+  private TIntProcedureVisit visitProc = new TIntProcedureVisit();
+  
+  /**
+   * Constructor. Use init() method to initialize parameters of the RTree.
+   */
+  public RTree() {  
+    return; // NOP    
+  }
+  
+  //-------------------------------------------------------------------------
+  // public implementation of SpatialIndex interface:
+  //  init(Properties)
+  //  add(Rectangle, int)
+  //  delete(Rectangle, int)
+  //  nearest(Point, IntProcedure, float)
+  //  intersects(Rectangle, IntProcedure)
+  //  contains(Rectangle, IntProcedure)
+  //  size()
+  //-------------------------------------------------------------------------
+  /**
+   * <p>Initialize implementation dependent properties of the RTree.
+   * Currently implemented properties are:
+   * <ul>
+   * <li>MaxNodeEntries</li> This specifies the maximum number of entries
+   * in a node. The default value is 10, which is used if the property is
+   * not specified, or is less than 2.
+   * <li>MinNodeEntries</li> This specifies the minimum number of entries
+   * in a node. The default value is half of the MaxNodeEntries value (rounded
+   * down), which is used if the property is not specified or is less than 1.
+   * </ul></p>
+   * 
+   * @see com.infomatiq.jsi.SpatialIndex#init(Properties)
+   */
+  public void init(Properties props) {
+    maxNodeEntries = Integer.parseInt(props.getProperty("MaxNodeEntries", "0"));
+    minNodeEntries = Integer.parseInt(props.getProperty("MinNodeEntries", "0"));
+    
+    // Obviously a node with less than 2 entries cannot be split.
+    // The node splitting algorithm will work with only 2 entries
+    // per node, but will be inefficient.
+    if (maxNodeEntries < 2) { 
+      System.err.println("Invalid MaxNodeEntries = " + maxNodeEntries + " Resetting to default value of " + DEFAULT_MAX_NODE_ENTRIES);
+      maxNodeEntries = DEFAULT_MAX_NODE_ENTRIES;
+    }
+    
+    // The MinNodeEntries must be less than or equal to (int) (MaxNodeEntries / 2)
+    if (minNodeEntries < 1 || minNodeEntries > maxNodeEntries / 2) {
+      System.err.println("MinNodeEntries must be between 1 and MaxNodeEntries / 2");
+      minNodeEntries = maxNodeEntries / 2;
+    }
+    
+    entryStatus = new byte[maxNodeEntries];  
+    initialEntryStatus = new byte[maxNodeEntries];
+    
+    for (int i = 0; i < maxNodeEntries; i++) {
+      initialEntryStatus[i] = ENTRY_STATUS_UNASSIGNED;
+    }
+    
+    Node root = new Node(rootNodeId, 1, maxNodeEntries);
+    nodeMap.put(rootNodeId, root);
+    
+  }
+  
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#add(Rectangle, int)
+   */
+  public void add(Rectangle r, int id) {
+      /*
+    if (log.isDebugEnabled()) {
+    //////log.debug("Adding rectangle " + r + ", id " + id);
+    }
+      */
+    add(r.copy(), id, 1); 
+    
+    size++; 
+  }
+  
+  /**
+   * Adds a new entry at a specified level in the tree
+   */
+  private void add(Rectangle r, int id, int level) {
+    // I1 [Find position for new record] Invoke ChooseLeaf to select a 
+    // leaf node L in which to place r
+    Node n = chooseNode(r, level);
+    Node newLeaf = null;
+    
+    // I2 [Add record to leaf node] If L has room for another entry, 
+    // install E. Otherwise invoke SplitNode to obtain L and LL containing
+    // E and all the old entries of L
+    if (n.entryCount < maxNodeEntries) {
+      n.addEntryNoCopy(r, id);
+    } else {
+      newLeaf = splitNode(n, r, id);  
+    }
+    
+    // I3 [Propagate changes upwards] Invoke AdjustTree on L, also passing LL
+    // if a split was performed
+    Node newNode = adjustTree(n, newLeaf); 
+
+    // I4 [Grow tree taller] If node split propagation caused the root to 
+    // split, create a new root whose children are the two resulting nodes.
+    if (newNode != null) {
+      int oldRootNodeId = rootNodeId;
+      Node oldRoot = getNode(oldRootNodeId);
+      
+      rootNodeId = getNextNodeId();
+      treeHeight++;
+      Node root = new Node(rootNodeId, treeHeight, maxNodeEntries);
+      root.addEntry(newNode.mbr, newNode.nodeId);
+      root.addEntry(oldRoot.mbr, oldRoot.nodeId);
+      nodeMap.put(rootNodeId, root);
+    }
+    
+    if (INTERNAL_CONSISTENCY_CHECKING) {
+      checkConsistency(rootNodeId, treeHeight, null);
+    }
+  } 
+  
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#delete(Rectangle, int)
+   */
+  public boolean delete(Rectangle r, int id) {
+   // FindLeaf algorithm inlined here. Note the "official" algorithm 
+   // searches all overlapping entries. This seems inefficient to me, 
+   // as an entry is only worth searching if it contains (NOT overlaps)
+   // the rectangle we are searching for.
+   //
+   // Also the algorithm has been changed so that it is not recursive.
+    
+    // FL1 [Search subtrees] If root is not a leaf, check each entry 
+    // to determine if it contains r. For each entry found, invoke
+    // findLeaf on the node pointed to by the entry, until r is found or
+    // all entries have been checked.
+       parents.clear();
+       parents.push(rootNodeId);
+       
+       parentsEntry.clear();
+       parentsEntry.push(-1);
+       Node n = null;
+       int foundIndex = -1;  // index of entry to be deleted in leaf
+       
+       while (foundIndex == -1 && parents.size() > 0) {
+         n = getNode(parents.peek());
+         int startIndex = parentsEntry.peek() + 1;
+      
+      if (!n.isLeaf()) {
+          /*
+            //deleteLog.debug("searching node " + n.nodeId + ", from index " + startIndex);
+          */
+                 boolean contains = false;
+        for (int i = startIndex; i < n.entryCount; i++) {
+                   if (n.entries[i].contains(r)) {
+                     parents.push(n.ids[i]);
+                     parentsEntry.pop();
+                     parentsEntry.push(i); // this becomes the start index when the child has been searched
+                     parentsEntry.push(-1);
+                     contains = true;
+            break; // ie go to next iteration of while()
+                   }
+                 }
+        if (contains) {
+          continue;
+        }
+      } else {
+        foundIndex = n.findEntry(r, id);        
+      }
+      
+      parents.pop();
+      parentsEntry.pop();
+       } // while not found
+       
+       if (foundIndex != -1) {
+         n.deleteEntry(foundIndex, minNodeEntries);
+      condenseTree(n);
+      size--;
+       }
+    
+    return (foundIndex != -1);
+  }
+  
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#nearest(Point, IntProcedure, float)
+   */
+  public void nearest(Point p, IntProcedure v, float furthestDistance) {
+    Node rootNode = getNode(rootNodeId);
+   
+    nearest(p, rootNode, furthestDistance);
+   
+    visitProc.setProcedure(v);
+    nearestIds.forEach(visitProc);
+    nearestIds.clear();
+  }
+   
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#intersects(Rectangle, IntProcedure)
+   */
+  public void intersects(Rectangle r, IntProcedure v) {
+    Node rootNode = getNode(rootNodeId);
+    intersects(r, v, rootNode);
+  }
+
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#contains(Rectangle, IntProcedure)
+   */
+  public void contains(Rectangle r, IntProcedure v) {
+    // find all rectangles in the tree that are contained by the passed rectangle
+    // written to be non-recursive (should model other searches on this?)
+        
+    parents.clear();
+    parents.push(rootNodeId);
+    
+    parentsEntry.clear();
+    parentsEntry.push(-1);
+    
+    // TODO: possible shortcut here - could test for intersection with the 
+    // MBR of the root node. If no intersection, return immediately.
+    
+    while (parents.size() > 0) {
+      Node n = getNode(parents.peek());
+      int startIndex = parentsEntry.peek() + 1;
+      
+      if (!n.isLeaf()) {
+        // go through every entry in the index node to check
+        // if it intersects the passed rectangle. If so, it 
+        // could contain entries that are contained.
+        boolean intersects = false;
+        for (int i = startIndex; i < n.entryCount; i++) {
+          if (r.intersects(n.entries[i])) {
+            parents.push(n.ids[i]);
+            parentsEntry.pop();
+            parentsEntry.push(i); // this becomes the start index when the child has been searched
+            parentsEntry.push(-1);
+            intersects = true;
+            break; // ie go to next iteration of while()
+          }
+        }
+        if (intersects) {
+          continue;
+        }
+      } else {
+        // go through every entry in the leaf to check if 
+        // it is contained by the passed rectangle
+        for (int i = 0; i < n.entryCount; i++) {
+          if (r.contains(n.entries[i])) {
+            v.execute(n.ids[i]);
+          } 
+        }                       
+      }
+      parents.pop();
+      parentsEntry.pop();  
+    }
+  }
+
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#size()
+   */
+  public int size() {
+    return size;
+  }
+
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#getBounds()
+   */
+  public Rectangle getBounds() {
+    Rectangle bounds = null;
+    
+    Node n = getNode(getRootNodeId());
+    if (n != null && n.getMBR() != null) {
+      bounds = n.getMBR().copy();
+    }
+    return bounds;
+  }
+    
+  /**
+   * @see com.infomatiq.jsi.SpatialIndex#getVersion()
+   */
+  public String getVersion() {
+    return "RTree-" + version;
+  }
+  //-------------------------------------------------------------------------
+  // end of SpatialIndex methods
+  //-------------------------------------------------------------------------
+  
+  /**
+   * Get the next available node ID. Reuse deleted node IDs if
+   * possible
+   */
+  private int getNextNodeId() {
+    int nextNodeId = 0;
+    if (deletedNodeIds.size() > 0) {
+      nextNodeId = deletedNodeIds.pop();
+    } else {
+      nextNodeId = 1 + highestUsedNodeId++;
+    }
+    return nextNodeId;
+  }
+
+  /**
+   * Get a node object, given the ID of the node.
+   */
+  public Node getNode(int index) {
+    return (Node) nodeMap.get(index);
+  }
+
+  /**
+   * Get the highest used node ID
+   */  
+  public int getHighestUsedNodeId() {
+    return highestUsedNodeId;
+  }
+
+  /**
+   * Get the root node ID
+   */
+  public int getRootNodeId() {
+    return rootNodeId; 
+  }
+      
+  /**
+   * Split a node. Algorithm is taken pretty much verbatim from
+   * Guttman's original paper.
+   * 
+   * @return new node object.
+   */
+  private Node splitNode(Node n, Rectangle newRect, int newId) {
+    // [Pick first entry for each group] Apply algorithm pickSeeds to 
+    // choose two entries to be the first elements of the groups. Assign
+    // each to a group.
+    
+    // debug code
+    float initialArea = 0;
+    /*
+    if (log.isDebugEnabled()) {
+      Rectangle union = n.mbr.union(newRect);
+      initialArea = union.area();
+    }
+    */ 
+    System.arraycopy(initialEntryStatus, 0, entryStatus, 0, maxNodeEntries);
+    
+    Node newNode = null;
+    newNode = new Node(getNextNodeId(), n.level, maxNodeEntries);
+    nodeMap.put(newNode.nodeId, newNode);
+    
+    pickSeeds(n, newRect, newId, newNode); // this also sets the entryCount to 1
+    
+    // [Check if done] If all entries have been assigned, stop. If one
+    // group has so few entries that all the rest must be assigned to it in 
+    // order for it to have the minimum number m, assign them and stop. 
+    while (n.entryCount + newNode.entryCount < maxNodeEntries + 1) {
+      if (maxNodeEntries + 1 - newNode.entryCount == minNodeEntries) {
+        // assign all remaining entries to original node
+        for (int i = 0; i < maxNodeEntries; i++) {
+          if (entryStatus[i] == ENTRY_STATUS_UNASSIGNED) {
+            entryStatus[i] = ENTRY_STATUS_ASSIGNED;
+            n.mbr.add(n.entries[i]);
+            n.entryCount++;
+          }
+        }
+        break;
+      }   
+      if (maxNodeEntries + 1 - n.entryCount == minNodeEntries) {
+        // assign all remaining entries to new node
+        for (int i = 0; i < maxNodeEntries; i++) {
+          if (entryStatus[i] == ENTRY_STATUS_UNASSIGNED) {
+            entryStatus[i] = ENTRY_STATUS_ASSIGNED;
+            newNode.addEntryNoCopy(n.entries[i], n.ids[i]);
+            n.entries[i] = null;
+          }
+        }
+        break;
+      }
+      
+      // [Select entry to assign] Invoke algorithm pickNext to choose the
+      // next entry to assign. Add it to the group whose covering rectangle 
+      // will have to be enlarged least to accommodate it. Resolve ties
+      // by adding the entry to the group with smaller area, then to the 
+      // the one with fewer entries, then to either. Repeat from S2
+      pickNext(n, newNode);   
+    }
+      
+    n.reorganize(this);
+    
+    // check that the MBR stored for each node is correct.
+    if (INTERNAL_CONSISTENCY_CHECKING) {
+      if (!n.mbr.equals(calculateMBR(n))) {
+          throw new Error("Error: splitNode old node MBR wrong");
+      }
+      
+      if (!newNode.mbr.equals(calculateMBR(newNode))) {
+        throw new Error("Error: splitNode new node MBR wrong");
+      }
+    }
+    
+    // debug code
+    /*
+    if (log.isDebugEnabled()) {
+      float newArea = n.mbr.area() + newNode.mbr.area();
+      float percentageIncrease = (100 * (newArea - initialArea)) / initialArea;
+      //log.debug("Node " + n.nodeId + " split. New area increased by " + percentageIncrease + "%");   
+    }
+    */
+    return newNode;
+  }
+  
+  /**
+   * Pick the seeds used to split a node.
+   * Select two entries to be the first elements of the groups
+   */
+  private void pickSeeds(Node n, Rectangle newRect, int newId, Node newNode) {
+    // Find extreme rectangles along all dimension. Along each dimension,
+    // find the entry whose rectangle has the highest low side, and the one 
+    // with the lowest high side. Record the separation.
+    float maxNormalizedSeparation = 0;
+    int highestLowIndex = 0;
+    int lowestHighIndex = 0;
+    
+    // for the purposes of picking seeds, take the MBR of the node to include
+    // the new rectangle aswell.
+    n.mbr.add(newRect);
+    
+    /*
+    if (log.isDebugEnabled()) {
+        //log.debug("pickSeeds(): NodeId = " + n.nodeId + ", newRect = " + newRect);
+    }
+    */
+    
+    for (int d = 0; d < Rectangle.DIMENSIONS; d++) {
+      float tempHighestLow = newRect.min[d];
+      int tempHighestLowIndex = -1; // -1 indicates the new rectangle is the seed
+      
+      float tempLowestHigh = newRect.max[d];
+      int tempLowestHighIndex = -1;
+      
+      for (int i = 0; i < n.entryCount; i++) {
+        float tempLow = n.entries[i].min[d];
+        if (tempLow >= tempHighestLow) {
+           tempHighestLow = tempLow;
+           tempHighestLowIndex = i;
+        } else {  // ensure that the same index cannot be both lowestHigh and highestLow
+          float tempHigh = n.entries[i].max[d];
+          if (tempHigh <= tempLowestHigh) {
+            tempLowestHigh = tempHigh;
+            tempLowestHighIndex = i;
+          }
+        }
+      
+        // PS2 [Adjust for shape of the rectangle cluster] Normalize the separations
+        // by dividing by the widths of the entire set along the corresponding
+        // dimension
+        float normalizedSeparation = (tempHighestLow - tempLowestHigh) / (n.mbr.max[d] - n.mbr.min[d]);
+        
+        if (normalizedSeparation > 1 || normalizedSeparation < -1) {
+          throw new Error("Invalid normalized separation");
+        }
+      
+        /*
+        if (log.isDebugEnabled()) {
+            log.debug("Entry " + i + ", dimension " + d + ": HighestLow = " + tempHighestLow + 
+                    " (index " + tempHighestLowIndex + ")" + ", LowestHigh = " +
+                    tempLowestHigh + " (index " + tempLowestHighIndex + ", NormalizedSeparation = " + normalizedSeparation);
+        }
+        */
+        
+        // PS3 [Select the most extreme pair] Choose the pair with the greatest
+        // normalized separation along any dimension.
+        if (normalizedSeparation > maxNormalizedSeparation) {
+          maxNormalizedSeparation = normalizedSeparation;
+          highestLowIndex = tempHighestLowIndex;
+          lowestHighIndex = tempLowestHighIndex;
+        }
+      }
+    }
+      
+    // highestLowIndex is the seed for the new node.
+    if (highestLowIndex == -1) {
+      newNode.addEntry(newRect, newId);
+    } else {
+      newNode.addEntryNoCopy(n.entries[highestLowIndex], n.ids[highestLowIndex]);
+      n.entries[highestLowIndex] = null;
+      
+      // move the new rectangle into the space vacated by the seed for the new node
+      n.entries[highestLowIndex] = newRect;
+      n.ids[highestLowIndex] = newId;
+    }
+    
+    // lowestHighIndex is the seed for the original node. 
+    if (lowestHighIndex == -1) {
+      lowestHighIndex = highestLowIndex;
+    }
+    
+    entryStatus[lowestHighIndex] = ENTRY_STATUS_ASSIGNED;
+    n.entryCount = 1;
+    n.mbr.set(n.entries[lowestHighIndex].min, n.entries[lowestHighIndex].max);
+  }
+
+  /** 
+   * Pick the next entry to be assigned to a group during a node split.
+   * 
+   * [Determine cost of putting each entry in each group] For each 
+   * entry not yet in a group, calculate the area increase required
+   * in the covering rectangles of each group  
+   */
+  private int pickNext(Node n, Node newNode) {
+    float maxDifference = Float.NEGATIVE_INFINITY;
+    int next = 0;
+    int nextGroup = 0;
+    
+    maxDifference = Float.NEGATIVE_INFINITY;
+   
+    /*
+    if (log.isDebugEnabled()) {
+        //log.debug("pickNext()");
+    }
+    */
+   
+    for (int i = 0; i < maxNodeEntries; i++) {
+      if (entryStatus[i] == ENTRY_STATUS_UNASSIGNED) {
+        
+        if (n.entries[i] == null) {
+          throw new Error("Error: Node " + n.nodeId + ", entry " + i + " is null");
+        }
+        
+        float nIncrease = n.mbr.enlargement(n.entries[i]);
+        float newNodeIncrease = newNode.mbr.enlargement(n.entries[i]);
+        float difference = Math.abs(nIncrease - newNodeIncrease);
+         
+        if (difference > maxDifference) {
+          next = i;
+          
+          if (nIncrease < newNodeIncrease) {
+            nextGroup = 0; 
+          } else if (newNodeIncrease < nIncrease) {
+            nextGroup = 1;
+          } else if (n.mbr.area() < newNode.mbr.area()) {
+            nextGroup = 0;
+          } else if (newNode.mbr.area() < n.mbr.area()) {
+            nextGroup = 1;
+          } else if (newNode.entryCount < maxNodeEntries / 2) {
+            nextGroup = 0;
+          } else {
+            nextGroup = 1;
+          }
+          maxDifference = difference; 
+        }
+        /*
+        if (log.isDebugEnabled()) {
+            log.debug("Entry " + i + " group0 increase = " + nIncrease + ", group1 increase = " + newNodeIncrease +
+                    ", diff = " + difference + ", MaxDiff = " + maxDifference + " (entry " + next + ")");
+        }
+        */
+      }
+    }
+    
+    entryStatus[next] = ENTRY_STATUS_ASSIGNED;
+      
+    if (nextGroup == 0) {
+      n.mbr.add(n.entries[next]);
+      n.entryCount++;
+    } else {
+      // move to new node.
+      newNode.addEntryNoCopy(n.entries[next], n.ids[next]);
+      n.entries[next] = null;
+    }
+    
+    return next; 
+  }
+
+  /**
+   * Recursively searches the tree for the nearest entry. Other queries
+   * call execute() on an IntProcedure when a matching entry is found; 
+   * however nearest() must store the entry Ids as it searches the tree,
+   * in case a nearer entry is found.
+   * Uses the member variable nearestIds to store the nearest
+   * entry IDs.
+   * 
+   * [x] TODO rewrite this to be non-recursive?
+   */
+  private float nearest(Point p, Node n, float nearestDistance) {
+    for (int i = 0; i < n.entryCount; i++) {
+      float tempDistance = n.entries[i].distance(p);  
+      if (n.isLeaf()) { // for leaves, the distance is an actual nearest distance 
+        if (tempDistance < nearestDistance) {
+          nearestDistance = tempDistance;
+          nearestIds.clear();
+        }
+        if (tempDistance <= nearestDistance) {
+          nearestIds.add(n.ids[i]);
+        }     
+      } else { // for index nodes, only go into them if they potentially could have
+               // a rectangle nearer than actualNearest
+         if (tempDistance <= nearestDistance) {
+           // search the child node
+           nearestDistance = nearest(p, getNode(n.ids[i]), nearestDistance);
+         }
+      }
+    }
+    return nearestDistance;
+  }
+  
+  /** 
+   * Recursively searches the tree for all intersecting entries.
+   * Immediately calls execute() on the passed IntProcedure when 
+   * a matching entry is found.
+   * 
+   * [x] TODO rewrite this to be non-recursive? Make sure it
+   * doesn't slow it down.
+   */
+  private void intersects(Rectangle r, IntProcedure v, Node n) {
+    for (int i = 0; i < n.entryCount; i++) {
+      if (r.intersects(n.entries[i])) {
+        if (n.isLeaf()) {
+          v.execute(n.ids[i]);
+        } else {
+          Node childNode = getNode(n.ids[i]);
+          intersects(r, v, childNode);
+        }
+      }
+    }
+  }
+
+  /**
+   * Used by delete(). Ensures that all nodes from the passed node
+   * up to the root have the minimum number of entries.
+   * 
+   * Note that the parent and parentEntry stacks are expected to
+   * contain the nodeIds of all parents up to the root.
+   */
+    private Rectangle oldRectangle = new Rectangle(0, 0, 0, 0, 0, 0);
+  private void condenseTree(Node l) {
+    // CT1 [Initialize] Set n=l. Set the list of eliminated
+    // nodes to be empty.
+    Node n = l;
+    Node parent = null;
+    int parentEntry = 0;
+    
+    TIntStack eliminatedNodeIds = new TIntStack();
+  
+    // CT2 [Find parent entry] If N is the root, go to CT6. Otherwise 
+    // let P be the parent of N, and let En be N's entry in P  
+    while (n.level != treeHeight) {
+      parent = getNode(parents.pop());
+      parentEntry = parentsEntry.pop();
+      
+      // CT3 [Eliminiate under-full node] If N has too few entries,
+      // delete En from P and add N to the list of eliminated nodes
+      if (n.entryCount < minNodeEntries) {
+        parent.deleteEntry(parentEntry, minNodeEntries);
+        eliminatedNodeIds.push(n.nodeId);
+      } else {
+        // CT4 [Adjust covering rectangle] If N has not been eliminated,
+        // adjust EnI to tightly contain all entries in N
+        if (!n.mbr.equals(parent.entries[parentEntry])) {
+          oldRectangle.set(parent.entries[parentEntry].min, parent.entries[parentEntry].max);
+          parent.entries[parentEntry].set(n.mbr.min, n.mbr.max);
+          parent.recalculateMBR(oldRectangle);
+        }
+      }
+      // CT5 [Move up one level in tree] Set N=P and repeat from CT2
+      n = parent;
+    }
+    
+    // CT6 [Reinsert orphaned entries] Reinsert all entries of nodes in set Q.
+    // Entries from eliminated leaf nodes are reinserted in tree leaves as in 
+    // Insert(), but entries from higher level nodes must be placed higher in 
+    // the tree, so that leaves of their dependent subtrees will be on the same
+    // level as leaves of the main tree
+    while (eliminatedNodeIds.size() > 0) {
+      Node e = getNode(eliminatedNodeIds.pop());
+      for (int j = 0; j < e.entryCount; j++) {
+        add(e.entries[j], e.ids[j], e.level); 
+        e.entries[j] = null;
+      }
+      e.entryCount = 0;
+      deletedNodeIds.push(e.nodeId);
+    }
+  }
+
+  /**
+   *  Used by add(). Chooses a leaf to add the rectangle to.
+   */
+  private Node chooseNode(Rectangle r, int level) {
+    // CL1 [Initialize] Set N to be the root node
+    Node n = getNode(rootNodeId);
+    parents.clear();
+    parentsEntry.clear();
+     
+    // CL2 [Leaf check] If N is a leaf, return N
+    while (true) {
+      if (n == null) {
+        throw new Error("Could not get root node (" + rootNodeId + ")");  
+      }
+   
+      if (n.level == level) {
+        return n;
+      }
+      
+      // CL3 [Choose subtree] If N is not at the desired level, let F be the entry in N 
+      // whose rectangle FI needs least enlargement to include EI. Resolve
+      // ties by choosing the entry with the rectangle of smaller area.
+      float leastEnlargement = n.getEntry(0).enlargement(r);
+      int index = 0; // index of rectangle in subtree
+      for (int i = 1; i < n.entryCount; i++) {
+        Rectangle tempRectangle = n.getEntry(i);
+        float tempEnlargement = tempRectangle.enlargement(r);
+        if ((tempEnlargement < leastEnlargement) ||
+            ((tempEnlargement == leastEnlargement) && 
+             (tempRectangle.area() < n.getEntry(index).area()))) {
+          index = i;
+          leastEnlargement = tempEnlargement;
+        }
+      }
+      
+      parents.push(n.nodeId);
+      parentsEntry.push(index);
+    
+      // CL4 [Descend until a leaf is reached] Set N to be the child node 
+      // pointed to by Fp and repeat from CL2
+      n = getNode(n.ids[index]);
+    }
+  }
+  
+  /**
+   * Ascend from a leaf node L to the root, adjusting covering rectangles and
+   * propagating node splits as necessary.
+   */
+  private Node adjustTree(Node n, Node nn) {
+    // AT1 [Initialize] Set N=L. If L was split previously, set NN to be 
+    // the resulting second node.
+    
+    // AT2 [Check if done] If N is the root, stop
+    while (n.level != treeHeight) {
+    
+      // AT3 [Adjust covering rectangle in parent entry] Let P be the parent 
+      // node of N, and let En be N's entry in P. Adjust EnI so that it tightly
+      // encloses all entry rectangles in N.
+      Node parent = getNode(parents.pop());
+      int entry = parentsEntry.pop(); 
+      
+      if (parent.ids[entry] != n.nodeId) {
+        throw new Error("Error: entry " + entry + " in node " + 
+             parent.nodeId + " should point to node " + 
+             n.nodeId + "; actually points to node " + parent.ids[entry]);
+      }
+      
+      if (!parent.entries[entry].equals(n.mbr)) {
+        parent.entries[entry].set(n.mbr.min, n.mbr.max);
+        parent.mbr.set(parent.entries[0].min, parent.entries[0].max);
+        for (int i = 1; i < parent.entryCount; i++) {
+          parent.mbr.add(parent.entries[i]);
+        }
+      }
+      
+      // AT4 [Propagate node split upward] If N has a partner NN resulting from 
+      // an earlier split, create a new entry Enn with Ennp pointing to NN and 
+      // Enni enclosing all rectangles in NN. Add Enn to P if there is room. 
+      // Otherwise, invoke splitNode to produce P and PP containing Enn and
+      // all P's old entries.
+      Node newNode = null;
+      if (nn != null) {
+        if (parent.entryCount < maxNodeEntries) {
+          parent.addEntry(nn.mbr, nn.nodeId);
+        } else {
+          newNode = splitNode(parent, nn.mbr.copy(), nn.nodeId);
+        }
+      }
+      
+      // AT5 [Move up to next level] Set N = P and set NN = PP if a split 
+      // occurred. Repeat from AT2
+      n = parent;
+      nn = newNode;
+      
+      parent = null;
+      newNode = null;
+    }
+    
+    return nn;
+  }
+  
+  /**
+   * Check the consistency of the tree.
+   */
+  private void checkConsistency(int nodeId, int expectedLevel, Rectangle expectedMBR) {
+    // go through the tree, and check that the internal data structures of 
+    // the tree are not corrupted.    
+    Node n = getNode(nodeId);
+    
+    if (n == null) {
+      throw new Error("Error: Could not read node " + nodeId);
+    }
+    
+    if (n.level != expectedLevel) {
+      throw new Error("Error: Node " + nodeId + ", expected level " + expectedLevel + ", actual level " + n.level);
+    }
+    
+    Rectangle calculatedMBR = calculateMBR(n);
+    
+    if (!n.mbr.equals(calculatedMBR)) {
+      throw new Error("Error: Node " + nodeId + ", calculated MBR does not equal stored MBR");
+    }
+    
+    if (expectedMBR != null && !n.mbr.equals(expectedMBR)) {
+      throw new Error("Error: Node " + nodeId + ", expected MBR (from parent) does not equal stored MBR");
+    }
+    
+    // Check for corruption where a parent entry is the same object as the child MBR
+    if (expectedMBR != null && n.mbr.sameObject(expectedMBR)) {
+      throw new Error("Error: Node " + nodeId + " MBR using same rectangle object as parent's entry");
+    }
+    
+    for (int i = 0; i < n.entryCount; i++) {
+      if (n.entries[i] == null) {
+        throw new Error("Error: Node " + nodeId + ", Entry " + i + " is null");
+      }     
+      
+      if (n.level > 1) { // if not a leaf
+        checkConsistency(n.ids[i], n.level - 1, n.entries[i]); 
+      }   
+    } 
+  }
+  
+  /**
+   * Given a node object, calculate the node MBR from it's entries.
+   * Used in consistency checking
+   */
+  private Rectangle calculateMBR(Node n) {
+    Rectangle mbr = new Rectangle(n.entries[0].min, n.entries[0].max);
+    
+    for (int i = 1; i < n.entryCount; i++) {
+      mbr.add(n.entries[i]);
+    }
+    return mbr; 
+  }
+}
diff --git a/src/edu/berkeley/qfat/geom/HasBoundingBox.java b/src/edu/berkeley/qfat/geom/HasBoundingBox.java
new file mode 100644 (file)
index 0000000..3c239fe
--- /dev/null
@@ -0,0 +1,11 @@
+package edu.berkeley.qfat.geom;
+import javax.media.opengl.*;
+
+public interface HasBoundingBox {
+    public float getMaxX();
+    public float getMinX();
+    public float getMaxY();
+    public float getMinY();
+    public float getMaxZ();
+    public float getMinZ();
+}
index e72b0d0..f1914e3 100644 (file)
@@ -2,6 +2,12 @@ package edu.berkeley.qfat.geom;
 import javax.media.opengl.*;
 
 /** point in 3-space; immutable */
 import javax.media.opengl.*;
 
 /** point in 3-space; immutable */
-public abstract class HasPoint {
+public abstract class HasPoint implements HasBoundingBox {
     public abstract Point getPoint();
     public abstract Point getPoint();
+    public float getMaxX() { return getPoint().getMaxX(); }
+    public float getMinX() { return getPoint().getMinX(); }
+    public float getMaxY() { return getPoint().getMaxY(); }
+    public float getMinY() { return getPoint().getMinY(); }
+    public float getMaxZ() { return getPoint().getMaxZ(); }
+    public float getMinZ() { return getPoint().getMinZ(); }
 }
 }
diff --git a/src/edu/berkeley/qfat/geom/IntervalTree.java b/src/edu/berkeley/qfat/geom/IntervalTree.java
new file mode 100644 (file)
index 0000000..9d7ed75
--- /dev/null
@@ -0,0 +1,90 @@
+package edu.berkeley.qfat.geom;
+import javax.media.opengl.*;
+
+public class IntervalTree<V extends HasBoundingBox> {
+
+
+    private static final int X_AXIS = 0;
+    private static final int Y_AXIS = 1;
+    private static final int Z_AXIS = 2;
+
+    private class Node {
+        int axis; /* 0, 1, 2 */
+        V vthis;
+        float maxmax, minmin;
+        Node small, large;
+        float bias = 0;
+
+        public Node(V vthis, int axis) {
+            this.vthis = vthis;
+            this.axis = axis;
+            maxmax = getMax(vthis);
+            minmin = getMin(vthis);
+        }
+
+        private float getMin(V v) {
+            switch(axis) {
+                case X_AXIS: return v.getMinX();
+                case Y_AXIS: return v.getMinY();
+                case Z_AXIS: return v.getMinZ();
+            }
+            throw new Error();
+        }
+        private float getMax(V v) {
+            switch(axis) {
+                case X_AXIS: return v.getMaxX();
+                case Y_AXIS: return v.getMaxY();
+                case Z_AXIS: return v.getMaxZ();
+            }
+            throw new Error();
+        }
+        public Node insert(V vnew) {
+            // FIXME bias
+            if (getMax(vnew) > maxmax && getMin(vnew) < minmin) {
+                float diff = Math.abs(maxmax - getMax(vnew)) - Math.abs(minmin - getMin(vnew));
+                if (diff < 0) {
+                    maxmax = getMax(vnew);
+                } else {
+                    minmin = getMin(vnew);
+                }
+            }
+            if (getMax(vnew) <= maxmax) {
+                right++;
+                small = small == null ? new Node(vnew, (axis+1)%3) : small.insert(vnew);
+                bias = (left - right) / (float)(left + right);
+                if (Math.abs(bias) > Math.abs(maxbias)) maxbias = bias;
+                return this;
+
+            } else if (getMin(vnew) >= minmin) {
+                left++;
+                large = large == null ? new Node(vnew, (axis+1)%3) : large.insert(vnew);
+                bias = (left - right) / (float)(left + right);
+                if (Math.abs(bias) > Math.abs(maxbias)) maxbias = bias;
+                return this;
+            }
+
+            throw new Error();
+        }
+        int left=0, right=0;
+
+    }
+
+    private Node root;
+
+    public void insert(V v) {
+        if (root==null) {
+            root = new Node(v, X_AXIS);
+        } else {
+            root.insert(v);
+        }
+        System.out.println(maxbias + " " + size);
+        size++;
+    }
+    int size = 0;
+    float maxbias = 0;
+
+    public V nearest(Point p) {
+        return null;
+    }
+
+}
index d54e631..5511e44 100644 (file)
@@ -2,7 +2,7 @@ package edu.berkeley.qfat.geom;
 import javax.media.opengl.*;
 
 /** point in 3-space; immutable */
 import javax.media.opengl.*;
 
 /** point in 3-space; immutable */
-public final class Point {
+public final class Point implements HasBoundingBox {
     public final float x, y, z;
     public Point(double x, double y, double z) { this((float)x, (float)y, (float)z); }
     public Point(float x, float y, float z) { this.x = x; this.y = y; this.z = z; }
     public final float x, y, z;
     public Point(double x, double y, double z) { this((float)x, (float)y, (float)z); }
     public Point(float x, float y, float z) { this.x = x; this.y = y; this.z = z; }
@@ -16,4 +16,11 @@ public final class Point {
     private void _glVertex(GL gl) { gl.glVertex3f(x, y, z); }
     public String toString() { return "("+x+","+y+","+z+")"; }
     public int hashCode() { return Float.floatToIntBits(x) ^ Float.floatToIntBits(y) ^ Float.floatToIntBits(z); }
     private void _glVertex(GL gl) { gl.glVertex3f(x, y, z); }
     public String toString() { return "("+x+","+y+","+z+")"; }
     public int hashCode() { return Float.floatToIntBits(x) ^ Float.floatToIntBits(y) ^ Float.floatToIntBits(z); }
+
+    public float getMaxX() { return x; }
+    public float getMinX() { return x; }
+    public float getMaxY() { return y; }
+    public float getMinY() { return y; }
+    public float getMaxZ() { return z; }
+    public float getMinZ() { return z; }
 }
 }
index cd3edb7..fe44993 100644 (file)
@@ -4,6 +4,8 @@ import java.util.*;
 
 public class PointSet<V extends HasPoint> implements Iterable<V> {
 
 
 public class PointSet<V extends HasPoint> implements Iterable<V> {
 
+    private final RTree<V> rtree = new RTree<V>();
+
     private /*final*/ KDTree kd = new KDTree(3);
     private final double[] doubles = new double[3];
 
     private /*final*/ KDTree kd = new KDTree(3);
     private final double[] doubles = new double[3];
 
@@ -23,10 +25,12 @@ public class PointSet<V extends HasPoint> implements Iterable<V> {
     }
 
     public void rebuild() {
     }
 
     public void rebuild() {
+        /*
         HashMap<Point,V> old_exact = exact;
         exact = new HashMap<Point,V>();
         kd = new KDTree(3);
         for(V v : old_exact.values()) add(v);
         HashMap<Point,V> old_exact = exact;
         exact = new HashMap<Point,V>();
         kd = new KDTree(3);
         for(V v : old_exact.values()) add(v);
+        */
     }
 
     public void add(V v) {
     }
 
     public void add(V v) {
@@ -34,6 +38,7 @@ public class PointSet<V extends HasPoint> implements Iterable<V> {
         if (x != null && x.equals(v)) return;
         if (x != null) throw new Error("duplicates!");
         Point p = v.getPoint();
         if (x != null && x.equals(v)) return;
         if (x != null) throw new Error("duplicates!");
         Point p = v.getPoint();
+        /*
         doubles[0] = p.x;
         doubles[1] = p.y;
         doubles[2] = p.z;
         doubles[0] = p.x;
         doubles[1] = p.y;
         doubles[2] = p.z;
@@ -42,22 +47,28 @@ public class PointSet<V extends HasPoint> implements Iterable<V> {
         } catch (Exception e) {
             throw new Error(e);
         }
         } catch (Exception e) {
             throw new Error(e);
         }
+        */
+        rtree.insert(v);
         exact.put(p, v);
     }
 
         exact.put(p, v);
     }
 
-    public void remove(HasPoint v) { remove(v.getPoint()); }
-    public void remove(Point p) {
+    public void remove(V v) {
+        Point p = v.getPoint();
+        /*
         doubles[0] = p.x;
         doubles[1] = p.y;
         doubles[2] = p.z;
         try {
             kd.delete(doubles);
         } catch (Exception e) { }
         doubles[0] = p.x;
         doubles[1] = p.y;
         doubles[2] = p.z;
         try {
             kd.delete(doubles);
         } catch (Exception e) { }
+        */
+        rtree.remove(v);
         exact.remove(p);
     }
 
     public V nearest(Point p) {
         if (exact.size()==0) return null;
         exact.remove(p);
     }
 
     public V nearest(Point p) {
         if (exact.size()==0) return null;
+        /*
         Object[] results;
         try {
             doubles[0] = p.x;
         Object[] results;
         try {
             doubles[0] = p.x;
@@ -67,7 +78,11 @@ public class PointSet<V extends HasPoint> implements Iterable<V> {
         } catch (Exception e) {
             throw new Error(e);
         }
         } catch (Exception e) {
             throw new Error(e);
         }
-        return (V)results[0];
+        V kd_says = (V)results[0];
+        */
+        V rt_says = rtree.nearest(p);
+        //if (kd_says != rt_says) System.err.println("disagree: " + p + " " + kd_says + " " + rt_says);
+        return rt_says;
     }
 
 
     }
 
 
diff --git a/src/edu/berkeley/qfat/geom/RTree.java b/src/edu/berkeley/qfat/geom/RTree.java
new file mode 100644 (file)
index 0000000..48fade2
--- /dev/null
@@ -0,0 +1,57 @@
+package edu.berkeley.qfat.geom;
+import javax.media.opengl.*;
+import java.util.*;
+import com.infomatiq.jsi.*;
+import com.infomatiq.jsi.rtree.*;
+
+public class RTree<V extends HasBoundingBox> {
+
+    private com.infomatiq.jsi.rtree.RTree rtree =
+        new com.infomatiq.jsi.rtree.RTree();
+
+    int lowid = 0;
+    HashMap<Integer, V> idToV = new HashMap<Integer, V>();
+    HashMap<V, Integer> vToId = new HashMap<V, Integer>();
+
+    public RTree() {
+        Properties props = new Properties();
+        props.put("MinNodeEntries", "1");
+        props.put("MaxNodeEntries", "5");
+        rtree.init(props);
+    }
+
+    public void insert(V v) {
+        int id = lowid++;
+        idToV.put(id, v);
+        vToId.put(v, id);
+        rtree.add(new com.infomatiq.jsi.Rectangle(v.getMinX(), v.getMinY(), v.getMinZ(),
+                                                  v.getMaxX(), v.getMaxY(), v.getMaxZ()),
+                  id);
+    }
+
+    public void remove(V v) {
+        int id = vToId.get(v);
+        idToV.remove(id);
+        vToId.remove(v);
+        rtree.delete(new com.infomatiq.jsi.Rectangle(v.getMinX(), v.getMinY(), v.getMinZ(),
+                                                     v.getMaxX(), v.getMaxY(), v.getMaxZ()),
+                     id);
+    }
+
+
+    V found = null;
+
+    private IntProcedure finder = new IntProcedure() {
+            public boolean execute(int id) {
+                found = idToV.get(id);
+                return false;
+            }
+        };
+
+    public V nearest(Point p) {
+        rtree.nearest(new com.infomatiq.jsi.Point(p.x, p.y, p.z), finder, Float.POSITIVE_INFINITY);
+        V ret = found;
+        found = null;
+        return ret;
+    }
+}
index 6a938a3..df7c463 100644 (file)
@@ -1,7 +1,7 @@
 package edu.berkeley.qfat.geom;
 import javax.media.opengl.*;
 
 package edu.berkeley.qfat.geom;
 import javax.media.opengl.*;
 
-public abstract class Triangle {
+public abstract class Triangle implements HasBoundingBox {
     public abstract Point p1();
     public abstract Point p2();
     public abstract Point p3();
     public abstract Point p1();
     public abstract Point p2();
     public abstract Point p3();
@@ -44,4 +44,10 @@ public abstract class Triangle {
         return (area()/(max*max));
     }
 
         return (area()/(max*max));
     }
 
+    public float getMaxX() { return Math.max(p1().x, Math.max(p2().x, p3().x)); }
+    public float getMinX() { return Math.min(p1().x, Math.min(p2().x, p3().x)); }
+    public float getMaxY() { return Math.max(p1().y, Math.max(p2().y, p3().y)); }
+    public float getMinY() { return Math.min(p1().y, Math.min(p2().y, p3().y)); }
+    public float getMaxZ() { return Math.max(p1().z, Math.max(p2().z, p3().z)); }
+    public float getMinZ() { return Math.min(p1().z, Math.min(p2().z, p3().z)); }
 }
\ No newline at end of file
 }
\ No newline at end of file