Skip to content

Decoding Morton Codes

December 13, 2009

There’s lots of material on the web about computing Morton codes (also called Morton keys or Morton numbers) efficiently – bitwise interleaving of two or more numbers. This may sound esoteric, but it’s surprisingly useful in some applications. If you haven’t heard of Morton codes yet, step by Wikipedia or look into a book like Real-Time Collision Detection to learn more about them. Anyway, the subject of this post is not to introduce you to Morton codes, but rather to fill a rather curious gap: As I discovered a few months ago, there’s lots of material on the web and in books about how to generate morton codes from 2 or 3 numbers, but I didn’t find a single site explaining how to de-interleave the bits again to get the original numbers back from a morton code! I figured it’s time to change that.

The “classic” algorithms to generate Morton codes look like this:

uint32 EncodeMorton2(uint32 x, uint32 y)
{
  return (Part1By1(y) << 1) + Part1By1(x);
}

uint32 EncodeMorton3(uint32 x, uint32 y, uint32 z)
{
  return (Part1By2(z) << 2) + (Part1By2(y) << 1) + Part1By2(x);
}

// "Insert" a 0 bit after each of the 16 low bits of x
uint32 Part1By1(uint32 x)
{
  x &= 0x0000ffff;                  // x = ---- ---- ---- ---- fedc ba98 7654 3210
  x = (x ^ (x <<  8)) & 0x00ff00ff; // x = ---- ---- fedc ba98 ---- ---- 7654 3210
  x = (x ^ (x <<  4)) & 0x0f0f0f0f; // x = ---- fedc ---- ba98 ---- 7654 ---- 3210
  x = (x ^ (x <<  2)) & 0x33333333; // x = --fe --dc --ba --98 --76 --54 --32 --10
  x = (x ^ (x <<  1)) & 0x55555555; // x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
  return x;
}

// "Insert" two 0 bits after each of the 10 low bits of x
uint32 Part1By2(uint32 x)
{
  x &= 0x000003ff;                  // x = ---- ---- ---- ---- ---- --98 7654 3210
  x = (x ^ (x << 16)) & 0xff0000ff; // x = ---- --98 ---- ---- ---- ---- 7654 3210
  x = (x ^ (x <<  8)) & 0x0300f00f; // x = ---- --98 ---- ---- 7654 ---- ---- 3210
  x = (x ^ (x <<  4)) & 0x030c30c3; // x = ---- --98 ---- 76-- --54 ---- 32-- --10
  x = (x ^ (x <<  2)) & 0x09249249; // x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
  return x;
}

The meat is in the Part1By1 and Part1By2 functions, which separate the bits from each other. Reversing this function is actually very easy – it mainly boils down to executing the function in reverse, turning left shifts into right shifts and moving the masks aroundd a bit. All of this is pretty easy to work out by taking one single step and figuring out how to reverse it. So, without further ado, here’s the inverses of Part1By1 and Part1By2:

// Inverse of Part1By1 - "delete" all odd-indexed bits
uint32 Compact1By1(uint32 x)
{
  x &= 0x55555555;                  // x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
  x = (x ^ (x >>  1)) & 0x33333333; // x = --fe --dc --ba --98 --76 --54 --32 --10
  x = (x ^ (x >>  2)) & 0x0f0f0f0f; // x = ---- fedc ---- ba98 ---- 7654 ---- 3210
  x = (x ^ (x >>  4)) & 0x00ff00ff; // x = ---- ---- fedc ba98 ---- ---- 7654 3210
  x = (x ^ (x >>  8)) & 0x0000ffff; // x = ---- ---- ---- ---- fedc ba98 7654 3210
  return x;
}

// Inverse of Part1By2 - "delete" all bits not at positions divisible by 3
uint32 Compact1By2(uint32 x)
{
  x &= 0x09249249;                  // x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
  x = (x ^ (x >>  2)) & 0x030c30c3; // x = ---- --98 ---- 76-- --54 ---- 32-- --10
  x = (x ^ (x >>  4)) & 0x0300f00f; // x = ---- --98 ---- ---- 7654 ---- ---- 3210
  x = (x ^ (x >>  8)) & 0xff0000ff; // x = ---- --98 ---- ---- ---- ---- 7654 3210
  x = (x ^ (x >> 16)) & 0x000003ff; // x = ---- ---- ---- ---- ---- --98 7654 3210
  return x;
}

