Class HilbertCode

1	/*
2	* Copyright (c) 2019 Martin Davis.
3	*
4	* All rights reserved. This program and the accompanying materials
5	* are made available under the terms of the Eclipse Public License 2.0
6	* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
7	* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
8	* and the Eclipse Distribution License is available at
9	*
10	* http://www.eclipse.org/org/documents/edl-v10.php.
11	*/
12
13	package org.locationtech.jts.shape.fractal;
14
15	import org.locationtech.jts.geom.Coordinate;
16	import org.locationtech.jts.geom.Geometry;
17	import org.locationtech.jts.geom.GeometryFactory;
18	import org.locationtech.jts.geom.LineSegment;
19	import org.locationtech.jts.shape.GeometricShapeBuilder;
20
21	/**
22	* Encodes points as the index along finite planar Hilbert curves.
23	* <p>
24	* The planar Hilbert Curve is a continuous space-filling curve.
25	* In the limit the Hilbert curve has infinitely many vertices and fills
26	* the space of the unit square.
27	* A sequence of finite approximations to the infinite Hilbert curve
28	* is defined by the level number.
29	* The finite Hilbert curve at level n H<sub>n</sub> contains 2<sup>n + 1</sup> points.
30	* Each finite Hilbert curve defines an ordering of the
31	* points in the 2-dimensional range square containing the curve.
32	* Curves fills the range square of side 2<sup>level</sup>.
33	* Curve points have ordinates in the range [0, 2<sup>level</sup> - 1].
34	* The index of a point along a Hilbert curve is called the Hilbert code.
35	* The code for a given point is specific to the level chosen.
36	* <p>
37	* This implementation represents codes using 32-bit integers.
38	* This allows levels 0 to 16 to be handled.
39	* The class supports encoding points in the range of a given level curve
40	* and decoding the point for a given code value.
41	* <p>
42	* The Hilbert order has the property that it tends to preserve locality.
43	* This means that codes which are near in value will have spatially proximate
44	* points. The converse is not always true - the delta between
45	* codes for nearby points is not always small. But the average delta
46	* is small enough that the Hilbert order is an effective way of linearizing space
47	* to support range queries.
48	*
49	* @author Martin Davis
50	*
51	* @see HilbertCurveBuilder
52	* @see MortonCode
53	*/
54	public class HilbertCode
55	{
56	/**
57	* The maximum curve level that can be represented.
58	*/
59	public static final int MAX_LEVEL = 16;
60
61	/**
62	* The number of points in the curve for the given level.
63	* The number of points is 2<sup>2 * level</sup>.
64	*
65	* @param level the level of the curve
66	* @return the number of points
67	*/
68	public static int size(int level) {
69	checkLevel(level);
70	return (int) Math.pow(2, 2 *level);
71	}
72
73	/**
74	* The maximum ordinate value for points
75	* in the curve for the given level.
76	* The maximum ordinate is 2<sup>level</sup> - 1.
77	*
78	* @param level the level of the curve
79	* @return the maximum ordinate value
80	*/
81	public static int maxOrdinate(int level) {
82	checkLevel(level);
83	return (int) Math.pow(2, level) - 1;
84	}
85
86	/**
87	* The level of the finite Hilbert curve which contains at least
88	* the given number of points.
89	*
90	* @param numPoints the number of points required
91	* @return the level of the curve
92	*/
93	public static int level(int numPoints) {
94	int pow2 = (int) ( (Math.log(numPoints)/Math.log(2)));
95	int level = pow2 / 2;
96	int size = size(level);
97	if (size < numPoints) level += 1;
98	return level;
99	}
100
101	private static void checkLevel(int level) {
102	if (level > MAX_LEVEL) {
103	throw new IllegalArgumentException("Level must be in range 0 to " + MAX_LEVEL);
104	}
105	}
106
107	/**
108	* Encodes a point (x,y)
109	* in the range of the the Hilbert curve at a given level
110	* as the index of the point along the curve.
111	* The index will lie in the range [0, 2<sup>level + 1</sup>].
112	*
113	* @param level the level of the Hilbert curve
114	* @param x the x ordinate of the point
115	* @param y the y ordinate of the point
116	* @return the index of the point along the Hilbert curve
117	*/
118	public static int encode(int level, int x, int y) {
119	// Fast Hilbert curve algorithm by http://threadlocalmutex.com/
120	// Ported from C++ https://github.com/rawrunprotected/hilbert_curves (public
121	// domain)
122
123	int lvl = levelClamp(level);
124
125	x = x << (16 - lvl);
126	y = y << (16 - lvl);
127
128	long a = x ^ y;
129	long b = 0xFFFF ^ a;
130	long c = 0xFFFF ^ (x \| y);
131	long d = x & (y ^ 0xFFFF);
132
133	long A = a \| (b >> 1);
134	long B = (a >> 1) ^ a;
135	long C = ((c >> 1) ^ (b & (d >> 1))) ^ c;
136	long D = ((a & (c >> 1)) ^ (d >> 1)) ^ d;
137
138	a = A;
139	b = B;
140	c = C;
141	d = D;
142	A = ((a & (a >> 2)) ^ (b & (b >> 2)));
143	B = ((a & (b >> 2)) ^ (b & ((a ^ b) >> 2)));
144	C ^= ((a & (c >> 2)) ^ (b & (d >> 2)));
145	D ^= ((b & (c >> 2)) ^ ((a ^ b) & (d >> 2)));
146
147	a = A;
148	b = B;
149	c = C;
150	d = D;
151	A = ((a & (a >> 4)) ^ (b & (b >> 4)));
152	B = ((a & (b >> 4)) ^ (b & ((a ^ b) >> 4)));
153	C ^= ((a & (c >> 4)) ^ (b & (d >> 4)));
154	D ^= ((b & (c >> 4)) ^ ((a ^ b) & (d >> 4)));
155
156	a = A;
157	b = B;
158	c = C;
159	d = D;
160	C ^= ((a & (c >> 8)) ^ (b & (d >> 8)));
161	D ^= ((b & (c >> 8)) ^ ((a ^ b) & (d >> 8)));
162
163	a = C ^ (C >> 1);
164	b = D ^ (D >> 1);
165
166	long i0 = x ^ y;
167	long i1 = b \| (0xFFFF ^ (i0 \| a));
168
169	i0 = (i0 \| (i0 << 8)) & 0x00FF00FF;
170	i0 = (i0 \| (i0 << 4)) & 0x0F0F0F0F;
171	i0 = (i0 \| (i0 << 2)) & 0x33333333;
172	i0 = (i0 \| (i0 << 1)) & 0x55555555;
173
174	i1 = (i1 \| (i1 << 8)) & 0x00FF00FF;
175	i1 = (i1 \| (i1 << 4)) & 0x0F0F0F0F;
176	i1 = (i1 \| (i1 << 2)) & 0x33333333;
177	i1 = (i1 \| (i1 << 1)) & 0x55555555;
178
179	long index = ((i1 << 1) \| i0) >> (32 - 2 * lvl);
180	return (int) index;
181	}
182
183	/**
184	* Clamps a level to the range valid for
185	* the index algorithm used.
186	*
187	* @param level the level of a Hilbert curve
188	* @return a valid level
189	*/
190	private static int levelClamp(int level) {
191	// clamp order to [1, 16]
192	int lvl = level < 1 ? 1 : level;
193	lvl = lvl > MAX_LEVEL ? MAX_LEVEL : lvl;
194	return lvl;
195	}
196
197	/**
198	* Computes the point on a Hilbert curve
199	* of given level for a given code index.
200	* The point ordinates will lie in the range [0, 2<sup>level</sup></i> - 1].
201	*
202	* @param level the Hilbert curve level
203	* @param index the index of the point on the curve
204	* @return the point on the Hilbert curve
205	*/
206	public static Coordinate decode(int level, int index) {
207	checkLevel(level);
208	int lvl = levelClamp(level);
209
210	index = index << (32 - 2 * lvl);
211
212	long i0 = deinterleave(index);
213	long i1 = deinterleave(index >> 1);
214
215	long t0 = (i0 \| i1) ^ 0xFFFF;
216	long t1 = i0 & i1;
217
218	long prefixT0 = prefixScan(t0);
219	long prefixT1 = prefixScan(t1);
220
221	long a = (((i0 ^ 0xFFFF) & prefixT1) \| (i0 & prefixT0));
222
223	long x = (a ^ i1) >> (16 - lvl);
224	long y = (a ^ i0 ^ i1) >> (16 - lvl);
225
226	return new Coordinate(x, y);
227	}
228
229	private static long prefixScan(long x) {
230	x = (x >> 8) ^ x;
231	x = (x >> 4) ^ x;
232	x = (x >> 2) ^ x;
233	x = (x >> 1) ^ x;
234	return x;
235	}
236
237	private static long deinterleave(int x) {
238	x = x & 0x55555555;
239	x = (x \| (x >> 1)) & 0x33333333;
240	x = (x \| (x >> 2)) & 0x0F0F0F0F;
241	x = (x \| (x >> 4)) & 0x00FF00FF;
242	x = (x \| (x >> 8)) & 0x0000FFFF;
243	return x;
244	}
245	}
246