HilbertCoordinateMapper.cpp

Go to the documentation of this file.

/*
 *  Copyright (C) 2007  A.N. Yzelman
 *  Released under LGPL, see license.txt
 *
 *  Last modified at 23th of May, 2007, by A.N. Yzelman
 *
 *  HilbertCoordinateMapper.cpp: implementation of HilbertCoordinateMapper.h
 *
 */

#include "HilbertCoordinateMapper.h"
#include <iostream>

void HilbertCoordinateMapper::makeVec( const double x, const double y, vector< double >* where ) {
        (*where)[ 0 ] = x;
        (*where)[ 1 ] = y;
}

void HilbertCoordinateMapper::makeVec( const double x, const double y, const double z, vector< double >* where ) {
        (*where)[ 0 ] = x;
        (*where)[ 1 ] = y;
        (*where)[ 2 ] = z;
}

double HilbertCoordinateMapper::dist( const vector< double > *vx, const vector< double > *vy, const unsigned int dim ) {

        double ret = ( (*vx)[ 0 ] - (*vy)[ 0 ] ); ret *= ret;
        for ( unsigned int i=1; i<dim; i++ ) {
                const double cur = ( (*vx)[ i ] - (*vy)[ i ] );
                ret += cur*cur;
        }
        
        return sqrt( ret );
}

unsigned int HilbertCoordinateMapper::closestTo( const vector< double > *coor, const vector< vector< double > > *corners, const unsigned int dim ) {

        if (corners->size() == 0 )
                cerr << "Empty corners vector passed to HilbertCoordinateMapper::closestTo()!" << endl;
        
        unsigned int ret = 0;
        double min = dist( coor, &( corners->at( 0 ) ), dim );
        
        
        for ( unsigned int i=1 ; i<corners->size(); i++ ) {
                const double curcost = dist( coor, & ( corners->at( i ) ), dim );
                if ( curcost < min ) {
                        min = curcost;
                        ret = i;
                }
        }
        
        return ret;
}

void HilbertCoordinateMapper::mid( const vector< double > *x, const vector< double > *y, vector< double > *to ) {
        
        double* toArr = & ( (*to)[ 0 ] );
        for( unsigned int i=0; i<x->size(); i++ )
                toArr[ i ] = ( (*x)[ i ] + (*y)[ i ] ) / 2;

}

inline
void HilbertCoordinateMapper::mid3D( unsigned int pos, unsigned int i, unsigned int j ) {
        //multi_array<double, 1> ret(extents[3]);
        for( unsigned int d=0; d<3; d++ )
                Box_3Dtemp[ pos ][ d ] = ( Box_3D[ i ][ d ] + Box_3D[ j ][ d ] ) / 2.0;
//return ret;
}

void HilbertCoordinateMapper::refineCoors2D( vector< vector< double > > *corners, const unsigned int cell ) {

        vector< double > a( 2 );
        vector< double > b( 2 );
        vector< double > c( 2 );
        vector< double > d( 2 );

        vector< double > *in = &( ( *corners )[ 0 ] );

        //assume order of the vector to be alphabetical
        switch( cell ) {
                case 1:
                        a = in[ 0 ];
                        mid( &in[ 0 ], &in[ 3 ], &b );
                        mid( &in[ 0 ], &in[ 2 ], &c );
                        mid( &in[ 0 ], &in[ 1 ], &d );
                        break;
                case 2:
                        mid( &in[ 0 ], &in[ 1 ], &a );
                        b = in[ 1 ];
                        mid( &in[ 1 ], &in[ 2 ], &c );
                        mid( &in[ 1 ], &in[ 3 ], &d );
                        break;
                case 3:
                        mid( &in[ 0 ], &in[ 2 ], &a );
                        mid( &in[ 1 ], &in[ 2 ], &b );
                        c = in[ 2 ];
                        mid( &in[ 2 ], &in[ 3 ], &d );
                        break;
                case 4:
                        mid( &in[ 3 ], &in[ 2 ], &a );
                        mid( &in[ 1 ], &in[ 3 ], &b );
                        mid( &in[ 3 ], &in[ 0 ], &c );
                        d = in[ 3 ];
                        break;
                default:
                        cerr << "Invalid cell no. given in refineHilbert2D!" << endl;
                        exit( 1 );
        }
        
        /*corners->clear();
        corners->push_back( a );
        corners->push_back( b );
        corners->push_back( c );
        corners->push_back( d );*/
        (*corners)[ 0 ] = a;
        (*corners)[ 1 ] = b;
        (*corners)[ 2 ] = c;
        (*corners)[ 3 ] = d;

}