Using these, getting the original x, y and z coordinates from Morton codes is then trivial:

uint32 DecodeMorton2X(uint32 code)
{
  return Compact1By1(code >> 0);
}

uint32 DecodeMorton2Y(uint32 code)
{
  return Compact1By1(code >> 1);
}

uint32 DecodeMorton3X(uint32 code)
{
  return Compact1By2(code >> 0);
}

uint32 DecodeMorton3Y(uint32 code)
{
  return Compact1By2(code >> 1);
}

uint32 DecodeMorton3Z(uint32 code)
{
  return Compact1By2(code >> 2);
}

There you go, hope it comes in handy.

From → Coding

21 Comments
  1. nikitos permalink

    Hi, there is a wonderful resource on such codes, and many others:
    http://www.jjj.de/fxt/fxtpage.html, look at bitwizardy section in book.

  2. Do I’m confused or what is this Morton-Code not a Z-order leads to Z-curve, which is the Lebesgue space filling curve to subdivide the plane.

  3. Correct. It’s known under a lot of different names.

  4. Thank you for this statement!

  5. Joerg permalink

    Saying that I used 16-entry-LUTs for the nibbles is a bit lame I suppose…thanks for sharing that code.

  6. Works perfectly. Thanks!

  7. John permalink

    What would this look like for 64 bit numbers? (or 128bit.. etc) I can’t quite see how to extend the code for those types.. i think I’ve got the Part1By1 and Compact1By1 bits right, but it’s the 1By2 ones I’m confused by.

    uint32 Part1By1( uint64 X )
    {
    X &= 0x00000000ffffffff;
    X = ( X ^ (X << 16 ) ) & 0x0000ffff0000ffff;
    X = ( X ^ (X << 8 ) ) & 0x00ff00ff00ff00ff;
    X = ( X ^ (X << 4 ) ) & 0x0f0f0f0f0f0f0f0f;
    X = ( X ^ (X << 2 ) ) & 0x3333333333333333;
    X = ( X ^ (X <> 1 ) ) & 0x3333333333333333;
    X = ( X ^ ( X >> 2 ) ) & 0x0f0f0f0f0f0f0f0f;
    X = ( X ^ ( X >> 4 ) ) & 0x00ff00ff00ff00ff;
    X = ( X ^ ( X >> 8 ) ) & 0x0000ffff0000ffff;
    X = ( X ^ ( X >> 16 ) ) & 0x00000000ffffffff;
    return X;
    }

    How should the 1By2 variants looks?

    • Yeah, the 1:1 interleaves have regular enough structure to just type it in, but 1:2 requires some more work. I’d generally suggest writing a small program to generate the bitmasks :). Though your code looks wrong – it should end after the X ^ (X << 2)... with X = (X ^ (X << 1)) & 0x555555...; (16 ‘5’s).

      The basic idea used here is this: For every input bit, you keep track of the distance it has to move. For the 1:2 ratio interleaves, bit N (I number bits in increasing order starting from the LSB) has to travel 2*N positions to the left.

      You then build the moves out of successive power-of-2 distance shifts (effectively looking at the individual bits for the movement distances). For example, in my original Part1By2, the first line is this:
      x = (x ^ (x << 16)) & 0xff0000ff;
      but it might as well be
      x = (x ^ (x << 16)) & 0x030000ff;
      since we masked out the high-level bits of x in the line before.

      So where does that mask come from? Only bits 8 and 9 in the original value need to travel a distance >=16 positions. So what that mask really is is (0x300 << 16) | 0xff) (keep the shifted-by-16 versions of bits 8 and 9, stay with the unshifted versions for the remaining bits).

      Here’s the total list of move distances:

      bit 0 - distance  0 = 0b00000
      bit 1 - distance  2 = 0b00010
      bit 2 - distance  4 = 0b00100
      bit 3 - distance  6 = 0b00110
      bit 4 - distance  8 = 0b01000
      bit 5 - distance 10 = 0b01010
      bit 6 - distance 12 = 0b01100
      bit 7 - distance 14 = 0b01110
      bit 8 - distance 16 = 0b10000
      bit 9 - distance 18 = 0b10010
      

      So in the second pass, we’re going to have to move source bits 4, 5, 6, and 7 (the ones with bit 3 set in their move distance), while the rest remain where they are (note you need to compute the masks from their positions after the first pass!). Thus the mask is (0xf0 << 8) | 0x0300000f == 0x0300f00f. In the third pass, we’re going to move source bits 2, 3, 6 and 7 (which have bit 2 set in their move distance). Again, we need to work from the current position of the bits, and we get (0xc00c << 4) | 0x03003003 == 0x030c30c3.(Note that figuring out these masks manually is much easier if you write the positions of the bits after each line, as I do in the code – you can just read it off from there, basically). And so on. Does that help? As said, I wrote a small program to do this for me; it’s not hard to do manually, just annoying and somewhat error-prone.

      This process (moving by successive power-of-2 distances) will always work and give a result with <=log2(N) passes for N-bit output values, but it won't always give you minimum number of passes. For example, if you want to do a 1:3 interleave with a sequence of zeros, you want your move distances to be of the form 3*2^n, not 2^n (the latter will give you more passes than necessary). The xor trick (actually plain or would be fine too) only works if you move a bit to a position that was unoccupied before. In the general case you need to use something like (x & bits_that_stay) | ((x & ~bits_that_stay) << distance) – a bit more work for a pass.

      • John permalink

        Thanks for the tips, and somehow I’d pasted code for Part1By1 and Compact1By1 into one function. not sure how i’d managed that!

        Thanks again,

  8. Gabriel permalink

    I posted a stackoverflow addition with your inputs and hints! posted also the python script!
    http://stackoverflow.com/a/18528775/293195

  9. Just to let you know that in July 2016, almost 7 years from your original post, the information you presented on Morton numbers is still the best around. I originally did an 8 bit Morton encoding/decoding using a table lookup. Done on an ATmega328. I put the table in program memory to save dynamic memory but the lookup was pretty slow. I found your code and converted my program to compute the Morton and inverse Morton using your method. That function then ran about 30 times faster.

  10. Okey, sorry for the noise Ryg, I’m trying to write some code but HTML entities get broken. This is my last attempt:

    I used this other methoid in 2007 for a 4k intro, wrote some tiny bit about it here: http://www.iquilezles.org/www/articles/wavelet/wavelet.htm

    ivec2 imorton( int i, int level)
    {
        ivec2 res = ivec2(0,0);
        for( int k=0; k < level; k++ )
        {
            res += (ivec2(i,i>>1) & 1) << k;
            i >>= 2;
        }
        return res;
    }
  11. I’m wondering how I can deal wit not square 2d grid? when performing decode can I just discard points those are outside or is there a better way?

    • You round up both dimensions to the nearest power of 2 and just pack the major axis linearly. (I.e. stop interleaving bits once you run out of bits in the minor axis.)

  12. Hi, shouldn’t the bit shifting be done in number of bits not number of hex digits?

    E.g.
    x = (x ^ (x << 8)) & 0x00ff00ff; // x = —- —- fedc ba98 —- —- 7654 3210

    shouldn't it be

    x = (x ^ (x << 24)) & 0x00ff00ff;

    since 8 places left is 8×4=24 bits

    • The things in the comments are not hex digits! They’re bits. Which is why there’s 32 of them for a 32-bit number.

      It’s not showing an example number. It’s showing, for each position in the current 32-bit word, which bit in the original word it came from.

Trackbacks & Pingbacks

  1. Out-Of-Core construction of Sparse Voxel Octrees
  2. Morton Codes | Asger Hoedt
  3. Morton encoding/decoding through bit interleaving: Implementations
  4. Implementing Morton ordering for chunked voxel data - Volumes Of Fun

Leave a reply to kannaiah Cancel reply