Fix root-node precision problems

README.md

@@ -41,6 +41,9 @@ Indexes can be created via factories in the interfaces, e.g. `PointMap.Factory.c
 Note:
  - **STR-Trees** are simply R-Trees that are preloaded using the STR algorithm. THis can be done with
    the factory methods `....Factory.createAndLoadStrRTree(...)`.
+ - **Quadtree precision problems**. Many quadtree implementations are vulnerable to precision problem due to repeated 
+   division by 2.0 of radius and subnode centers. Our quadtree implementations avoid this problem by aligning node center
+   coordinates and radius to a power of two. Power of two values are mostly immune to precision problems when dividing by 2.0.
  - `PointArray` and `BoxArray` are simple array based implementations. They scale badly with size, their only use is for verifying correctness of other indexes. 
 
 ## Changelog

src/main/java/org/tinspin/index/PointMap.java

@@ -185,6 +185,7 @@ static <T> PointMap<T> createQuadtree(int dims) {
         }
 
         /**
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
          * Create a plain Quadtree.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
@@ -195,11 +196,31 @@ static <T> PointMap<T> createQuadtree(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New Quadtree
+         * @deprecated Please use {@link #createAlignedQuadtree(int, int, double[], double)}
          */
+        @Deprecated
         static <T> PointMap<T> createQuadtree(int dims, int maxNodeCapacity, double[] center, double radius) {
             return QuadTreeKD0.create(dims, maxNodeCapacity, center, radius);
         }
 
+        /**
+         * Create a plain Quadtree.
+         * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param dims            Number of dimensions.
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param center          Estimated center of all coordinates.
+         * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
+         * @param <T>             Value type
+         * @return New Quadtree
+         */
+        static <T> PointMap<T> createAlignedQuadtree(int dims, int maxNodeCapacity, double[] center, double radius) {
+            return QuadTreeKD0.createAligned(dims, maxNodeCapacity, center, radius);
+        }
+
         /**
          * Create a Quadtree with hypercube navigation.
          *
@@ -212,7 +233,8 @@ static <T> PointMap<T> createQuadtreeHC(int dims) {
         }
 
         /**
-         * Create a plain Quadtree.
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
+         * Create a Quadtree with hypercube navigation.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
          *
@@ -222,13 +244,33 @@ static <T> PointMap<T> createQuadtreeHC(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC
+         * @deprecated Please use {@link #createAlignedQuadtreeHC(int, int, double[], double)}
          */
+        @Deprecated
         static <T> PointMap<T> createQuadtreeHC(int dims, int maxNodeCapacity, double[] center, double radius) {
             return QuadTreeKD.create(dims, maxNodeCapacity, center, radius);
         }
 
         /**
          * Create a Quadtree with hypercube navigation.
+         * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param dims            Number of dimensions.
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param center          Estimated center of all coordinates.
+         * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
+         * @param <T>             Value type
+         * @return New QuadtreeHC
+         */
+        static <T> PointMap<T> createAlignedQuadtreeHC(int dims, int maxNodeCapacity, double[] center, double radius) {
+            return QuadTreeKD.createAligned(dims, maxNodeCapacity, center, radius);
+        }
+
+        /**
+         * Create a Quadtree with extended hypercube navigation.
          *
          * @param dims Number of dimensions.
          * @param <T>  Value type
@@ -239,7 +281,7 @@ static <T> PointMap<T> createQuadtreeHC2(int dims) {
         }
 
         /**
-         * Create a plain Quadtree.
+         * Create a Quadtree with extended hypercube navigation.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
          *
@@ -249,11 +291,31 @@ static <T> PointMap<T> createQuadtreeHC2(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC2
+         * @deprecated PLease use {@link #createAlignedQuadtreeHC2(int, int, double[], double)}
          */
+        @Deprecated
         static <T> PointMap<T> createQuadtreeHC2(int dims, int maxNodeCapacity, double[] center, double radius) {
             return QuadTreeKD2.create(dims, maxNodeCapacity, center, radius);
         }
 
+        /**
+         * Create a Quadtree with extended hypercube navigation.
+         * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param dims            Number of dimensions.
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param center          Estimated center of all coordinates.
+         * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
+         * @param <T>             Value type
+         * @return New QuadtreeHC2
+         */
+        static <T> PointMap<T> createAlignedQuadtreeHC2(int dims, int maxNodeCapacity, double[] center, double radius) {
+            return QuadTreeKD2.createAligned(dims, maxNodeCapacity, center, radius);
+        }
+
         /**
          * Create an R*Tree.
          *

src/main/java/org/tinspin/index/PointMultimap.java

@@ -199,6 +199,7 @@ static <T> PointMultimap<T> createQuadtree(int dims) {
         }
 
         /**
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
          * Create a plain Quadtree.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
@@ -209,11 +210,31 @@ static <T> PointMultimap<T> createQuadtree(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New Quadtree
+         * @deprecated Please use {@link #createAlignedQuadtree(int, int, double[], double)}
          */
+        @Deprecated
         static <T> PointMultimap<T> createQuadtree(int dims, int maxNodeCapacity, double[] center, double radius) {
             return QuadTreeKD0.create(dims, maxNodeCapacity, center, radius);
         }
 
+        /**
+         * Create a plain Quadtree.
+         * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param dims            Number of dimensions.
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param center          Estimated center of all coordinates.
+         * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
+         * @param <T>             Value type
+         * @return New Quadtree
+         */
+        static <T> PointMultimap<T> createAlignedQuadtree(int dims, int maxNodeCapacity, double[] center, double radius) {
+            return QuadTreeKD0.createAligned(dims, maxNodeCapacity, center, radius);
+        }
+
         /**
          * Create a Quadtree with hypercube navigation.
          *
@@ -226,7 +247,8 @@ static <T> PointMultimap<T> createQuadtreeHC(int dims) {
         }
 
         /**
-         * Create a plain Quadtree.
+         * WARNING: Unaligned center and radius can cause precision problems, see README.
+         * Create a Quadtree with hypercube navigation.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
          *
@@ -236,13 +258,33 @@ static <T> PointMultimap<T> createQuadtreeHC(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC
+         * @deprecated Please use {@link #createAlignedQuadtreeHC(int, int, double[], double)}
          */
+        @Deprecated
         static <T> PointMultimap<T> createQuadtreeHC(int dims, int maxNodeCapacity, double[] center, double radius) {
             return QuadTreeKD.create(dims, maxNodeCapacity, center, radius);
         }
 
         /**
          * Create a Quadtree with hypercube navigation.
+         * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param dims            Number of dimensions.
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param center          Estimated center of all coordinates.
+         * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
+         * @param <T>             Value type
+         * @return New QuadtreeHC
+         */
+        static <T> PointMultimap<T> createAlignedQuadtreeHC(int dims, int maxNodeCapacity, double[] center, double radius) {
+            return QuadTreeKD.createAligned(dims, maxNodeCapacity, center, radius);
+        }
+
+        /**
+         * Create a Quadtree with extended hypercube navigation.
          *
          * @param dims Number of dimensions.
          * @param <T>  Value type
@@ -253,7 +295,7 @@ static <T> PointMultimap<T> createQuadtreeHC2(int dims) {
         }
 
         /**
-         * Create a plain Quadtree.
+         * Create a Quadtree with extended hypercube navigation.
          * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
          * span an area that is somewhat larger rather than smaller than the actual data.
          *
@@ -263,11 +305,31 @@ static <T> PointMultimap<T> createQuadtreeHC2(int dims) {
          * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
          * @param <T>             Value type
          * @return New QuadtreeHC2
+         * @deprecated PLease use {@link #createAlignedQuadtreeHC2(int, int, double[], double)}
          */
+        @Deprecated
         static <T> PointMultimap<T> createQuadtreeHC2(int dims, int maxNodeCapacity, double[] center, double radius) {
             return QuadTreeKD2.create(dims, maxNodeCapacity, center, radius);
         }
 
+        /**
+         * Create a Quadtree with extended hypercube navigation.
+         * Center/radius are used to find a good initial root. They do not need to be exact. If possible, they should
+         * span an area that is somewhat larger rather than smaller than the actual data.
+         * <p>
+         * Center and radius will be aligned with powers of two to avoid precision problems.
+         *
+         * @param dims            Number of dimensions.
+         * @param maxNodeCapacity Maximum entries in a node before the node is split. The default is 10.
+         * @param center          Estimated center of all coordinates.
+         * @param radius          Estimated maximum orthogonal distance from center for all coordinates.
+         * @param <T>             Value type
+         * @return New QuadtreeHC2
+         */
+        static <T> PointMultimap<T> createAlignedQuadtreeHC2(int dims, int maxNodeCapacity, double[] center, double radius) {
+            return QuadTreeKD2.createAligned(dims, maxNodeCapacity, center, radius);
+        }
+
         /**
          * Create an R*Tree.
          *

src/main/java/org/tinspin/index/qthypercube/QuadTreeKD.java

@@ -73,9 +73,41 @@ public static <T> QuadTreeKD<T> create(int dims) {
 	public static <T> QuadTreeKD<T> create(int dims, int maxNodeSize) {
 		return new QuadTreeKD<>(dims, maxNodeSize);
 	}
-
-	public static <T> QuadTreeKD<T> create(int dims, int maxNodeSize, 
-			double[] center, double radius) {
+
+	/**
+	 * Note: This will align center and radius to a power of two before creating a tree.
+	 * @param dims dimensions, usually 2 or 3
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @param center center of initial root node
+	 * @param radius radius of initial root node
+	 * @return New quadtree
+	 * @param <T> Value type
+	 */
+	public static <T> QuadTreeKD<T> createAligned(int dims, int maxNodeSize,
+										   double[] center, double radius) {
+		QuadTreeKD<T> t = new QuadTreeKD<>(dims, maxNodeSize);
+		if (radius <= 0) {
+			throw new IllegalArgumentException("Radius must be > 0 but was " + radius);
+		}
+		double[] alignedCenter = MathTools.floorPowerOfTwoCopy(center);
+		double alignedRadius = MathTools.ceilPowerOfTwo(radius);
+		t.root = new QNode<>(Arrays.copyOf(alignedCenter, alignedCenter.length), alignedRadius);
+		return t;
+	}
+
+	/**
+	 * WARNING: Unaligned center and radius can cause precision problems.
+	 * @param dims dimensions, usually 2 or 3
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @param center center of initial root node
+	 * @param radius radius of initial root node
+	 * @return New quadtree
+	 * @param <T> Value type
+	 * @deprecated Please use {@link #createAligned(int, int, double[], double)}
+	 */
+	@Deprecated
+	public static <T> QuadTreeKD<T> create(int dims, int maxNodeSize,
+										   double[] center, double radius) {
 		QuadTreeKD<T> t = new QuadTreeKD<>(dims, maxNodeSize);
 		if (radius <= 0) {
 			throw new IllegalArgumentException("Radius must be > 0 but was " + radius);

src/main/java/org/tinspin/index/qthypercube2/QuadTreeKD2.java

@@ -65,8 +65,8 @@ public class QuadTreeKD2<T> implements PointMap<T>, PointMultimap<T> {
 	private final int dims;
 	private final int maxNodeSize;
 	private QNode<T> root = null;
-	private int size = 0; 
-	
+	private int size = 0;
+
 
 	private QuadTreeKD2(int dims, int maxNodeSize) {
 		if (DEBUG) {
@@ -76,18 +76,61 @@ private QuadTreeKD2(int dims, int maxNodeSize) {
 		this.maxNodeSize = maxNodeSize;
 	}
 
+	/**
+	 * @param dims dimensions, usually 2 or 3
+	 * @return New quadtree
+	 * @param <T> Value type
+	 */
 	public static <T> QuadTreeKD2<T> create(int dims) {
 		int maxNodeSize = DEFAULT_MAX_NODE_SIZE;
 		if (2 * dims > DEFAULT_MAX_NODE_SIZE) {
 			maxNodeSize = 2*dims;
 		}
 		return new QuadTreeKD2<>(dims, maxNodeSize);
 	}
-
+
+	/**
+	 * @param dims dimensions, usually 2 or 3
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @return New quadtree
+	 * @param <T> Value type
+	 */
 	public static <T> QuadTreeKD2<T> create(int dims, int maxNodeSize) {
 		return new QuadTreeKD2<>(dims, maxNodeSize);
 	}
-
+
+	/**
+	 * Note: This will align center and radius to a power of two before creating a tree.
+	 * @param dims dimensions, usually 2 or 3
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @param center center of initial root node
+	 * @param radius radius of initial root node
+	 * @return New quadtree
+	 * @param <T> Value type
+	 */
+	public static <T> QuadTreeKD2<T> createAligned(int dims, int maxNodeSize,
+											double[] center, double radius) {
+		QuadTreeKD2<T> t = new QuadTreeKD2<>(dims, maxNodeSize);
+		if (radius <= 0) {
+			throw new IllegalArgumentException("Radius must be > 0 but was " + radius);
+		}
+		double[] alignedCenter = MathTools.floorPowerOfTwoCopy(center);
+		double alignedRadius = MathTools.ceilPowerOfTwo(radius);
+		t.root = new QNode<>(Arrays.copyOf(alignedCenter, alignedCenter.length), alignedRadius);
+		return t;
+	}
+
+	/**
+	 * WARNING: Unaligned center and radius can cause precision problems.
+	 * @param dims dimensions, usually 2 or 3
+	 * @param maxNodeSize maximum entries per node, default is 10
+	 * @param center center of initial root node
+	 * @param radius radius of initial root node
+	 * @return New quadtree
+	 * @param <T> Value type
+	 * @deprecated Please use {@link #createAligned(int, int, double[], double)}
+	 */
+	@Deprecated
 	public static <T> QuadTreeKD2<T> create(int dims, int maxNodeSize, 
 			double[] center, double radius) {
 		QuadTreeKD2<T> t = new QuadTreeKD2<>(dims, maxNodeSize);