view hilbert.cpp @ 35:790939017078

V0.5. Now uses new continuous hilbert curve. HOORAY.

//
//  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);
    }
    return bin;
    
    
}
//--------------------------------------------------------------------
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 = 982635920;
    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;
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
unsigned long long Hilbert::calculateIndexFromParams(vector<int> params){
    // set  %set start and direction of biggest cube

    unsigned long long h = 0;
    
    unsigned int i;
    unsigned int entryPoint = 0;
    int direction = P - 1;
    unsigned int vertex=0,  newe=0, newd=0;
    unsigned long long subindex; // needs to be long cos we lewt shift by alot
    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;
}
//--------------------------------------------------------------------