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