From: adam
Date: Fri, 14 Dec 2007 01:53:17 +0000 (-0800)
Subject: add RTree classes
X-Git-Url: http://git.megacz.com/?p=anneal.git;a=commitdiff_plain;h=0f9ce20a060db6537a47b549cbf24fd268699ac6;hp=57376a862c00fa1c8731f9989085fcfceeee0370
add RTree classes
darcs-hash:20071214015317-5007d-31f9de5bc2afdc9cade724105114fec68e940c5b.gz
---
diff --git a/Makefile b/Makefile
index dec73e8..d42a2b7 100644
--- 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
- 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
index 0000000..6081136
--- /dev/null
+++ b/src/com/infomatiq/jsi/IntProcedure.java
@@ -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
index 0000000..c608971
--- /dev/null
+++ b/src/com/infomatiq/jsi/Point.java
@@ -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
index 0000000..737e6cf
--- /dev/null
+++ b/src/com/infomatiq/jsi/Rectangle.java
@@ -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
index 0000000..a385cf5
--- /dev/null
+++ b/src/com/infomatiq/jsi/SpatialIndex.java
@@ -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
index 0000000..ac00dc0
--- /dev/null
+++ b/src/com/infomatiq/jsi/rtree/Node.java
@@ -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;
+
+/**
+ * Used by RTree. There are no public methods in this class.
+ *
+ * @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
index 0000000..ba9b130
--- /dev/null
+++ b/src/com/infomatiq/jsi/rtree/RTree.java
@@ -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;
+
+/**
+ * This is a lightweight RTree implementation, specifically designed
+ * for the following features (in order of importance):
+ *
+ * - 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.
+ * - Low memory requirements.
+ * - Fast add performance.
+ *
+ *
+ * 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.
+ *
+ * @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()
+ //-------------------------------------------------------------------------
+ /**
+ * Initialize implementation dependent properties of the RTree.
+ * Currently implemented properties are:
+ *
+ * - MaxNodeEntries
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.
+ * - MinNodeEntries
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.
+ *
+ *
+ * @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
index 0000000..3c239fe
--- /dev/null
+++ b/src/edu/berkeley/qfat/geom/HasBoundingBox.java
@@ -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();
+}
diff --git a/src/edu/berkeley/qfat/geom/HasPoint.java b/src/edu/berkeley/qfat/geom/HasPoint.java
index e72b0d0..f1914e3 100644
--- a/src/edu/berkeley/qfat/geom/HasPoint.java
+++ b/src/edu/berkeley/qfat/geom/HasPoint.java
@@ -2,6 +2,12 @@ package edu.berkeley.qfat.geom;
import javax.media.opengl.*;
/** point in 3-space; immutable */
-public abstract class HasPoint {
+public abstract class HasPoint implements HasBoundingBox {
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
index 0000000..9d7ed75
--- /dev/null
+++ b/src/edu/berkeley/qfat/geom/IntervalTree.java
@@ -0,0 +1,90 @@
+package edu.berkeley.qfat.geom;
+import javax.media.opengl.*;
+
+public class IntervalTree {
+
+
+ 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;
+ }
+
+}
diff --git a/src/edu/berkeley/qfat/geom/Point.java b/src/edu/berkeley/qfat/geom/Point.java
index d54e631..5511e44 100644
--- a/src/edu/berkeley/qfat/geom/Point.java
+++ b/src/edu/berkeley/qfat/geom/Point.java
@@ -2,7 +2,7 @@ package edu.berkeley.qfat.geom;
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; }
@@ -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); }
+
+ 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; }
}
diff --git a/src/edu/berkeley/qfat/geom/PointSet.java b/src/edu/berkeley/qfat/geom/PointSet.java
index cd3edb7..fe44993 100644
--- a/src/edu/berkeley/qfat/geom/PointSet.java
+++ b/src/edu/berkeley/qfat/geom/PointSet.java
@@ -4,6 +4,8 @@ import java.util.*;
public class PointSet implements Iterable {
+ private final RTree rtree = new RTree();
+
private /*final*/ KDTree kd = new KDTree(3);
private final double[] doubles = new double[3];
@@ -23,10 +25,12 @@ public class PointSet implements Iterable {
}
public void rebuild() {
+ /*
HashMap old_exact = exact;
exact = new HashMap();
kd = new KDTree(3);
for(V v : old_exact.values()) add(v);
+ */
}
public void add(V v) {
@@ -34,6 +38,7 @@ public class PointSet implements Iterable {
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;
@@ -42,22 +47,28 @@ public class PointSet implements Iterable {
} catch (Exception e) {
throw new Error(e);
}
+ */
+ rtree.insert(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) { }
+ */
+ rtree.remove(v);
exact.remove(p);
}
public V nearest(Point p) {
if (exact.size()==0) return null;
+ /*
Object[] results;
try {
doubles[0] = p.x;
@@ -67,7 +78,11 @@ public class PointSet implements Iterable {
} 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
index 0000000..48fade2
--- /dev/null
+++ b/src/edu/berkeley/qfat/geom/RTree.java
@@ -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 {
+
+ private com.infomatiq.jsi.rtree.RTree rtree =
+ new com.infomatiq.jsi.rtree.RTree();
+
+ int lowid = 0;
+ HashMap idToV = new HashMap();
+ HashMap vToId = new HashMap();
+
+ 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;
+ }
+}
diff --git a/src/edu/berkeley/qfat/geom/Triangle.java b/src/edu/berkeley/qfat/geom/Triangle.java
index 6a938a3..df7c463 100644
--- a/src/edu/berkeley/qfat/geom/Triangle.java
+++ b/src/edu/berkeley/qfat/geom/Triangle.java
@@ -1,7 +1,7 @@
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();
@@ -44,4 +44,10 @@ public abstract class Triangle {
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