annotate hilbert.cpp @ 49:178642d134a7 tip

xtra files
author Robert Tubb <rt300@eecs.qmul.ac.uk>
date Wed, 01 May 2013 17:34:33 +0100
parents b91a1859829a
children
rev   line source
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
rt300@43 11 // UTILS
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@43 77 //====================================================================
rt300@43 78 //--------------------------------------------------------------------
rt300@43 79
rt300@43 80 Hilbert::Hilbert(){
rt300@43 81 N = 7;
rt300@43 82 P = 5;
rt300@43 83 codeLength = (int)pow(2.0,P);
rt300@43 84
rt300@43 85
rt300@43 86 makeCode();
rt300@43 87 makeRotationRules();
rt300@43 88
rt300@43 89 }
rt300@43 90
rt300@43 91 //--------------------------------------------------------------------
rt300@43 92
rt300@43 93 Hilbert::Hilbert(int aN, int aP){
rt300@34 94 N = aN;
rt300@34 95 P = aP;
rt300@34 96 codeLength = (int)pow(2.0,P);
rt300@43 97
rt300@43 98
rt300@34 99 makeCode();
rt300@34 100 makeRotationRules();
rt300@34 101
rt300@43 102 }
rt300@43 103 //--------------------------------------------------------------------
rt300@43 104
rt300@43 105 void Hilbert::changeCurve(int aN, int aP){
rt300@43 106 N = aN;
rt300@43 107 P = aP;
rt300@43 108 codeLength = (int)pow(2.0,P);
rt300@34 109
rt300@34 110
rt300@43 111 makeCode();
rt300@43 112 makeRotationRules();
rt300@34 113
rt300@34 114 }
rt300@34 115 //--------------------------------------------------------------
rt300@34 116 void Hilbert::makeRotationRules(){
rt300@43 117 theRotations.clear();
rt300@43 118 theEntryVertices.clear();
rt300@34 119
rt300@43 120 const int *rp;
rt300@43 121 const int *ep;
rt300@43 122 if(P == 5){
rt300@43 123 rp = rotations5;
rt300@43 124 ep = entryVertices5;
rt300@43 125 }else if(P == 4){
rt300@43 126 rp = rotations4;
rt300@43 127 ep = entryVertices4;
rt300@43 128 }else if(P == 3){
rt300@43 129 rp = rotations3;
rt300@43 130 ep = entryVertices3;
rt300@43 131 }else if(P == 2){
rt300@43 132 rp = rotations2;
rt300@43 133 ep = entryVertices2;
rt300@43 134 }else if(P == 1){ // will this even work?
rt300@43 135 rp = rotations1;
rt300@43 136 ep = entryVertices1;
rt300@43 137 }else{
rt300@43 138 cout << "ERROR: makeRotationRules: bad numb of MIDI params";
rt300@43 139 return;
rt300@43 140 }
rt300@43 141
rt300@43 142 for(int i=0; i < codeLength; i++){ // fill vector
rt300@43 143
rt300@43 144 theEntryVertices.push_back(ep[i]);
rt300@43 145 theRotations.push_back(rp[i]);
rt300@43 146
rt300@43 147 }
rt300@34 148
rt300@34 149 }
rt300@34 150 //--------------------------------------------------------------
rt300@34 151 void Hilbert::makeCode(){
rt300@34 152
rt300@34 153 ////////////////////////////////~~~,,,,,,,,.........
rt300@34 154 ////////// gray code generation
rt300@34 155 //////
rt300@34 156 ///
rt300@43 157 // choose the gray code sequence and rotation based on the numnber of dimensions P
rt300@34 158
rt300@34 159 // ALL SEQUENCES SHOULD END WITH 1 0 0 ... i.e. MSB changes
rt300@34 160
rt300@34 161 // max MRL 5 digit
rt300@43 162
rt300@43 163 theGrayCode.clear();
rt300@43 164 theGrayCodeD.clear();
rt300@43 165
rt300@43 166 int trans5[] = {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 167
rt300@34 168 // other ones
rt300@34 169 int trans4[] = {1,2,0,3,0,2,1,2,3,0,3,2,1,3,1,0};
rt300@34 170 int trans3[] = {1,2,1,0,1,2,1,0};
rt300@34 171 int trans2[] = {1,0,1,0};
rt300@43 172 int trans1[] = {0,0};
rt300@34 173 // balanced 5
rt300@34 174 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 175
rt300@43 176 int *tp;
rt300@43 177
rt300@34 178
rt300@34 179 if(P == 5){
rt300@43 180 tp = trans5;
rt300@43 181 }else if(P == 4){
rt300@43 182 tp = trans4;
rt300@43 183 }else if(P == 3){
rt300@43 184 tp = trans3;
rt300@43 185 }else if(P == 2){
rt300@43 186 tp = trans2;
rt300@43 187 }else if(P == 1){ // will this even work?
rt300@43 188 tp = trans1;
rt300@34 189 }
rt300@43 190
rt300@34 191
rt300@34 192 int code[codeLength][P]; // start with normal array
rt300@34 193 for(int j=0; j<P; j++){
rt300@34 194 code[0][j] = 0;
rt300@34 195 }
rt300@34 196
rt300@34 197 for(int i=0; i < codeLength-1; i++){ // don't need last 3
rt300@34 198 for(int j=0; j<P; j++){
rt300@43 199 if (j == abs(tp[i])){
rt300@34 200 code[i+1][j] = !code[i][j];
rt300@34 201 }else{
rt300@34 202 code[i+1][j] = code[i][j];
rt300@34 203 }
rt300@34 204 }
rt300@34 205
rt300@34 206 }
rt300@34 207
rt300@34 208 for(int i=0; i < codeLength; i++){ // fill vector
rt300@34 209
rt300@34 210 theGrayCode.push_back(vector<bool>());
rt300@34 211
rt300@34 212 for(int j=0; j<P; j++){
rt300@34 213 theGrayCode[i].push_back(code[i][j]);
rt300@34 214
rt300@34 215 }
rt300@34 216
rt300@34 217 }
rt300@34 218
rt300@34 219 for(int i=0;i<codeLength;i++){
rt300@34 220 theGrayCodeD.push_back(bin2dec(theGrayCode[i]));
rt300@34 221 //cout << bin2dec(theGrayCode[i]) << "\n";
rt300@34 222 //cout << theGrayCodeD[i] << "\n";
rt300@34 223 }
rt300@34 224
rt300@34 225 }
rt300@34 226
rt300@34 227 //--------------------------------------------------------------------
rt300@34 228 //--------------------------------------------------------------------
rt300@34 229 vector<int> Hilbert::calculateParamsFromIndex(unsigned long long index){
rt300@34 230 // set %set start and direction of biggest cube
rt300@34 231 unsigned int i;
rt300@34 232 unsigned int mask = (1 << P) - 1;
rt300@34 233 unsigned int entryPoint = 0;
rt300@34 234 int direction = P - 1, directionInv = P-1;
rt300@34 235 unsigned int vertex=0, subindex=0, newe=0, newd=0, entryPointInv = 0;
rt300@34 236
rt300@34 237 unsigned int pbin[N];
rt300@34 238
rt300@34 239 vector<int> params;
rt300@34 240 for(i=0; i<P;i++){
rt300@34 241 params.push_back(0);
rt300@34 242 }
rt300@34 243 for(i=0;i<N;i++){
rt300@34 244 pbin[i] = 0;
rt300@34 245
rt300@34 246 }
rt300@34 247
rt300@34 248 // for loop thru bit levels
rt300@36 249
rt300@34 250 for(int blev = N-1; blev >= 0; blev--){
rt300@34 251 // get next highest bits of index
rt300@34 252
rt300@34 253 subindex = index >> blev*P;
rt300@34 254 subindex = mask & subindex;
rt300@34 255
rt300@35 256 //cout << "subindex: " << subindex << "\n";
rt300@34 257 // which vertex corresponds to this index?
rt300@34 258 if(subindex < theGrayCode.size()) vertex = theGrayCodeD[subindex];
rt300@34 259
rt300@35 260 //cout << "rotated: " << vertex << "\n";
rt300@34 261
rt300@34 262
rt300@34 263 //% inverse rotate this T_e,d
rt300@34 264 entryPointInv = rotr(entryPoint, direction+1,P);
rt300@34 265 directionInv = (P - direction - 2);
rt300@34 266 vertex = rotate(vertex,entryPointInv,directionInv);
rt300@34 267
rt300@35 268 //cout << "vertex: " << vertex << "\n";
rt300@34 269
rt300@34 270 //% build up bits of params
rt300@34 271
rt300@34 272 for(int j=0;j<P;j++){
rt300@34 273
rt300@34 274 params[j] += ((vertex >> P-j-1) & 1)*(1<<blev);
rt300@34 275 }
rt300@43 276 if(subindex > theEntryVertices.size()){
rt300@43 277 cout << "ERROR: subindex too large for theEntryVertices\n";
rt300@34 278 }
rt300@43 279 newe = theEntryVertices[subindex];
rt300@43 280 newd = theRotations[subindex];
rt300@34 281 // next rotations...
rt300@34 282 int lse2 = rotr(newe, P-direction-1,P);
rt300@34 283 entryPoint = entryPoint ^ lse2;
rt300@34 284 direction = (direction + newd + 1) % P;
rt300@35 285 //cout << "------------\n";
rt300@34 286 }
rt300@34 287
rt300@34 288 return params;
rt300@34 289 }
rt300@34 290 //--------------------------------------------------------------------
rt300@34 291 int Hilbert::rotate(int vertex, int entryPoint, int direction){
rt300@34 292 vertex = vertex ^ entryPoint;
rt300@34 293 vertex = rotr(vertex,direction+1,P);
rt300@34 294 return vertex;
rt300@34 295 }
rt300@34 296 //--------------------------------------------------------------------
rt300@34 297
rt300@34 298
rt300@34 299 //--------------------------------------------------------------------
rt300@35 300 unsigned long long Hilbert::calculateIndexFromParams(vector<int> params){
rt300@34 301 // set %set start and direction of biggest cube
rt300@34 302
rt300@35 303 unsigned long long h = 0;
rt300@34 304
rt300@34 305 unsigned int i;
rt300@34 306 unsigned int entryPoint = 0;
rt300@34 307 int direction = P - 1;
rt300@35 308 unsigned int vertex=0, newe=0, newd=0;
rt300@35 309 unsigned long long subindex; // needs to be long cos we lewt shift by alot
rt300@34 310 vector<int> subindices;
rt300@34 311
rt300@34 312
rt300@34 313 // for loop thru bit levels
rt300@36 314
rt300@34 315 for(int blev = N-1; blev >= 0; blev--){
rt300@34 316 //% get next highest bit of param
rt300@34 317
rt300@34 318 vector<bool> vertexb;
rt300@34 319 for(i=0;i<P;i++){
rt300@34 320 vertexb.push_back((params[i]/(1 << blev)) % 2);
rt300@34 321
rt300@34 322 }
rt300@35 323 //printBool(vertexb);
rt300@34 324 vertex = bin2dec(vertexb);
rt300@35 325 //cout << "INV vertex: " << vertex << "\n";
rt300@34 326 // rotate it
rt300@34 327 vertex = rotate(vertex,entryPoint,direction);
rt300@34 328
rt300@34 329 // % get new index of thisun
rt300@35 330 //cout << "INV rotated: " << vertex << "\n";
rt300@34 331
rt300@34 332
rt300@34 333 for(i=0;i<codeLength;i++){
rt300@34 334 if(vertex == theGrayCodeD[i]){
rt300@34 335 subindex = i;
rt300@34 336 }
rt300@34 337 }
rt300@35 338 //cout << "INV subindex: " << subindex << "\n";
rt300@34 339 // work out next rotations
rt300@43 340 if(subindex > theEntryVertices.size()){
rt300@43 341 cout << "ERROR: subindex too large for theEntryVertices\n";
rt300@43 342 }
rt300@43 343 newe = theEntryVertices[subindex];
rt300@43 344 newd = theRotations[subindex];
rt300@34 345
rt300@34 346 // next rotations...
rt300@34 347 int lse2 = rotr(newe, P-direction-1,P);
rt300@34 348 entryPoint = entryPoint ^ lse2;
rt300@34 349 direction = (direction + newd + 1) % P;
rt300@34 350
rt300@34 351 // now build up h from subindices
rt300@34 352 h += subindex << (blev*P);
rt300@35 353 //cout << "------------\n";
rt300@34 354 }
rt300@34 355
rt300@34 356 return h;
rt300@34 357 }
rt300@34 358 //--------------------------------------------------------------------