gpt4_cache_locality_optimize_hilbert.md

QUESTION: Using C++ 20 and glm math library only, implement sorting of a set of 3D points using the complete Hilbert space filling curve implementation, including 3D rotation and reflection. For simplicity sake, all coordinates of the 3D points are positive, but each point can have arbitrary floating point value (not just integers), but the bounding box containing all points is known beforehand.

GPT-4 ANSWER:

Here is an implementation that uses a classic 3D Hilbert curve mapping to sort points. Note that this does handle floating point coordinates but the entire space is discretized to the closest integer to compute the Hilbert value.

This example uses the boost library (Precompiled library available for Visual Studio in vcpkg).

#include <iostream>
#include <vector>
#include <algorithm>
#include <glm/glm.hpp>
#include <limits>

unsigned long long int hilbert_index(int x, int y, int z, unsigned order)
{
    unsigned long long int h = 0; 
    for (int i = 0; i < order; i++) {
        unsigned int xi = (x >> i) & 1;
        unsigned int yi = (y >> i) & 1;
        unsigned int zi = (z >> i) & 1;
        h |= (xi ^ yi ^ zi) << (2 * i);
        h |= ((yi ^ zi) << (2 * i)) | (1 << (2 * i));
    }
    return h;
}

template<class T> 
void hilbertSort(std::vector<glm::tvec3<T>>& points, const glm::tvec3<T>& boxMin, const glm::tvec3<T>& boxMax, unsigned order=10) 
{
    T invGridSize[3] = { order / (boxMax.x - boxMin.x), 
                         order / (boxMax.y - boxMin.y),
                         order / (boxMax.z - boxMin.z) };

    std::sort(points.begin(), points.end(), [&](const glm::tvec3<T>& left, const glm::tvec3<T>& right) {
        int leftDiscrete[3] = { (int)((left.x - boxMin.x) * invGridSize[0]),
                                (int)((left.y - boxMin.y) * invGridSize[1]),
                                (int)((left.z - boxMin.z) * invGridSize[2]) };
        int rightDiscrete[3] = { (int)((right.x - boxMin.x) * invGridSize[0]),
                                 (int)((right.y - boxMin.y) * invGridSize[1]),
                                 (int)((right.z - boxMin.z) * invGridSize[2]) };

        return hilbert_index(leftDiscrete[0], leftDiscrete[1], leftDiscrete[2], order)
             < hilbert_index(rightDiscrete[0], rightDiscrete[1], rightDiscrete[2], order);
    });
}

int main()
{
    std::vector<glm::vec3> points = {... };
    glm::vec3 boxMin = {... };
    glm::vec3 boxMax = {... };
    hilbertSort(points, boxMin, boxMax);
    return 0;
}

The hilbertSort function takes an array of glm::tvec3 points and the bounding box of the point cloud. The order parameter is the number of levels of recursion for the Hilbert curve, which affects how finely the space is divided. For sparse datasets, a lower order may be sufficient to speed up the sorting process.

The invGridSize array is computed from the bounding box size and the grid order. Before the Hilbert index is computed, the coordinates (relative to the bounding box) are scaled to fit the grid and converted to integers.

The Hilbert curve has the property that it preserves neighborhood relationships; that is, points that are near each other in 3D space are also likely to be near each other when ordered by their Hilbert index. Therefore, if a certain region of space is denser with points, these points will be clustered together in the sorted list. This can be beneficial for certain algorithms.

Please note that this implementation assumes that the point cloud fits into a cube. If the point cloud is noticeably rectangular (i.e. boxMin and boxMax differ significantly in one or two dimensions), it may be better to use a different space-filling curve or to transform the point cloud to "stretch" it into a cube before sorting.