Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
185 lines (159 sloc) 4.57 KB
/*
Copyright 2019 jtejido
Package hilbert is created to support multidimensional encoding/decoding to/from a Space-Filling Curve
http://en.wikipedia.org/wiki/Hilbert_curve
*/
package hilbert
import (
"errors"
"math/big"
)
var (
ErrNNotPositive = errors.New("Dimension must be greater than zero")
ErrBNotPositive = errors.New("Number of bits must be greater than zero")
)
// This algorithm is derived from work done by John Skilling and published
// in "Programming the Hilbert curve". (c) 2004 American Institute of Physics.
// https://doi.org/10.1063/1.1751381
type Hilbert struct {
bits, dimension, length uint32
}
func New(b, n uint32) (*Hilbert, error) {
if n <= 0 {
return nil, ErrNNotPositive
}
if b <= 0 {
return nil, ErrBNotPositive
}
return &Hilbert{
bits: b,
dimension: n,
length: b * n,
}, nil
}
// The number of dimensions
func (s *Hilbert) Dimension() uint32 {
return s.dimension
}
// The number of bits
func (s *Hilbert) Bits() uint32 {
return s.bits
}
// The number of length
func (s *Hilbert) Len() uint32 {
return s.length
}
// Converts points to its Hilbert curve index.
func (s *Hilbert) Encode(x ...uint64) *big.Int {
return s.untranspose(s.axesToTranspose(x...))
}
// Converts an index (distance along the Hilbert Curve from 0)
// to a point of dimensions defined
func (s *Hilbert) Decode(index *big.Int) []uint64 {
return s.transposedToAxes(s.transpose(index))
}
// The high-order bit from the last number in vector becomes the high-order bit of last byte in the generated byte array.
// The high-order bit of the next to last number becomes the second highest-ordered bit in the last byte in the generated byte array.
// The low-order bit of the first number becomes the low order bit of the first byte in the new array.
func (s *Hilbert) untranspose(x []uint64) *big.Int {
t := make([]byte, s.length)
bIndex := s.length - 1
mask := uint64(1 << (s.bits - 1))
for i := 0; i < int(s.bits); i++ {
for j := 0; j < len(x); j++ {
if (x[j] & mask) != 0 {
t[s.length-1-bIndex/8] |= 1 << (bIndex % 8)
//b |= 1 << (bIndex % 8) // small int
}
bIndex--
}
mask >>= 1
}
return new(big.Int).SetBytes(t)
}
// Returns the transposed representation of the Hilbert curve index.
// The Hilbert index is expressed internally as an array of transposed bits.
// Example: 5 bits for each of n=3 coordinates.
// 15-bit Hilbert integer = A B C D E F G H I J K L M N O is stored
// as its Transpose ^
// X[0] = A D G J M X[2]| 7
// X[1] = B E H K N <-------> | /X[1]
// X[2] = C F I L O axes |/
// high low 0------> X[0]
func (s *Hilbert) transpose(index *big.Int) []uint64 {
x := make([]uint64, s.dimension)
b := index.Bytes()
for idx := 0; idx < 8*len(b); idx++ {
if (b[len(b)-1-idx/8] & (1 << (uint32(idx) % 8))) != 0 {
dim := (s.length - uint32(idx) - 1) % s.dimension
shift := (uint32(idx) / s.dimension) % s.bits
x[dim] |= 1 << shift
}
}
return x
}
// Converts the Hilbert transposed index into an N-dimensional point expressed
// as a vector of int.
func (s *Hilbert) transposedToAxes(x []uint64) []uint64 {
N := uint64(2 << (s.bits - 1))
// Note that x is mutated by this method (as a performance improvement
// to avoid allocation)
n := len(x)
// Gray decode by H ^ (H/2)
t := x[n-1] >> 1
// Corrected error in Skilling's paper on the following line. The
// appendix had i >= 0 leading to negative array index.
for i := n - 1; i > 0; i-- {
x[i] ^= x[i-1]
}
x[0] ^= t
// Undo excess work
for q := uint64(2); q != N; q <<= 1 {
p := q - 1
for i := n - 1; i >= 0; i-- {
if (x[i] & q) != 0 {
x[0] ^= p // invert
} else {
t = (x[0] ^ x[i]) & p
x[0] ^= t
x[i] ^= t
}
}
} // exchange
return x
}
// Given the axes (coordinates) of a point in N-Dimensional space, find the
// distance to that point along the Hilbert curve. That distance will be
// transposed; broken into pieces and distributed into an array.
// The number of dimensions is the length of the hilbertAxes array.
func (s *Hilbert) axesToTranspose(x ...uint64) []uint64 {
M := uint64(1 << (s.bits - 1))
n := len(x)
var t uint64
for q := M; q > 1; q >>= 1 {
p := q - 1
for i := 0; i < n; i++ {
if (x[i] & q) != 0 {
x[0] ^= p // invert
} else {
t = (x[0] ^ x[i]) & p
x[0] ^= t
x[i] ^= t
}
}
}
// Gray encode
for i := 1; i < n; i++ {
x[i] ^= x[i-1]
}
t = 0
for q := M; q > 1; q >>= 1 {
if (x[n-1] & q) != 0 {
t ^= q - 1
}
}
for i := 0; i < n; i++ {
x[i] ^= t
}
return x
}
You can’t perform that action at this time.