Morton numbers

Long time no posting, but I have excuses (also I’m posting some at openstreetmap user diaries).

So anyway here’s a cheap trick I came up with but which you might already know.  If you’re indexing any georeferenced data, such as when doing any fun stuff with OpenStreetMap data, you’ve probably wanted to index by location among other things, and location is two or three dimensional (without loss of generality assume two as in GIS).  So obviously you can combine latitude and longitude as one key and index by that but that’s only good for searching for exact pairs of values.  If your index is for a hash table then you can’t hope for anything more but if it’s for sorting of an array you can do a little better (well, here’s my trick): convert the two numbers to fixed point and interleave their bits to make one number.  This is better because two positions that are close to each other in an array sorted by this number probably are close to each other on the map.  You could probably use floating point too if you stuff the exponent in the most significant bits and get a result similar to some degree.  With fixed point you can then compare only the top couple of bits when searching in the array to locate something with a desired accuracy.

Converting to and from the interleaved bits form is straight forward and you can easily come up with a O(log(number of bits)) procedure (5 steps for 32 bit lat / lon) or use lookup tables as suggested by the Bit Twiddling Hacks page, where I learnt they’re called Morton numbers.  32-bit lat/lon will give you a 64-bit number and that should be accurate enough for most uses if you map the whole -90 – 90 / -180 – 180 deg range to integers.  Even 20-bit lat/lon (5 bytes for the index) gives you 0.0003 deg accuracy.

So what else can you do with this notation? Obviously you can compare two numbers and use bisection search in arrays or the different kinds of trees.  You  can not add or subtract them directly (or rather, you won’t get useful results) but you can add / subtract individual coordinates without converting to normal notation and back, here’s how:

First separate latitude from longitude by masking:

uint64_t x = a & 0x5555555555555555;
uint64_t y = a & 0xaaaaaaaaaaaaaaaa;

Now you can subtract two numbers directly, you’ll notice that the carry flags are correctly carried over the unused bits, you’ll just need to mask them out of the result:

uint64_t difference(uint64_t a, uint64_t b)
	return ((ax - bx) & 0x5555555555555555) |
		((ay - by) & 0xaaaaaaaaaaaaaaaa);

(you can also use this Masked Merge hack from the page I linked earlier).

The result is signed two’s complement with two sign bits in the top bits.

Now something much less obvious is that if you want to calculate absolute difference, you can call abs() directly on the result of subtraction and only mask out the unused bits afterwards.  How does this work?  The top bit in (ax - bx) always equals the sign bit even if ax and bx only use even bits (top bit is odd), so this part is ok.  Now, if the number is positive then there’s nothing to do with it.  If it’s negative, then abs negates it again (strips the minus).  Conveniently -x equals ~(x - 1) in two’s complement, so let’s see what these two operations do to a negative (ax - bx)~ or bitwise negation just works because it inverts all bits including the ones we’re interested in.  The x – 1 part also works because it flips all the bits until the first 1 bit starting from lowest bit, and you’ll find, although it may be tricky to see, that the first bit set in (ax - bx) is always even (or always odd).

uint64_t distance(uint64_t a, uint64_t b)
	return (llabs(ax - bx) & 0x5555555555555555) |
		(llabs(ay - by) & 0xaaaaaaaaaaaaaaaa);

Addition requires a little trick for the carry flags to work: just set all unused bits in either ax or bx:

uint64_t sum(uint64_t a, uint64_t b)
	uint64_t ax = a & 0x5555555555555555;
	uint64_t ay = a & 0xaaaaaaaaaaaaaaaa;
	uint64_t bx = b | 0xaaaaaaaaaaaaaaaa;
	uint64_t by = b | 0x5555555555555555;

	return ((ax + bx) & 0x5555555555555555) |
		((ay + by) & 0xaaaaaaaaaaaaaaaa);

5 Responses to “Morton numbers”

  1. ephemient Says:

    IEEE-754 floating point numbers are designed such that they can be compared/ordered simply by treating their binary representations as 2’s-complement integers.

    So you don’t need to switch to fixed-point to get “two positions that are close to each other in an array sorted by this number probably are close to each other on the map”; it works with even when mixing floats. Looking at the high bits to “locate something with a desired accuracy” does have a pretty warped scale, though, and the arithmetic tricks here won’t work either.

  2. Eric Says:

    I have a webapp that uses Google Maps and a database with lots of sites. I wanted a way of getting lat-lon into some 1D order. I hit upon the idea of scaling the lat & lon to 16-bit integers and interleaving the bits. I should have known I wasn’t the only one! Today I googled and found out that my idea is quite widely used.

    I used (lon+180)*180 and (lat+90)*180 to get 16 & 15 bit integers. I’m using JavaScript and PHP/MySQL.

    So I would use this number to filter which sites I load from my DB. I get the two opposite corners, convert both to the integer to see what range I need. A small window will always give you a small integer range unless your window happens to cross a major division.

    Alas, that’s exactly what happened! Adelaide (my city) is 138.6E.
    (138.6 + 180) * 180 = 0xE004. So the 7/8 of the way around the world longitude line ran right through the area I was looking at! This rendered the filter useless, and I’ve since refined it to get around such problems, but it was just funny that I chose an easy trick with a scaling system that would work almost anywhere but here!

    btw Morton number in SQL: conv(bin(n1),4,10) + 2*conv(bin(n2),4,10)

    • balrog Says:

      Hi Eric,
      [plug]I can’t help asking if you’ve considered switching over to OpenStreetMap for your web app?

      Re looking up sites inside a rectangular region, I pondered the same issue and my solution (although not implemented yet) is for every pair of opposite corners to be converted to four (or fewer) integer ranges mapping to one, two or four db queries, depending on where exactly the two corners fall in that “grid”. Complexity wise making it four lookups instead of one makes no difference. So if the corners were at, say (-1,-1), (3, 3), then you’d split it up into four square lookups:

      (-1, -1)-(0, 0),
      (-4, 0)-(0, 4),
      (0, -4)-(4, 0),
      (0, 0)-(4,4),

      each power-of-two-sized, then, in my case, I will just have the server return all results that matched and have the javascript filter out unwanted results. Assuming the original query is usually for objects visible on the screen, the query will almost always have a screen like width:height ratio, so nearly square, and if it’s square, the area returned from the four queries will never sum up to more than twice the original area.

      Note that you can do (lat+90)*360 in the conversion to integers to get a better resolution, unless you have other use for that 16th bit there. Cheers

  3. Eric Says:

    Mappage (my webapp) is here:

    The features I described are not yet on the live version though.

    OSM looks OK, and I’m generally a fan of free/open stuff. Adelaide is mostly covered now, and I suspect an old friend of mine has mapped most of the northern suburbs. GM was the best option when I started Mappage: for things like geocoding and translucent polygon overlay there was nothing else. Later it might be possible to decouple from GM so it can sit on something else, but unlikely.

    My solution to getting the rectangular range was much like what you suggested.

    Scale-wise I kept the scale the same in both directions, keeping the 16th bit clear so I could do comparisons without worrying about sign. Accuracy isn’t important as it’s only for rough sorting, and being in Australia the graticules are close to square.

  4. Reid Cantero Says:

    Sick of obtaining low amounts of useless traffic for your website? Well i wish to let you know about a new underground tactic which makes myself $900 on a daily basis on 100% AUTOPILOT. I could be here all day and going into detail but why dont you merely check their site out? There is a great video that explains everything. So if your serious about making easy cash this is the website for you. Auto Traffic Avalanche

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: