1 package org.davidmoten.hilbert;
2
3 import java.math.BigInteger;
4 import java.util.Arrays;
5
6 import com.github.davidmoten.guavamini.Preconditions;
7 import com.github.davidmoten.guavamini.annotations.VisibleForTesting;
8
9 /**
10 * Converts between Hilbert index ({@code BigInteger}) and N-dimensional points.
11 *
12 * <p>
13 * Note: This algorithm is derived from work done by John Skilling and published
14 * in "Programming the Hilbert curve". (c) 2004 American Institute of Physics.
15 * With thanks also to Paul Chernoch who published a C# algorithm for Skilling's
16 * work on StackOverflow and
17 * <a href="https://github.com/paulchernoch/HilbertTransformation">GitHub</a>).
18 */
19 public final class HilbertCurve {
20
21 private final int bits;
22 private final int dimensions;
23 // cached calculations
24 private final int length;
25
26 private HilbertCurve(int bits, int dimensions) {
27 this.bits = bits;
28 this.dimensions = dimensions;
29 // cache a calculated values for small perf improvements
30 this.length = bits * dimensions;
31 }
32
33 /**
34 * Returns a builder for and object that performs transformations for a
35 * Hilbert curve with the given number of bits.
36 *
37 * @param bits
38 * depth of the Hilbert curve. If bits is one, this is the
39 * top-level Hilbert curve
40 * @return builder for object to do transformations with the Hilbert Curve
41 */
42 public static Builder bits(int bits) {
43 return new Builder(bits);
44 }
45
46 public static SmallHilbertCurve.Builder small() {
47 return new SmallHilbertCurve.Builder();
48 }
49
50 /**
51 * Builds a {@link HilbertCurve} instance.
52 */
53 public static final class Builder {
54 final int bits;
55
56 private Builder(int bits) {
57 Preconditions.checkArgument(bits > 0, "bits must be greater than zero");
58 Preconditions.checkArgument(bits < 64, "bits must be 63 or less");
59 this.bits = bits;
60 }
61
62 public HilbertCurve dimensions(int dimensions) {
63 Preconditions.checkArgument(dimensions > 1, "dimensions must be at least 2");
64 return new HilbertCurve(bits, dimensions);
65 }
66 }
67
68 /**
69 * Converts a point to its Hilbert curve index.
70 *
71 * @param point
72 * an array of {@code long}. Each coordinate can be between 0 and
73 * 2<sup>bits</sup>-1.
74 * @return index (nonnegative {@link BigInteger})
75 * @throws IllegalArgumentException
76 * if length of point array is not equal to the number of
77 * dimensions.
78 */
79 public BigInteger index(long... point) {
80 Preconditions.checkArgument(point.length == dimensions);
81 return toIndex(transposedIndex(bits, point));
82 }
83
84 /**
85 * Converts a {@link BigInteger} index (distance along the Hilbert Curve
86 * from 0) to a point of dimensions defined in the constructor of
87 * {@code this}.
88 *
89 * @param index
90 * index along the Hilbert Curve from 0. Maximum value 2
91 * <sup>bits * dimensions</sup>-1.
92 * @return array of longs being the point
93 * @throws NullPointerException
94 * if index is null
95 * @throws IllegalArgumentException
96 * if index is negative
97 */
98 public long[] point(BigInteger index) {
99 Preconditions.checkNotNull(index);
100 Preconditions.checkArgument(index.signum() != -1, "index cannot be negative");
101 return transposedIndexToPoint(bits, transpose(index));
102 }
103
104 /**
105 * Converts a {@code long} index (distance along the Hilbert Curve from 0)
106 * to a point of dimensions defined in the constructor of {@code this}.
107 *
108 * @param index
109 * index along the Hilbert Curve from 0. Maximum value 2
110 * <sup>bits+1</sup>-1.
111 * @return array of longs being the point
112 * @throws IllegalArgumentException
113 * if index is negative
114 */
115 public long[] point(long index) {
116 return point(BigInteger.valueOf(index));
117 }
118
119 /**
120 * Returns the transposed representation of the Hilbert curve index.
121 *
122 * <p>
123 * The Hilbert index is expressed internally as an array of transposed bits.
124 *
125 * <pre>
126 Example: 5 bits for each of n=3 coordinates.
127 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
128 as its Transpose ^
129 X[0] = A D G J M X[2]| 7
130 X[1] = B E H K N <-------> | /X[1]
131 X[2] = C F I L O axes |/
132 high low 0------> X[0]
133 * </pre>
134 *
135 * @param index
136 * index to be tranposed
137 * @return transposed index
138 */
139 @VisibleForTesting
140 long[] transpose(BigInteger index) {
141 byte[] b = index.toByteArray();
142 long[] x = new long[dimensions];
143 for (int idx = 0; idx < 8 * b.length; idx++) {
144 if ((b[b.length - 1 - idx / 8] & (1L << (idx % 8))) != 0) {
145 int dim = (length - idx - 1) % dimensions;
146 int shift = (idx / dimensions) % bits;
147 x[dim] |= 1L << shift;
148 }
149 }
150 return x;
151 }
152
153 /**
154 * <p>
155 * Given the axes (coordinates) of a point in N-Dimensional space, find the
156 * distance to that point along the Hilbert curve. That distance will be
157 * transposed; broken into pieces and distributed into an array.
158 *
159 * <p>
160 * The number of dimensions is the length of the hilbertAxes array.
161 *
162 * <p>
163 * Note: In Skilling's paper, this function is called AxestoTranspose.
164 *
165 * @param mutate
166 *
167 * @param point
168 * Point in N-space
169 * @return The Hilbert distance (or index) as a transposed Hilbert index
170 */
171 @VisibleForTesting
172 static long[] transposedIndex(int bits, long... point) {
173 final long M = 1L << (bits - 1);
174 final int n = point.length; // n: Number of dimensions
175 final long[] x = Arrays.copyOf(point, n);
176 long p, q, t;
177 int i;
178 // Inverse undo
179 for (q = M; q > 1; q >>= 1) {
180 p = q - 1;
181 for (i = 0; i < n; i++)
182 if ((x[i] & q) != 0)
183 x[0] ^= p; // invert
184 else {
185 t = (x[0] ^ x[i]) & p;
186 x[0] ^= t;
187 x[i] ^= t;
188 }
189 } // exchange
190 // Gray encode
191 for (i = 1; i < n; i++)
192 x[i] ^= x[i - 1];
193 t = 0;
194 for (q = M; q > 1; q >>= 1)
195 if ((x[n - 1] & q) != 0)
196 t ^= q - 1;
197 for (i = 0; i < n; i++)
198 x[i] ^= t;
199
200 return x;
201 }
202
203 /**
204 * Converts the Hilbert transposed index into an N-dimensional point
205 * expressed as a vector of {@code long}.
206 *
207 * In Skilling's paper this function is named {@code TransposeToAxes}
208 *
209 * @param transposedIndex
210 * distance along the Hilbert curve in transposed form
211 * @return the coordinates of the point represented by the transposed index
212 * on the Hilbert curve
213 */
214 static long[] transposedIndexToPoint(int bits, long... x) {
215 final long N = 2L << (bits - 1);
216 // Note that x is mutated by this method (as a performance improvement
217 // to avoid allocation)
218 int n = x.length; // number of dimensions
219 long p, q, t;
220 int i;
221 // Gray decode by H ^ (H/2)
222 t = x[n - 1] >> 1;
223 // Corrected error in Skilling's paper on the following line. The
224 // appendix had i >= 0 leading to negative array index.
225 for (i = n - 1; i > 0; i--)
226 x[i] ^= x[i - 1];
227 x[0] ^= t;
228 // Undo excess work
229 for (q = 2; q != N; q <<= 1) {
230 p = q - 1;
231 for (i = n - 1; i >= 0; i--)
232 if ((x[i] & q) != 0L)
233 x[0] ^= p; // invert
234 else {
235 t = (x[0] ^ x[i]) & p;
236 x[0] ^= t;
237 x[i] ^= t;
238 }
239 } // exchange
240 return x;
241 }
242
243 // Quote from Paul Chernoch
244 // Interleaving means take one bit from the first matrix element, one bit
245 // from the next, etc, then take the second bit from the first matrix
246 // element, second bit from the second, all the way to the last bit of the
247 // last element. Combine those bits in that order into a single BigInteger,
248 // which can have as many bits as necessary. This converts the array into a
249 // single number.
250 @VisibleForTesting
251 BigInteger toIndex(long... transposedIndex) {
252 byte[] b = new byte[length];
253 int bIndex = length - 1;
254 long mask = 1L << (bits - 1);
255 for (int i = 0; i < bits; i++) {
256 for (int j = 0; j < transposedIndex.length; j++) {
257 if ((transposedIndex[j] & mask) != 0) {
258 b[length - 1 - bIndex / 8] |= 1 << (bIndex % 8);
259 }
260 bIndex--;
261 }
262 mask >>= 1;
263 }
264 // b is expected to be BigEndian
265 return new BigInteger(1, b);
266 }
267
268 }