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 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): + *

+ * + *

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: + *

+ * + * @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