org.davidmoten.hilbert.SmallHilbertCurve Maven / Gradle / Ivy
package org.davidmoten.hilbert;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import com.github.davidmoten.guavamini.Preconditions;
/**
* Converts between Hilbert index ({@code BigInteger}) and N-dimensional points.
*
*
* Note: This algorithm is derived from work done by John Skilling and published
* in "Programming the Hilbert curve". (c) 2004 American Institute of Physics.
* With thanks also to Paul Chernoch who published a C# algorithm for Skilling's
* work on StackOverflow and
* GitHub).
*/
public final class SmallHilbertCurve {
private static final int DEFAULT_BUFFER_SIZE = 1024;
private final int bits;
private final int dimensions;
private final int length;
private SmallHilbertCurve(int bits, int dimensions) {
this.bits = bits;
this.dimensions = dimensions;
this.length = bits * dimensions;
}
public int bits() {
return bits;
}
public int dimensions() {
return dimensions;
}
/**
* Converts a point to its Hilbert curve index.
*
* @param point an array of {@code long}. Each coordinate can be between 0 and
* 2bits-1.
* @return index {@code long} in the range 0 to 2bits * dimensions -
* 1
* @throws IllegalArgumentException if length of point array is not equal to the
* number of dimensions.
*/
public long index(long... point) {
Preconditions.checkArgument(point.length == dimensions);
return toIndex(HilbertCurve.transposedIndex(bits, point));
}
/**
* Converts a {@code long} index (distance along the Hilbert Curve from 0) to a
* point of dimensions defined in the constructor of {@code this}.
*
* @param index index along the Hilbert Curve from 0. Maximum value 2 bits
* * dimensions-1.
* @return array of longs being the point
* @throws IllegalArgumentException if index is negative
*/
public long[] point(long index) {
return HilbertCurve.transposedIndexToPoint(bits, transposeLong(index));
}
public void point(long index, long[] x) {
Util.zero(x);
transposeLong(index, x);
HilbertCurve.transposedIndexToPoint(bits, x);
}
// untranspose
private long toIndex(long... transposedIndex) {
long b = 0;
int bIndex = length - 1;
long mask = 1L << (bits - 1);
for (int i = 0; i < bits; i++) {
for (int j = 0; j < transposedIndex.length; j++) {
if ((transposedIndex[j] & mask) != 0) {
b |= 1L << bIndex;
}
bIndex--;
}
mask >>= 1;
}
// b is expected to be BigEndian
return b;
}
private void transposeLong(long index, long[] x) {
for (int idx = 0; idx < 64; idx++) {
if ((index & (1L << idx)) != 0) {
int dim = (length - idx - 1) % dimensions;
int shift = (idx / dimensions) % bits;
x[dim] |= 1L << shift;
}
}
}
private long[] transposeLong(long index) {
long[] x = new long[dimensions];
transposeLong(index, x);
return x;
}
public long maxOrdinate() {
return (1 << bits) - 1;
}
public long maxIndex() {
return (1 << (bits * dimensions)) - 1;
}
/////////////////////////////////////////////////
// Query support
////////////////////////////////////////////////
/**
* Returns index ranges exactly covering the region bounded by {@code a} and
* {@code b}. The list will be in increasing order of the range bounds (there
* should be no overlaps).
*
* @param a one vertex of the region
* @param b the opposing vertex to a
* @return ranges
*/
public Ranges query(long[] a, long[] b) {
return query(a, b, 0, 0);
}
/**
* Returns index ranges covering the region bounded by {@code a} and {@code b}.
* The list will be in increasing order of the range bounds (there should be no
* overlaps). The index ranges may cover a larger region than the search box
* because the set of exact covering ranges will have been reduced by joining
* ranges with minimal gaps. The buffer size used by this method is 1024.
*
* @param a one vertex of the region
* @param b the opposing vertex to a
* @param maxRanges the maximum number of ranges to be returned. If 0 then all
* ranges are returned.
* @return ranges
*/
public Ranges query(long[] a, long[] b, int maxRanges) {
if (maxRanges == 0) {
return query(a, b, 0, 0);
} else {
return query(a, b, maxRanges, Math.max(DEFAULT_BUFFER_SIZE, maxRanges));
}
}
/**
* Returns index ranges covering the region bounded by {@code a} and {@code b}.
* The list will be in increasing order of the range bounds (there should be no
* overlaps). The index ranges may cover a larger region than the search box
* because the set of exact covering ranges will have been reduced by joining
* ranges with minimal gaps. A buffer of ranges is used and has a size that is
* generally speaking a fair bit larger than maxRanges to increase the
* probability that minimal extra coverage be returned.
*
* @param a one vertex of the region
* @param b the opposing vertex to a
* @param maxRanges the maximum number of ranges to be returned. If 0 then all
* ranges are returned.
* @param bufferSize the buffer size of ranges to use. A larger buffer size will
* increase the probability that ranges have been joined
* optimally (we want to join ranges that have minimal gaps to
* minimize the resultant extra coverage). If 0 is passed to
* bufferSize then all ranges will be buffered before
* shrinking to maxRanges.
* @return ranges
*/
public Ranges query(long[] a, long[] b, int maxRanges, int bufferSize) {
Preconditions.checkArgument(maxRanges >= 0);
Preconditions.checkArgument(bufferSize >= maxRanges,
"bufferSize must be greater than or equal to maxRanges");
if (maxRanges == 0) {
// unlimited
bufferSize = 0;
}
// this is the implementation of the Perimeter Algorithm mentioned in README.md
Box box = new Box(a, b);
List list = new ArrayList<>();
box.visitPerimeter(cell -> {
long n = index(cell);
list.add(n);
});
Collections.sort(list);
int i = 0;
Ranges ranges = new Ranges(bufferSize);
long rangeStart = -1;
while (true) {
if (i == list.size()) {
break;
}
if (rangeStart == -1) {
rangeStart = list.get(i);
}
while (i < list.size() - 1 && list.get(i + 1) == list.get(i) + 1) {
i++;
}
if (i == list.size() - 1) {
ranges.add(Range.create(rangeStart, list.get(i)));
break;
}
long[] point = point(list.get(i) + 1);
if (box.contains(point)) {
// is not on the perimeter (would have been caught in previous while loop)
// so is internal to the box which means the next value in the sorted hilbert
// curve indexes for the perimiter must be where it exits
i += 1;
} else {
ranges.add(Range.create(rangeStart, list.get(i)));
rangeStart = -1;
i++;
}
}
if (ranges.size() <= maxRanges) {
return ranges;
} else {
Ranges r = new Ranges(maxRanges);
for (Range range : ranges) {
r.add(range);
}
return r;
}
}
public static final class Builder {
private int bits;
Builder() {
// private instantiation
}
public Builder bits(int bits) {
this.bits = bits;
return this;
}
public SmallHilbertCurve dimensions(int dimensions) {
Preconditions.checkArgument(bits * dimensions <= 63,
"bits * dimensions must be less than or equal to 63");
return new SmallHilbertCurve(bits, dimensions);
}
}
}