Morton ordering in 3D usign unit64_t as input

Question

I am trying to use Morton code to produce unique encoding for a given (x,y,z) where x,y,z are double precision floating point numbers. I presume that I can use type cast to convert the floats to integers and run Morton ordering on those integers. For example consider the following C++ code. (I don't now how to do the same in C)

double x=-1.123456789123456E205;
int64_t i = reinterpret_cast<int64_t &>(x);
cout<<i<<endl;
output >>> i = -1548698869907112442

And ditto for the reaming x,y. Once I have the "reinterpreted" values ,I would like to use them as a subroutine for Morton coding.

I checked the above type cast and it worked fine in reverse

double y = reinterpret_cast<double &>(i);
cout<<setprecision(16)<<y<<endl;
output>>-1.123456789123456e+205

I managed to find some codes for Morton coding, even some of them on this forum, but non of them used int64_t in 3D. Hence I am going to need the help of the experts on the forum, how to encode and decode int64_t integers.

I managed to reverse engineer the following code. Unfortunately there is some bug, I am not getting the proper numbers when I run the decode part. I would appreciate any help to figure out what is wrong.

2D morton code encode/decode 64bits.

#include <iostream>
#include <stdint.h>
#include<iomanip>

using namespace std;

uint64_t code_2D_M(double xd,double yd){

uint64_t x = reinterpret_cast<uint64_t& >(xd);
uint64_t y = reinterpret_cast<uint64_t& >(yd);

x = (x | (x << 16)) & 0x0000FFFF0000FFFF;
x = (x | (x << 8)) & 0x00FF00FF00FF00FF;
x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0F;
x = (x | (x << 2)) & 0x3333333333333333;
x = (x | (x << 1)) & 0x5555555555555555;

y = (y | (y << 16)) & 0x0000FFFF0000FFFF;
y = (y | (y << 8)) & 0x00FF00FF00FF00FF;
y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0F;
y = (y | (y << 2)) & 0x3333333333333333;
y = (y | (y << 1)) & 0x5555555555555555;

return x | (y << 1);

}


uint64_t code_3D_M(double xd,double yd,double zd){


uint64_t x = reinterpret_cast<uint64_t& >(xd);
uint64_t y = reinterpret_cast<uint64_t& >(yd);
uint64_t z = reinterpret_cast<uint64_t& >(zd);

x = (x | (x << 16)) & 0x0000FFFF0000FFFF;
x = (x | (x << 8)) & 0x00FF00FF00FF00FF;
x = (x | (x << 4)) & 0x0F0F0F0F0F0F0F0F;
x = (x | (x << 2)) & 0x3333333333333333;
x = (x | (x << 1)) & 0x5555555555555555;

y = (y | (y << 16)) & 0x0000FFFF0000FFFF;
y = (y | (y << 8)) & 0x00FF00FF00FF00FF;
y = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0F;
y = (y | (y << 2)) & 0x3333333333333333;
y = (y | (y << 1)) & 0x5555555555555555;

z = (y | (y << 16)) & 0x0000FFFF0000FFFF;
z = (y | (y << 8)) & 0x00FF00FF00FF00FF;
z = (y | (y << 4)) & 0x0F0F0F0F0F0F0F0F;
z = (y | (y << 2)) & 0x3333333333333333;
z = (y | (y << 1)) & 0x5555555555555555;

return x | (y << 1) | (z << 2);

}

double decode_M(uint64_t x)
{
    x = x & 0x5555555555555555;
    x = (x | (x >> 1)) & 0x3333333333333333;
    x = (x | (x >> 2)) & 0x0F0F0F0F0F0F0F0F;
    x = (x | (x >> 4)) & 0x00FF00FF00FF00FF;
    x = (x | (x >> 8)) & 0x0000FFFF0000FFFF;
    x = (x | (x >> 16)) & 0xFFFFFFFFFFFFFFFF;
    return reinterpret_cast<double& >(x);
}


int main (void){

uint64_t  mort;
double  x,y,z;

// test input
x=2.123456789123459E205;
y=1.789789123456129E205;
z=9.999999912345779E205;

// echo the input
cout<<setprecision(17)<<x<<endl;
cout<<setprecision(17)<<y<<endl;
cout<<setprecision(17)<<z<<endl;

// encode 2D case
mort = code_2D_M(x,y);
//decode and print the results to see if all was fine
cout<<setprecision(17)<<decode_M(mort>>0)<<endl;
cout<<setprecision(17)<<decode_M(mort>>1)<<endl;

// encode 3D case
mort = code_3D_M(x,y,z);
//decode and print the results to see if all was fine
cout<<setprecision(17)<<decode_M(mort>>0)<<endl;
cout<<setprecision(17)<<decode_M(mort>>1)<<endl;
cout<<setprecision(17)<<decode_M(mort>>2)<<endl;



return 0;
}

I am doing this because I would like not storing the coordinates as a 3D point (x,y,z) but rather as a single long integer and decode them when needed. By doing so I will reduce the size of my coordinate storage array 3-fold.

" I presume that I can use type cast to convert the floats to integers and run Morton ordering on those integers" -- you *could*. I suppose you'll effectively be looking at the locality of exponent, mantissa, sign. Something makes me think that if your dynamic range is very large it may not give you the results you want. — Brian Cain, Oct 14 '15 at 17:29
If your dynamic range is small enough maybe you could just do a linear scale to get from `float` to `int64`. — Brian Cain, Oct 14 '15 at 17:30
I am basically trying to cram 3 numbers into one and then decode them back. I am doing this for memory saving purposes and I was thinking that Morton order might be the simplest solution. — Alexander Cska, Oct 14 '15 at 17:50
If you are attempting to do a whole lot of numbers, you might be better off using a more generic compressor on the whole data stream rather than a per triplet solution as you have here. — Michael Dorgan, Oct 14 '15 at 18:24
Well I need to compress 3 floats. I am not familiar with the alternative solutions. — Alexander Cska, Oct 14 '15 at 19:00
@Alexander Cska: Obviously you can't squeeze 192 bits of information into 64-bits of storage without loss of precision, and interleaving the floating point encoding like this will simply throw away the most important parts (and decode to denormalized noise). Effective compression will require knowledge of your data we lack to figure out where bits of precision can safely be omitted. Generally speaking though switching to 32-bit floats is probably a reasonable first attempt. For coordinates space fixed-point or shared exponents may be suitable, arithmetic invariants exploited, and so on — doynax, Oct 14 '15 at 19:57
Indeed, conversion to integers by `reinterpret_cast` produces problems with the storage. I am saving coordinates. The initial numbers are Fortran signed double precision, with 17 places after the coma and exponent range to 205. I use iso-c binding to call the c functions. This sort of precision is important because I am running a numerical simulation code and truncating the precision might result in divergence or machine dependent results. — Alexander Cska, Oct 14 '15 at 20:06
@Alexander Cska: Then I don't quite see the problem. As long as the down-conversion is carefully written and consistently performed there's no reason why it can't be deterministic. Beyond you'll have to try lossless compression, which is usually less than stellar at dealing with floating-point representations. Again, Morton codes simply rearrange the order of the bits. They don't somehow magically turn 192-bits of information into 64-bits of data. — doynax, Oct 14 '15 at 20:19
@doynax: Well in numerical analysis if you run single precision real arithmetic, you will get as many results as there are platforms your code is running on. Not to mention that in some cases, the lack of precision might lead to failed convergence. Therefore my (x,y,z) have to be double precision. My problem can be formulated as being how to encode 3 doubles into a unique integer. I think that due to the nature of my code this might be a lost cause. I can not say in advanced how my randomly generated (x,y,z) might look like. — Alexander Cska, Oct 14 '15 at 20:59
@Alexander Cska: You'll need a compiler (or suitable compiler switches) for strictly enforcing IEEE-754 arithmetic certainly, and only use well-defined basic operations. That's no different whether you're using single or double-precision though, and for a limited C function solely dedicated to conversion between the two it shouldn't be a problem. In any event you must have some knowledge of the data if you've determined that you can make do with only double precision, perhaps that can be taken further. — doynax, Oct 14 '15 at 21:18

Morton ordering in 3D usign unit64_t as input

0 Answers0