LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Namespaces | Functions
curve.h File Reference

This file contains functions for space-filling curves. More...

#include <cstdint>
#include <tuple>

Go to the source code of this file.

Namespaces

 lsst
 A base class for image defects.
 
 lsst::sphgeom
 

Functions

uint64_t lsst::sphgeom::mortonIndex (uint32_t x, uint32_t y)
 mortonIndex interleaves the bits of x and y. More...
 
std::tuple< uint32_t, uint32_t > lsst::sphgeom::mortonIndexInverse (uint64_t z)
 mortonIndexInverse separates the even and odd bits of z. More...
 
uint64_t lsst::sphgeom::mortonToHilbert (uint64_t z, int m)
 mortonToHilbert converts the 2m-bit Morton index z to the corresponding Hilbert index. More...
 
uint64_t lsst::sphgeom::hilbertToMorton (uint64_t h, int m)
 hilbertToMorton converts the 2m-bit Hilbert index h to the corresponding Morton index. More...
 
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. More...
 
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. More...
 
uint8_t lsst::sphgeom::log2 (uint64_t x)
 
uint8_t lsst::sphgeom::log2 (uint32_t x)
 

Detailed Description

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.