void HilbertCoordinateMapper::refineCoors3D( const unsigned int cell ) {

        //assume order of the vector to be alphabetical
        switch( cell ) {
                case 1:
                        Box_3Dtemp[ 0 ] = Box_3D[ 0 ];
                        mid3D( 1, 0, 3 );
                        mid3D( 2, 0, 4 );
                        mid3D( 3, 0, 7 );
                        mid3D( 4, 0, 6 );
                        mid3D( 5, 0, 5 );
                        mid3D( 6, 0, 2 );
                        mid3D( 7, 0, 1 );
                        break;
                case 2:
                        mid3D( 0, 0, 1 );
                        mid3D( 1, 0, 6 );
                        mid3D( 2, 1, 6 );
                        Box_3Dtemp[ 3 ] = Box_3D[ 1 ];
                        mid3D( 4, 1, 2 );
                        mid3D( 5, 1, 5 );
                        mid3D( 6, 0, 5 );
                        mid3D( 7, 0, 2 );
                        break;
                case 3:
                        mid3D( 0, 0, 2 );
                        mid3D( 1, 0, 5 );
                        mid3D( 2, 1, 5 );
                        mid3D( 3, 1, 2 );
                        Box_3Dtemp[ 4 ] = Box_3D[ 2 ];
                        mid3D( 5, 2, 5 );
                        mid3D( 6, 3, 5 );
                        mid3D( 7, 3, 2 );
                        break;
                case 4:
                        mid3D( 0, 3, 2 );
                        Box_3Dtemp[ 1 ] = Box_3D[ 3 ];
                        mid3D( 2, 0, 3 );
                        mid3D( 3, 0, 2 );
                        mid3D( 4, 0, 5 );
                        mid3D( 5, 0, 4 );
                        mid3D( 6, 3, 4 );
                        mid3D( 7, 3, 5 );
                        break;
                case 5:
                        mid3D( 0, 3, 5 );
                        mid3D( 1, 3, 4 );
                        mid3D( 2, 0, 4 );
                        mid3D( 3, 0, 5 );
                        mid3D( 4, 7, 5 );
                        mid3D( 5, 4, 7 );
                        Box_3Dtemp[ 6 ] = Box_3D[ 4 ];
                        mid3D( 7, 4, 5 );
                        break;
                case 6:
                        mid3D( 0, 4, 5 );
                        mid3D( 1, 3, 5 );
                        mid3D( 2, 2, 5 );
                        Box_3Dtemp[ 3 ] = Box_3D[ 5 ];
                        mid3D( 4, 6, 5 );
                        mid3D( 5, 1, 5 );
                        mid3D( 6, 0, 5 );
                        mid3D( 7, 7, 5 );
                        break;
                case 7:
                        mid3D( 0, 7, 5 );
                        mid3D( 1, 0, 5 );
                        mid3D( 2, 2, 5 );
                        mid3D( 3, 6, 5 );
                        Box_3Dtemp[ 4 ] = Box_3D[ 6 ];
                        mid3D( 5, 1, 6 );
                        mid3D( 6, 0, 6 );
                        mid3D( 7, 7, 6 );
                        break;
                case 8:
                        mid3D( 0, 7, 6 );
                        mid3D( 1, 7, 5 );
                        mid3D( 2, 0, 5 );
                        mid3D( 3, 6, 0 );
                        mid3D( 4, 0, 7 );
                        mid3D( 5, 0, 5 );
                        mid3D( 6, 7, 5 );
                        Box_3Dtemp[ 7 ] = Box_3D[ 7 ];
                        break;
                default:
                        cerr << "Invalid cell no. given in refineHilbert3D!" << endl;
                        exit( 1 );
        }
                
        Box_3D = Box_3Dtemp;
}

