comparison hilbert.cpp @ 34:94df2cd72d7b

New hilbert code - does rotations!

+//
+//  hilbert.cpp
+//  sonicZoom
+//
+//  Created by Robert Tubb on 04/04/2013.
+//
+//
+
+#include "hilbert.h"
+Hilbert hilbert;
+
+//--------------------------------------------------------------------
+unsigned int rotl(unsigned int value, int shift, int digits) {
+    unsigned int ones = (1 << digits) - 1;
+    if(shift >= digits)
+        shift -= digits;
+    
+    unsigned int out = (value << shift) | (value >> (digits - shift));
+    return out & ones;
+}
+//--------------------------------------------------------------------
+unsigned int rotr(unsigned int value, int shift, int digits) {
+    int ones = (1 << digits) - 1;
+    // can only handle -digits to 2*digits-1
+    if(shift >= digits) shift = shift - digits;
+    if(shift < 0) shift = digits + shift;
+    unsigned int out = (value >> shift) | (value << (digits - shift));
+    return out & ones;
+}
+//--------------------------------------------------------------------
+unsigned int bin2dec(vector<bool> binary){
+    // takes in a vector of booleans. we assume that the leftmost bool (highest index) is the MSB
+    unsigned int dec = 0;
+    int i,pwr;
+    int N = binary.size();
+    
+    for (i = 0; i<N; i++){
+        pwr = N-i-1;
+        dec += binary[i]*(1 << pwr);
+        
+    }
+    
+    return dec;
+}
+//--------------------------------------------------------------------
+vector<bool> dec2bin(unsigned int number, unsigned int size){
+    vector<bool> bin;
+    int i, pwr;
+    for(i=0;i<size;i++){
+        bin.push_back(0);
+    }
+    if(number >= (1 << size)){
+        cout << " bad size for Hilbert::dec2bin" << "\n";
+        return bin;
+    }
+    for(i=0;i<size;i++){
+        pwr = size-i;
+        bin[i] = (number/(1 << pwr)) % 2;
+        bin.push_back(0);
+    }
+    
+    
+    
+}
+void printBool(vector <bool> vcode){
+    vector<bool>::iterator bit;
+
+        for(bit=vcode.begin(); bit!=vcode.end() ; bit++){
+            
+            cout  << *bit;
+        }
+        cout << "\n";
+
+}
+//--------------------------------------------------------------------
+void Hilbert::init(int aN, int aP){
+    N = aN;
+    P = aP;
+    codeLength = (int)pow(2.0,P);
+    makeCode();
+    makeRotationRules();
+    
+    
+    // whole process unit test(?)
+    
+    unsigned long long inCoord = 888999777;
+    vector<int> checkParam;
+    checkParam = calculateParamsFromIndex(inCoord);
+    cout << "PARAMS: ";
+    for(int i=0; i<P;i++){
+        cout << checkParam[i] << " ";
+    }
+    cout << '\n';
+    unsigned long long outCoord = calculateIndexFromParams(checkParam);
+    cout << "Unit TEst coord = " << outCoord << "\n";
+    
+    
+}
+//--------------------------------------------------------------
+void Hilbert::makeRotationRules(){
+    
+    
+}
+//--------------------------------------------------------------
+void Hilbert::makeCode(){
+    
+    ////////////////////////////////~~~,,,,,,,,.........
+    ////////// gray code generation
+    //////
+    ///
+    // TODO 5bit specific
+    // transition sequence.... what a palaver! only need to do once though.
+    
+    // ALL SEQUENCES SHOULD END WITH  1 0 0 ... i.e. MSB changes
+    
+    // max MRL 5 digit
+    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};
+    
+    // other ones
+    int trans4[] = {1,2,0,3,0,2,1,2,3,0,3,2,1,3,1,0};
+    int trans3[] = {1,2,1,0,1,2,1,0};
+    int trans2[] = {1,0,1,0};
+    // balanced 5
+    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};
+    
+    
+    if(P == 5){
+        
+    }else{
+        cout << "Only 5 digit supported now\n";
+        return;
+    }
+    
+    
+    int code[codeLength][P]; // start with normal array
+    for(int j=0; j<P; j++){
+        code[0][j] = 0;
+    }
+    
+    for(int i=0; i < codeLength-1; i++){  // don't need last 3
+        for(int j=0; j<P; j++){
+            if (j == abs(trans[i])){
+                code[i+1][j] = !code[i][j];
+            }else{
+                code[i+1][j] = code[i][j];
+            }
+        }
+        
+    }
+    
+    for(int i=0; i < codeLength; i++){ // fill vector
+        
+        theGrayCode.push_back(vector<bool>());
+        
+        for(int j=0; j<P; j++){
+            theGrayCode[i].push_back(code[i][j]);
+
+        }
+
+    }
+
+    for(int i=0;i<codeLength;i++){
+        theGrayCodeD.push_back(bin2dec(theGrayCode[i]));
+        //cout << bin2dec(theGrayCode[i]) << "\n";
+        //cout << theGrayCodeD[i] << "\n";
+    }
+  
+}
+
+//--------------------------------------------------------------------
+//--------------------------------------------------------------------
+vector<int> Hilbert::calculateParamsFromIndex(unsigned long long index){
+    // set  %set start and direction of biggest cube
+    unsigned int i;
+    unsigned int mask = (1 << P) - 1;
+    unsigned int entryPoint = 0;
+    int direction = P - 1, directionInv = P-1;
+    unsigned int vertex=0, subindex=0, newe=0, newd=0, entryPointInv = 0;
+    
+    unsigned int pbin[N];
+
+    vector<int> params;
+    for(i=0; i<P;i++){
+        params.push_back(0);
+    }
+    for(i=0;i<N;i++){
+        pbin[i] = 0;
+
+    }
+    
+    // for loop thru bit levels
+    i=0;
+    for(int blev = N-1; blev >= 0; blev--){
+        // get next highest bits of index
+        
+        subindex = index >> blev*P;
+        subindex = mask & subindex;
+        
+        cout << "subindex: " << subindex << "\n";
+        // which vertex corresponds to this index?
+        if(subindex < theGrayCode.size()) vertex = theGrayCodeD[subindex];
+
+        cout << "rotated: " << vertex << "\n";
+ 
+        
+        //% inverse rotate this T_e,d
+        entryPointInv = rotr(entryPoint, direction+1,P);
+        directionInv = (P - direction - 2);
+        vertex = rotate(vertex,entryPointInv,directionInv);
+
+        cout << "vertex: " << vertex << "\n";
+        
+        //% build up bits of params
+     
+        for(int j=0;j<P;j++){
+            
+            params[j] += ((vertex >> P-j-1) & 1)*(1<<blev);
+        }
+        if(P==5){
+            newe = entryVertices5[subindex];
+            newd = rotations5[subindex];
+        }
+        // next rotations...
+        int lse2 = rotr(newe, P-direction-1,P);
+        entryPoint = entryPoint ^ lse2;
+        direction = (direction + newd + 1) % P;
+        cout << "------------\n";
+    }
+
+    return params;
+}
+//--------------------------------------------------------------------
+int Hilbert::rotate(int vertex, int entryPoint, int direction){
+    vertex = vertex ^ entryPoint;
+    vertex = rotr(vertex,direction+1,P);
+    return vertex;
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+long long Hilbert::calculateIndexFromParams(vector<int> params){
+    // set  %set start and direction of biggest cube
+
+    long long h = 0;
+    
+    unsigned int i;
+    unsigned int entryPoint = 0;
+    int direction = P - 1;
+    unsigned int vertex=0, subindex=0, newe=0, newd=0;
+    
+    vector<int> subindices;
+
+    
+    // for loop thru bit levels
+    i=0;
+    for(int blev = N-1; blev >= 0; blev--){
+        //% get next highest bit of param
+
+        vector<bool> vertexb;
+        for(i=0;i<P;i++){
+            vertexb.push_back((params[i]/(1 << blev)) % 2);
+            
+        }
+        printBool(vertexb);
+        vertex = bin2dec(vertexb);
+        cout << "INV vertex: " << vertex << "\n";
+        // rotate it
+        vertex = rotate(vertex,entryPoint,direction);
+        
+        //    % get new index of thisun
+        cout << "INV rotated: " << vertex << "\n";
+        
+        
+        for(i=0;i<codeLength;i++){
+            if(vertex == theGrayCodeD[i]){
+                subindex = i;
+            }
+        }
+        cout << "INV subindex: " << subindex << "\n";
+        // work out next rotations
+        newe = entryVertices5[subindex];
+        newd = rotations5[subindex];
+        
+        // next rotations...
+        int lse2 = rotr(newe, P-direction-1,P);
+        entryPoint = entryPoint ^ lse2;
+        direction = (direction + newd + 1) % P;
+        
+        // now build up h from subindices
+        h += subindex << (blev*P);
+        cout << "------------\n";
+    }
+
+    return h;
+}
+//--------------------------------------------------------------------