LSST Applications g0603fd7c41+022847dfd1,g0aad566f14+f45185db35,g180d380827+40e913b07a,g2079a07aa2+86d27d4dc4,g2305ad1205+696e5f3872,g2bbee38e9b+047b288a59,g337abbeb29+047b288a59,g33d1c0ed96+047b288a59,g3a166c0a6a+047b288a59,g3d1719c13e+f45185db35,g3de15ee5c7+5201731f0d,g487adcacf7+19f9b77d7d,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+248b16177b,g63cd9335cc+585e252eca,g858d7b2824+f45185db35,g88963caddf+0cb8e002cc,g991b906543+f45185db35,g99cad8db69+1747e75aa3,g9b9dfce982+78139cbddb,g9ddcbc5298+9a081db1e4,ga1e77700b3+a912195c07,gae0086650b+585e252eca,gb0e22166c9+60f28cb32d,gb3a676b8dc+b4feba26a1,gb4b16eec92+f82f04eb54,gba4ed39666+c2a2e4ac27,gbb8dafda3b+215b19b0ab,gc120e1dc64+b0284b5341,gc28159a63d+047b288a59,gc3e9b769f7+dcad4ace9a,gcf0d15dbbd+78139cbddb,gdaeeff99f8+f9a426f77a,ge79ae78c31+047b288a59,w.2024.19
LSST Data Management Base Package
|
This file contains functions for space-filling curves. More...
#include <cstdint>
#include <tuple>
Go to the source code of this file.
Namespaces | |
namespace | lsst |
namespace | lsst::sphgeom |
Functions | |
uint64_t | lsst::sphgeom::mortonIndex (uint32_t x, uint32_t y) |
mortonIndex interleaves the bits of x and y. | |
std::tuple< uint32_t, uint32_t > | lsst::sphgeom::mortonIndexInverse (uint64_t z) |
mortonIndexInverse separates the even and odd bits of z. | |
uint64_t | lsst::sphgeom::mortonToHilbert (uint64_t z, int m) |
mortonToHilbert converts the 2m-bit Morton index z to the corresponding Hilbert index. | |
uint64_t | lsst::sphgeom::hilbertToMorton (uint64_t h, int m) |
hilbertToMorton converts the 2m-bit Hilbert index h to the corresponding Morton index. | |
uint64_t | lsst::sphgeom::hilbertIndex (uint32_t x, uint32_t y, int m) |
hilbertIndex returns the index of (x, y) in a 2-D Hilbert curve. | |
std::tuple< uint32_t, uint32_t > | lsst::sphgeom::hilbertIndexInverse (uint64_t h, int m) |
hilbertIndexInverse returns the point (x, y) with Hilbert index h, where x and y are m bit integers. | |
uint8_t | lsst::sphgeom::log2 (uint64_t x) |
uint8_t | lsst::sphgeom::log2 (uint32_t x) |
This file contains functions for space-filling curves.
Mappings between 2-D points with non-negative integer coordinates and their corresponding Morton or Hilbert indexes are provided.
The Morton order implementation, mortonIndex, is straightforward. The Hilbert order implementation is derived from Algorithm 2 in:
C. Hamilton. Compact Hilbert indices. Technical Report CS-2006-07, Dalhousie University, Faculty of Computer Science, Jul 2006. https://www.cs.dal.ca/research/techreports/cs-2006-07
Using the variable names from that paper, n is fixed at 2. As a first step, the arithmetic in the loop over the bits of the input coordinates is replaced by a table lookup. In particular, the lookup maps the values of (e, d, l) at the beginning of a loop iteration to the values (e, d, w) at the end. Since e and d can both be represented by a single bit, and l and w are 2 bits wide, the lookup table has 16 4 bit entries and fits in a single 64 bit integer constant (0x8d3ec79a6b5021f4). The implementation then looks like:
inline uint64_t hilbertIndex(uint32_t x, uint32_t y, uint32_t m) { uint64_t const z = mortonIndex(x, y); uint64_t h = 0; uint64_t i = 0; for (m = 2 * m; m != 0;) { m -= 2; i = (i & 0xc) | ((z >> m) & 3); i = UINT64_C(0x8d3ec79a6b5021f4) >> (i * 4); h = (h << 2) | (i & 3); } return h; }
Note that interleaving x and y with mortonIndex beforehand allows the loop to extract 2 bits at a time from z, rather than extracting bits from x and y and then pasting them together. This lowers the total operation count.
Performance is further increased by executing j loop iterations at a time. This requires using a larger lookup table that maps the values of e and d at the beginning of a loop iteration, along with 2j input bits, to the values of e and d after j iterations, along with 2j output bits. In this implementation, j = 3, which corresponds to a 256 byte LUT. On recent Intel CPUs the LUT fits in 4 cache lines, and, because of adjacent cache line prefetch, should become cache resident after just 2 misses.
For a helpful presentation of the technical report, as well as a reference implementation of its algorithms in Python, see Pierre de Buyl's notebook. The Hilbert curve lookup tables below were generated by a modification of that code (available in makeHilbertLuts.py).
Definition in file curve.h.