double HilbertCoordinateMapper::toHilbert2D( const vector< double > *coordinates ) {
        vector< unsigned int > digits;
        vector< vector< double > > Box;
        
        //note: remember that the order matters!
        vector< double > a( 2 ); makeVec( 0, 0, &a); Box.push_back( a );
        vector< double > b( 2 ); makeVec( 0, 1, &b); Box.push_back( b );
        vector< double > c( 2 ); makeVec( 1, 1, &c); Box.push_back( c );
        vector< double > d( 2 ); makeVec( 1, 0, &d); Box.push_back( d );
        
        double curDist;
        unsigned int closest = closestTo( coordinates, &Box, 2 );
        digits.push_back( closest );
        
        for ( unsigned int i=0; i<RECURSE_LEVEL; i++ ) {
        
                refineCoors2D( &Box, closest+1 );
        
                closest = closestTo( coordinates, &Box, 2 );
                
                digits.push_back( closest );
                
                curDist = dist( coordinates, &Box.at( closest ), 2 );   
        
        }

        //const double divBy = pow( static_cast< double >( 4 ), static_cast< int >( RECURSE_LEVEL + 1 ) ) - 1;
        //return base4_to_base10( &digits ) / divBy;    
        return base4_to_base10( &digits );      
}

double HilbertCoordinateMapper::toHilbert3D( const vector< double > *coordinates ) {
        vector< unsigned int > digits;
        Box_3D = pattern;
        
        unsigned int closest = closestTo3D( coordinates );
        digits.push_back( closest );
        
        for ( unsigned int i=0; i<RECURSE_LEVEL; i++ ) {
        
                refineCoors3D( closest + 1 );
        
                closest = closestTo3D( coordinates );
                
                digits.push_back( closest );

        }
        
        //const double divBy = pow( static_cast< double >( 8 ), static_cast< int >( RECURSE_LEVEL + 1 ) ) - 1;
#ifdef _DEBUG
        //cout << "bs: " ;
        //for( unsigned int i=0; i<digits.size() ; i++ )
        //      cout << digits[ i ];
        //cout << endl;
        //cout << "up: " << base8_to_base10( &digits ) << endl;
        //cout << "dw: " << divBy << endl;
#endif
        //return base8_to_base10( &digits ) / divBy;
        return base8_to_base10( &digits );
}

unsigned int HilbertCoordinateMapper::closestTo3D( const vector< double > *coor ) {
        
        unsigned int ret = 0;
        double min = dist3D( coor, 0 );
        
        for ( unsigned int i=1; i<8; i++ ) {
                const double curcost = dist3D( coor, i );
                if ( curcost < min ) {
                        min = curcost;
                        ret = i;
                }
        }
        
        return ret;
}

inline double HilbertCoordinateMapper::dist3D( const vector< double > *vx, const unsigned int num ) {

        double ret = ( (*vx)[ 0 ] - Box_3D[ num ][ 0 ] ); ret *= ret;
        for ( unsigned int i=1; i<3; i++ ) {
                const double cur = ( (*vx)[ i ] - Box_3D[ num ][ i ] );
                ret += cur*cur;
        }
        
        //return sqrt( ret );
        return ret;
}

typedef boost::multi_array<double, 2 > array_type;   

array_type HilbertCoordinateMapper::pattern = HilbertCoordinateMapper::makePattern();

array_type HilbertCoordinateMapper::Box_3D = HilbertCoordinateMapper::pattern;

array_type HilbertCoordinateMapper::Box_3Dtemp = HilbertCoordinateMapper::pattern;

---