rt300@34
|
1 //
|
rt300@34
|
2 // hilbert.cpp
|
rt300@34
|
3 // sonicZoom
|
rt300@34
|
4 //
|
rt300@34
|
5 // Created by Robert Tubb on 04/04/2013.
|
rt300@34
|
6 //
|
rt300@34
|
7 //
|
rt300@34
|
8
|
rt300@34
|
9 #include "hilbert.h"
|
rt300@34
|
10 Hilbert hilbert;
|
rt300@34
|
11
|
rt300@34
|
12 //--------------------------------------------------------------------
|
rt300@34
|
13 unsigned int rotl(unsigned int value, int shift, int digits) {
|
rt300@34
|
14 unsigned int ones = (1 << digits) - 1;
|
rt300@34
|
15 if(shift >= digits)
|
rt300@34
|
16 shift -= digits;
|
rt300@34
|
17
|
rt300@34
|
18 unsigned int out = (value << shift) | (value >> (digits - shift));
|
rt300@34
|
19 return out & ones;
|
rt300@34
|
20 }
|
rt300@34
|
21 //--------------------------------------------------------------------
|
rt300@34
|
22 unsigned int rotr(unsigned int value, int shift, int digits) {
|
rt300@34
|
23 int ones = (1 << digits) - 1;
|
rt300@34
|
24 // can only handle -digits to 2*digits-1
|
rt300@34
|
25 if(shift >= digits) shift = shift - digits;
|
rt300@34
|
26 if(shift < 0) shift = digits + shift;
|
rt300@34
|
27 unsigned int out = (value >> shift) | (value << (digits - shift));
|
rt300@34
|
28 return out & ones;
|
rt300@34
|
29 }
|
rt300@34
|
30 //--------------------------------------------------------------------
|
rt300@34
|
31 unsigned int bin2dec(vector<bool> binary){
|
rt300@34
|
32 // takes in a vector of booleans. we assume that the leftmost bool (highest index) is the MSB
|
rt300@34
|
33 unsigned int dec = 0;
|
rt300@34
|
34 int i,pwr;
|
rt300@34
|
35 int N = binary.size();
|
rt300@34
|
36
|
rt300@34
|
37 for (i = 0; i<N; i++){
|
rt300@34
|
38 pwr = N-i-1;
|
rt300@34
|
39 dec += binary[i]*(1 << pwr);
|
rt300@34
|
40
|
rt300@34
|
41 }
|
rt300@34
|
42
|
rt300@34
|
43 return dec;
|
rt300@34
|
44 }
|
rt300@34
|
45 //--------------------------------------------------------------------
|
rt300@34
|
46 vector<bool> dec2bin(unsigned int number, unsigned int size){
|
rt300@34
|
47 vector<bool> bin;
|
rt300@34
|
48 int i, pwr;
|
rt300@34
|
49 for(i=0;i<size;i++){
|
rt300@34
|
50 bin.push_back(0);
|
rt300@34
|
51 }
|
rt300@34
|
52 if(number >= (1 << size)){
|
rt300@34
|
53 cout << " bad size for Hilbert::dec2bin" << "\n";
|
rt300@34
|
54 return bin;
|
rt300@34
|
55 }
|
rt300@34
|
56 for(i=0;i<size;i++){
|
rt300@34
|
57 pwr = size-i;
|
rt300@34
|
58 bin[i] = (number/(1 << pwr)) % 2;
|
rt300@34
|
59 bin.push_back(0);
|
rt300@34
|
60 }
|
rt300@35
|
61 return bin;
|
rt300@34
|
62
|
rt300@34
|
63
|
rt300@34
|
64 }
|
rt300@35
|
65 //--------------------------------------------------------------------
|
rt300@34
|
66 void printBool(vector <bool> vcode){
|
rt300@34
|
67 vector<bool>::iterator bit;
|
rt300@34
|
68
|
rt300@34
|
69 for(bit=vcode.begin(); bit!=vcode.end() ; bit++){
|
rt300@34
|
70
|
rt300@34
|
71 cout << *bit;
|
rt300@34
|
72 }
|
rt300@34
|
73 cout << "\n";
|
rt300@34
|
74
|
rt300@34
|
75 }
|
rt300@34
|
76 //--------------------------------------------------------------------
|
rt300@34
|
77 void Hilbert::init(int aN, int aP){
|
rt300@34
|
78 N = aN;
|
rt300@34
|
79 P = aP;
|
rt300@34
|
80 codeLength = (int)pow(2.0,P);
|
rt300@34
|
81 makeCode();
|
rt300@34
|
82 makeRotationRules();
|
rt300@34
|
83
|
rt300@34
|
84
|
rt300@34
|
85 // whole process unit test(?)
|
rt300@34
|
86
|
rt300@35
|
87 unsigned long long inCoord = 982635920;
|
rt300@34
|
88 vector<int> checkParam;
|
rt300@34
|
89 checkParam = calculateParamsFromIndex(inCoord);
|
rt300@34
|
90 cout << "PARAMS: ";
|
rt300@34
|
91 for(int i=0; i<P;i++){
|
rt300@34
|
92 cout << checkParam[i] << " ";
|
rt300@34
|
93 }
|
rt300@34
|
94 cout << '\n';
|
rt300@34
|
95 unsigned long long outCoord = calculateIndexFromParams(checkParam);
|
rt300@35
|
96 cout << "UNIT TEST COORD = " << outCoord << "\n";
|
rt300@34
|
97
|
rt300@34
|
98
|
rt300@34
|
99 }
|
rt300@34
|
100 //--------------------------------------------------------------
|
rt300@34
|
101 void Hilbert::makeRotationRules(){
|
rt300@34
|
102
|
rt300@34
|
103
|
rt300@34
|
104 }
|
rt300@34
|
105 //--------------------------------------------------------------
|
rt300@34
|
106 void Hilbert::makeCode(){
|
rt300@34
|
107
|
rt300@34
|
108 ////////////////////////////////~~~,,,,,,,,.........
|
rt300@34
|
109 ////////// gray code generation
|
rt300@34
|
110 //////
|
rt300@34
|
111 ///
|
rt300@34
|
112 // TODO 5bit specific
|
rt300@34
|
113 // transition sequence.... what a palaver! only need to do once though.
|
rt300@34
|
114
|
rt300@34
|
115 // ALL SEQUENCES SHOULD END WITH 1 0 0 ... i.e. MSB changes
|
rt300@34
|
116
|
rt300@34
|
117 // max MRL 5 digit
|
rt300@34
|
118 int trans[] = {2,3,4,0,2,1,4,3,2,0,4,3,2,1,4,0,2,3,4,0,2,1,4,3,2,0,4,3,2,1,4,0};
|
rt300@34
|
119
|
rt300@34
|
120 // other ones
|
rt300@34
|
121 int trans4[] = {1,2,0,3,0,2,1,2,3,0,3,2,1,3,1,0};
|
rt300@34
|
122 int trans3[] = {1,2,1,0,1,2,1,0};
|
rt300@34
|
123 int trans2[] = {1,0,1,0};
|
rt300@34
|
124 // balanced 5
|
rt300@34
|
125 int transB5[] = {1,2,3,4,5,1,5,2,3,5,2,4,2,3,4,1,4,3,2,3,1,5,3,4,1,5,2,5,3,4,1,3};
|
rt300@34
|
126
|
rt300@34
|
127
|
rt300@34
|
128 if(P == 5){
|
rt300@34
|
129
|
rt300@34
|
130 }else{
|
rt300@34
|
131 cout << "Only 5 digit supported now\n";
|
rt300@34
|
132 return;
|
rt300@34
|
133 }
|
rt300@34
|
134
|
rt300@34
|
135
|
rt300@34
|
136 int code[codeLength][P]; // start with normal array
|
rt300@34
|
137 for(int j=0; j<P; j++){
|
rt300@34
|
138 code[0][j] = 0;
|
rt300@34
|
139 }
|
rt300@34
|
140
|
rt300@34
|
141 for(int i=0; i < codeLength-1; i++){ // don't need last 3
|
rt300@34
|
142 for(int j=0; j<P; j++){
|
rt300@34
|
143 if (j == abs(trans[i])){
|
rt300@34
|
144 code[i+1][j] = !code[i][j];
|
rt300@34
|
145 }else{
|
rt300@34
|
146 code[i+1][j] = code[i][j];
|
rt300@34
|
147 }
|
rt300@34
|
148 }
|
rt300@34
|
149
|
rt300@34
|
150 }
|
rt300@34
|
151
|
rt300@34
|
152 for(int i=0; i < codeLength; i++){ // fill vector
|
rt300@34
|
153
|
rt300@34
|
154 theGrayCode.push_back(vector<bool>());
|
rt300@34
|
155
|
rt300@34
|
156 for(int j=0; j<P; j++){
|
rt300@34
|
157 theGrayCode[i].push_back(code[i][j]);
|
rt300@34
|
158
|
rt300@34
|
159 }
|
rt300@34
|
160
|
rt300@34
|
161 }
|
rt300@34
|
162
|
rt300@34
|
163 for(int i=0;i<codeLength;i++){
|
rt300@34
|
164 theGrayCodeD.push_back(bin2dec(theGrayCode[i]));
|
rt300@34
|
165 //cout << bin2dec(theGrayCode[i]) << "\n";
|
rt300@34
|
166 //cout << theGrayCodeD[i] << "\n";
|
rt300@34
|
167 }
|
rt300@34
|
168
|
rt300@34
|
169 }
|
rt300@34
|
170
|
rt300@34
|
171 //--------------------------------------------------------------------
|
rt300@34
|
172 //--------------------------------------------------------------------
|
rt300@34
|
173 vector<int> Hilbert::calculateParamsFromIndex(unsigned long long index){
|
rt300@34
|
174 // set %set start and direction of biggest cube
|
rt300@34
|
175 unsigned int i;
|
rt300@34
|
176 unsigned int mask = (1 << P) - 1;
|
rt300@34
|
177 unsigned int entryPoint = 0;
|
rt300@34
|
178 int direction = P - 1, directionInv = P-1;
|
rt300@34
|
179 unsigned int vertex=0, subindex=0, newe=0, newd=0, entryPointInv = 0;
|
rt300@34
|
180
|
rt300@34
|
181 unsigned int pbin[N];
|
rt300@34
|
182
|
rt300@34
|
183 vector<int> params;
|
rt300@34
|
184 for(i=0; i<P;i++){
|
rt300@34
|
185 params.push_back(0);
|
rt300@34
|
186 }
|
rt300@34
|
187 for(i=0;i<N;i++){
|
rt300@34
|
188 pbin[i] = 0;
|
rt300@34
|
189
|
rt300@34
|
190 }
|
rt300@34
|
191
|
rt300@34
|
192 // for loop thru bit levels
|
rt300@36
|
193
|
rt300@34
|
194 for(int blev = N-1; blev >= 0; blev--){
|
rt300@34
|
195 // get next highest bits of index
|
rt300@34
|
196
|
rt300@34
|
197 subindex = index >> blev*P;
|
rt300@34
|
198 subindex = mask & subindex;
|
rt300@34
|
199
|
rt300@35
|
200 //cout << "subindex: " << subindex << "\n";
|
rt300@34
|
201 // which vertex corresponds to this index?
|
rt300@34
|
202 if(subindex < theGrayCode.size()) vertex = theGrayCodeD[subindex];
|
rt300@34
|
203
|
rt300@35
|
204 //cout << "rotated: " << vertex << "\n";
|
rt300@34
|
205
|
rt300@34
|
206
|
rt300@34
|
207 //% inverse rotate this T_e,d
|
rt300@34
|
208 entryPointInv = rotr(entryPoint, direction+1,P);
|
rt300@34
|
209 directionInv = (P - direction - 2);
|
rt300@34
|
210 vertex = rotate(vertex,entryPointInv,directionInv);
|
rt300@34
|
211
|
rt300@35
|
212 //cout << "vertex: " << vertex << "\n";
|
rt300@34
|
213
|
rt300@34
|
214 //% build up bits of params
|
rt300@34
|
215
|
rt300@34
|
216 for(int j=0;j<P;j++){
|
rt300@34
|
217
|
rt300@34
|
218 params[j] += ((vertex >> P-j-1) & 1)*(1<<blev);
|
rt300@34
|
219 }
|
rt300@34
|
220 if(P==5){
|
rt300@34
|
221 newe = entryVertices5[subindex];
|
rt300@34
|
222 newd = rotations5[subindex];
|
rt300@34
|
223 }
|
rt300@34
|
224 // next rotations...
|
rt300@34
|
225 int lse2 = rotr(newe, P-direction-1,P);
|
rt300@34
|
226 entryPoint = entryPoint ^ lse2;
|
rt300@34
|
227 direction = (direction + newd + 1) % P;
|
rt300@35
|
228 //cout << "------------\n";
|
rt300@34
|
229 }
|
rt300@34
|
230
|
rt300@34
|
231 return params;
|
rt300@34
|
232 }
|
rt300@34
|
233 //--------------------------------------------------------------------
|
rt300@34
|
234 int Hilbert::rotate(int vertex, int entryPoint, int direction){
|
rt300@34
|
235 vertex = vertex ^ entryPoint;
|
rt300@34
|
236 vertex = rotr(vertex,direction+1,P);
|
rt300@34
|
237 return vertex;
|
rt300@34
|
238 }
|
rt300@34
|
239 //--------------------------------------------------------------------
|
rt300@34
|
240
|
rt300@34
|
241
|
rt300@34
|
242 //--------------------------------------------------------------------
|
rt300@35
|
243 unsigned long long Hilbert::calculateIndexFromParams(vector<int> params){
|
rt300@34
|
244 // set %set start and direction of biggest cube
|
rt300@34
|
245
|
rt300@35
|
246 unsigned long long h = 0;
|
rt300@34
|
247
|
rt300@34
|
248 unsigned int i;
|
rt300@34
|
249 unsigned int entryPoint = 0;
|
rt300@34
|
250 int direction = P - 1;
|
rt300@35
|
251 unsigned int vertex=0, newe=0, newd=0;
|
rt300@35
|
252 unsigned long long subindex; // needs to be long cos we lewt shift by alot
|
rt300@34
|
253 vector<int> subindices;
|
rt300@34
|
254
|
rt300@34
|
255
|
rt300@34
|
256 // for loop thru bit levels
|
rt300@36
|
257
|
rt300@34
|
258 for(int blev = N-1; blev >= 0; blev--){
|
rt300@34
|
259 //% get next highest bit of param
|
rt300@34
|
260
|
rt300@34
|
261 vector<bool> vertexb;
|
rt300@34
|
262 for(i=0;i<P;i++){
|
rt300@34
|
263 vertexb.push_back((params[i]/(1 << blev)) % 2);
|
rt300@34
|
264
|
rt300@34
|
265 }
|
rt300@35
|
266 //printBool(vertexb);
|
rt300@34
|
267 vertex = bin2dec(vertexb);
|
rt300@35
|
268 //cout << "INV vertex: " << vertex << "\n";
|
rt300@34
|
269 // rotate it
|
rt300@34
|
270 vertex = rotate(vertex,entryPoint,direction);
|
rt300@34
|
271
|
rt300@34
|
272 // % get new index of thisun
|
rt300@35
|
273 //cout << "INV rotated: " << vertex << "\n";
|
rt300@34
|
274
|
rt300@34
|
275
|
rt300@34
|
276 for(i=0;i<codeLength;i++){
|
rt300@34
|
277 if(vertex == theGrayCodeD[i]){
|
rt300@34
|
278 subindex = i;
|
rt300@34
|
279 }
|
rt300@34
|
280 }
|
rt300@35
|
281 //cout << "INV subindex: " << subindex << "\n";
|
rt300@34
|
282 // work out next rotations
|
rt300@34
|
283 newe = entryVertices5[subindex];
|
rt300@34
|
284 newd = rotations5[subindex];
|
rt300@34
|
285
|
rt300@34
|
286 // next rotations...
|
rt300@34
|
287 int lse2 = rotr(newe, P-direction-1,P);
|
rt300@34
|
288 entryPoint = entryPoint ^ lse2;
|
rt300@34
|
289 direction = (direction + newd + 1) % P;
|
rt300@34
|
290
|
rt300@34
|
291 // now build up h from subindices
|
rt300@34
|
292 h += subindex << (blev*P);
|
rt300@35
|
293 //cout << "------------\n";
|
rt300@34
|
294 }
|
rt300@34
|
295
|
rt300@34
|
296 return h;
|
rt300@34
|
297 }
|
rt300@34
|
298 //--------------------------------------------------------------------
|