1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
/*********************************************************************
* kd-forest *
* Copyright (C) 2015 Tavian Barnes <tavianator@tavianator.com> *
* *
* This program is free software. It comes without any warranty, to *
* the extent permitted by applicable law. You can redistribute it *
* and/or modify it under the terms of the Do What The Fuck You Want *
* To Public License, Version 2, as published by Sam Hocevar. See *
* the COPYING file or http://www.wtfpl.net/ for more details. *
*********************************************************************/
#include "hilbert.h"
#include <stdint.h>
// These algorithms are described in "Compact Hilbert Indices" by Chris Hamilton
// Right rotation of x by b bits out of n
static uint32_t
ror(uint32_t x, unsigned int b, unsigned int n)
{
uint32_t l = x & ((1 << b) - 1);
uint32_t r = x >> b;
return (l << (n - b)) | r;
}
// Left rotation of x by b bits out of n
static uint32_t
rol(uint32_t x, unsigned int b, unsigned int n)
{
return ror(x, n - b, n);
}
// Binary reflected Gray code
uint32_t
gray_code(uint32_t i)
{
return i ^ (i >> 1);
}
// e(i), the entry point for the ith sub-hypercube
uint32_t
entry_point(uint32_t i)
{
if (i == 0) {
return 0;
} else {
return gray_code((i - 1) & ~1U);
}
}
// g(i), the inter sub-hypercube direction
unsigned int
inter_direction(uint32_t i)
{
// g(i) counts the trailing set bits in i
unsigned int g = 0;
while (i & 1) {
++g;
i >>= 1;
}
return g;
}
// d(i), the intra sub-hypercube direction
unsigned int
intra_direction(uint32_t i)
{
if (i & 1) {
return inter_direction(i);
} else if (i > 0) {
return inter_direction(i - 1);
} else {
return 0;
}
}
// T transformation inverse
uint32_t
t_inverse(unsigned int dimensions, uint32_t e, unsigned int d, uint32_t a)
{
return rol(a, d, dimensions) ^ e;
}
// GrayCodeRankInverse
void
gray_code_rank_inverse(unsigned int dimensions, uint32_t mu, uint32_t pi, unsigned int r, unsigned int free_bits, uint32_t *i, uint32_t *g)
{
// *i is the inverse rank of r
// *g is gray_code(i)
*i = 0;
*g = 0;
for (unsigned int j = free_bits - 1, k = dimensions; k-- > 0;) {
if (mu & (1 << k)) {
*i |= ((r >> j) & 1) << k;
*g |= (*i ^ (*i >> 1)) & (1 << k);
--j;
} else {
*g |= pi & (1 << k);
*i |= (*g ^ (*i >> 1)) & (1 << k);
}
}
}
// ExtractMask
void
extract_mask(unsigned int dimensions, const unsigned int extents[], unsigned int i, uint32_t *mu, unsigned int *free_bits)
{
// *mu is the mask
// *free_bits is popcount(*mu)
*mu = 0;
*free_bits = 0;
for (unsigned int j = dimensions; j-- > 0;) {
*mu <<= 1;
if (extents[j] > i) {
*mu |= 1;
*free_bits += 1;
}
}
}
// CompactHilbertIndexInverse
void
hilbert_point(unsigned int dimensions, const unsigned int extents[], uint32_t index, uint32_t point[])
{
unsigned int max_extent = 0, total_extent = 0;
for (unsigned int i = 0; i < dimensions; ++i) {
if (extents[i] > max_extent) {
max_extent = extents[i];
}
total_extent += extents[i];
point[i] = 0;
}
uint32_t e = 0;
unsigned int k = 0;
// Next direction; we use d instead of d + 1 everywhere
unsigned int d = 1;
for (unsigned int i = max_extent; i-- > 0;) {
uint32_t mu;
unsigned int free_bits;
extract_mask(dimensions, extents, i, &mu, &free_bits);
mu = ror(mu, d, dimensions);
uint32_t pi = ror(e, d, dimensions) & ~mu;
unsigned int r = (index >> (total_extent - k - free_bits)) & ((1 << free_bits) - 1);
k += free_bits;
uint32_t w, l;
gray_code_rank_inverse(dimensions, mu, pi, r, free_bits, &w, &l);
l = t_inverse(dimensions, e, d, l);
for (unsigned int j = 0; j < 3; ++j) {
point[j] |= (l & 1) << i;
l >>= 1;
}
e = e ^ ror(entry_point(w), d, dimensions);
d = (d + intra_direction(w) + 1)%dimensions;
}
}
|
