[orx-shapes] Add Hilbert and Morton point ordering in 2d and 3d

orx-shapes/src/commonMain/kotlin/Shapes.kt

@@ -0,0 +1 @@
+package org.openrndr.extra.shapes

orx-shapes/src/commonMain/kotlin/ordering/Hilbert2d.kt

@@ -0,0 +1,26 @@
+package org.openrndr.extra.shapes.ordering
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert2dDecode16Bit(hilbert: UInt): UIntArray {
+    val morton = hilbertToMorton2d(hilbert, 10)
+    return morton2dDecode16Bit(morton)
+}
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert2dEncode16Bit(index1: UInt, index2: UInt): UInt {
+    val morton = morton2dEncode16Bit(index1, index2)
+    return mortonToHilbert2d(morton, 16)
+}
+
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert2dDecode5Bit(hilbert: UInt): UIntArray {
+    val morton = hilbertToMorton2d(hilbert, 5)
+    return morton2dDecode5Bit(morton)
+}
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert2dEncode5Bit(index1: UInt, index2: UInt): UInt {
+    val morton = morton2dEncode5Bit(index1, index2)
+    return mortonToHilbert3d(morton, 5)
+}

orx-shapes/src/commonMain/kotlin/ordering/Hilbert3d.kt

@@ -0,0 +1,26 @@
+package org.openrndr.extra.shapes.ordering
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert3dDecode10Bit(hilbert: UInt): UIntArray {
+    val morton = hilbertToMorton3d(hilbert, 10)
+    return morton3dDecode10Bit(morton)
+}
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert3dEncode10Bit(index1: UInt, index2: UInt, index3: UInt): UInt {
+    val morton = morton3dEncode10Bit(index1, index2, index3)
+    return mortonToHilbert3d(morton, 10)
+}
+
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert3dDecode5Bit(hilbert: UInt): UIntArray {
+    val morton = hilbertToMorton3d(hilbert, 5)
+    return morton3dDecode5Bit(morton)
+}
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbert3dEncode5Bit(index1: UInt, index2: UInt, index3: UInt): UInt {
+    val morton = morton3dEncode5Bit(index1, index2, index3)
+    return mortonToHilbert3d(morton, 5)
+}

orx-shapes/src/commonMain/kotlin/ordering/ListVector2Extensions.kt

@@ -0,0 +1,75 @@
+package org.openrndr.extra.shapes.ordering
+
+import org.openrndr.math.Vector2
+import org.openrndr.math.map
+import org.openrndr.shape.Rectangle
+import org.openrndr.shape.bounds
+import kotlin.math.max
+import kotlin.math.pow
+
+enum class Axis2DPermutation {
+    XY,
+    YX,
+}
+
+fun List<Vector2>.mortonOrder(
+    scale: Double = 1.0,
+    permutation: Axis2DPermutation = Axis2DPermutation.XY,
+    bits: Int = 16,
+): List<Vector2> {
+    val bounds = this.bounds
+    val md = max(bounds.width, bounds.height) * scale
+    val rbounds = Rectangle(bounds.corner.x, bounds.corner.y, md, md)
+    val extend = 2.0.pow(bits.toDouble()) - 1.0
+    val inputPoints = map {
+        it.map(
+            rbounds.position(0.0, 0.0),
+            rbounds.position(1.0, 1.0),
+            Vector2(0.0, 0.0),
+            Vector2(extend, extend)
+        )
+    }
+    val mortonCodes = when (bits) {
+        5 -> when (permutation) {
+            Axis2DPermutation.XY -> inputPoints.map { morton2dEncode5Bit(it.x.toUInt(), it.y.toUInt()) }
+            Axis2DPermutation.YX -> inputPoints.map { morton2dEncode5Bit(it.y.toUInt(), it.x.toUInt()) }
+        }
+        16 -> when (permutation) {
+            Axis2DPermutation.XY -> inputPoints.map { morton2dEncode16Bit(it.x.toUInt(), it.y.toUInt()) }
+            Axis2DPermutation.YX -> inputPoints.map { morton2dEncode16Bit(it.y.toUInt(), it.x.toUInt()) }
+        }
+        else -> error("Only 5 and 16 bit modes are supported.")
+    }
+    return (this zip mortonCodes).sortedBy { it.second }.map { it.first }
+}
+
+fun List<Vector2>.hilbertOrder(
+    scale: Double = 1.0,
+    permutation: Axis2DPermutation = Axis2DPermutation.XY,
+    bits: Int = 16,
+): List<Vector2> {
+    val bounds = this.bounds
+    val md = max(bounds.width, bounds.height) * scale
+    val rbounds = Rectangle(bounds.corner.x, bounds.corner.y, md, md)
+    val extend = 2.0.pow(bits.toDouble()) - 1.0
+    val inputPoints = map {
+        it.map(
+            rbounds.position(0.0, 0.0),
+            rbounds.position(1.0, 1.0),
+            Vector2(0.0, 0.0),
+            Vector2(extend, extend)
+        )
+    }
+    val hilbertCodes = when (bits) {
+        5 -> when (permutation) {
+            Axis2DPermutation.XY -> inputPoints.map { hilbert2dEncode5Bit(it.x.toUInt(), it.y.toUInt()) }
+            Axis2DPermutation.YX -> inputPoints.map { hilbert2dEncode5Bit(it.y.toUInt(), it.x.toUInt()) }
+        }
+        16 -> when (permutation) {
+            Axis2DPermutation.XY -> inputPoints.map { hilbert2dEncode16Bit(it.x.toUInt(), it.y.toUInt()) }
+            Axis2DPermutation.YX -> inputPoints.map { hilbert2dEncode16Bit(it.y.toUInt(), it.x.toUInt()) }
+        }
+        else -> error("Only 5 and 16 bit modes are supported.")
+    }
+    return (this zip hilbertCodes).sortedBy { it.second }.map { it.first }
+}

orx-shapes/src/commonMain/kotlin/ordering/ListVector3Extensions.kt

@@ -0,0 +1,92 @@
+package org.openrndr.extra.shapes.ordering
+
+import org.openrndr.math.Vector3
+import org.openrndr.math.map
+import org.openrndr.shape.Box
+import org.openrndr.shape.bounds
+import kotlin.math.max
+
+enum class Axis3DPermutation {
+    XYZ,
+    XZY,
+    YXZ,
+    YZX,
+    ZXY,
+    ZYX
+}
+
+fun List<Vector3>.mortonOrder(
+    scale: Double = 1.0,
+    permutation: Axis3DPermutation = Axis3DPermutation.XYZ,
+    bits: Int = 10,
+): List<Vector3> {
+    val bounds = this.bounds
+    val md = max(max(bounds.width, bounds.height), bounds.depth) * scale
+    val rbounds = Box(bounds.corner.x, bounds.corner.y, bounds.corner.z, md, md, md)
+    val inputPoints = map {
+        it.map(
+            rbounds.position(0.0, 0.0, 0.0),
+            rbounds.position(1.0, 1.0, 1.0),
+            Vector3(0.0, 0.0, 0.0),
+            Vector3(1023.0, 1023.0, 1023.0)
+        )
+    }
+    val mortonCodes = when (bits) {
+        5 -> when (permutation) {
+            Axis3DPermutation.XYZ -> inputPoints.map { morton3dEncode5Bit(it.x.toUInt(), it.y.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.XZY -> inputPoints.map { morton3dEncode5Bit(it.x.toUInt(), it.z.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.YXZ -> inputPoints.map { morton3dEncode5Bit(it.y.toUInt(), it.x.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.YZX -> inputPoints.map { morton3dEncode5Bit(it.y.toUInt(), it.z.toUInt(), it.x.toUInt()) }
+            Axis3DPermutation.ZXY -> inputPoints.map { morton3dEncode5Bit(it.z.toUInt(), it.x.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.ZYX -> inputPoints.map { morton3dEncode5Bit(it.z.toUInt(), it.y.toUInt(), it.x.toUInt()) }
+        }
+        10 -> when (permutation) {
+            Axis3DPermutation.XYZ -> inputPoints.map { morton3dEncode10Bit(it.x.toUInt(), it.y.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.XZY -> inputPoints.map { morton3dEncode10Bit(it.x.toUInt(), it.z.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.YXZ -> inputPoints.map { morton3dEncode10Bit(it.y.toUInt(), it.x.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.YZX -> inputPoints.map { morton3dEncode10Bit(it.y.toUInt(), it.z.toUInt(), it.x.toUInt()) }
+            Axis3DPermutation.ZXY -> inputPoints.map { morton3dEncode10Bit(it.z.toUInt(), it.x.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.ZYX -> inputPoints.map { morton3dEncode10Bit(it.z.toUInt(), it.y.toUInt(), it.x.toUInt()) }
+        }
+        else -> error("Only 5 and 10 bit modes are supported.")
+    }
+    return (this zip mortonCodes).sortedBy { it.second }.map { it.first }
+}
+
+fun List<Vector3>.hilbertOrder(
+    scale: Double = 1.0,
+    permutation: Axis3DPermutation = Axis3DPermutation.XYZ,
+    bits: Int
+): List<Vector3> {
+    val bounds = this.bounds
+    val md = max(max(bounds.width, bounds.height), bounds.depth) * scale
+    val rbounds = Box(bounds.corner.x, bounds.corner.y, bounds.corner.z, md, md, md)
+    val inputPoints = map {
+        it.map(
+            rbounds.position(0.0, 0.0, 0.0),
+            rbounds.position(1.0, 1.0, 1.0),
+            Vector3(0.0, 0.0, 0.0),
+            Vector3(1023.0, 1023.0, 1023.0)
+        )
+    }
+    val hilbertCodes = when(bits) {
+        5 -> when (permutation) {
+            Axis3DPermutation.XYZ -> inputPoints.map { hilbert3dEncode5Bit(it.x.toUInt(), it.y.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.XZY -> inputPoints.map { hilbert3dEncode5Bit(it.x.toUInt(), it.z.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.YXZ -> inputPoints.map { hilbert3dEncode5Bit(it.y.toUInt(), it.x.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.YZX -> inputPoints.map { hilbert3dEncode5Bit(it.y.toUInt(), it.z.toUInt(), it.x.toUInt()) }
+            Axis3DPermutation.ZXY -> inputPoints.map { hilbert3dEncode5Bit(it.z.toUInt(), it.x.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.ZYX -> inputPoints.map { hilbert3dEncode5Bit(it.z.toUInt(), it.y.toUInt(), it.x.toUInt()) }
+        }
+        10 -> when (permutation) {
+            Axis3DPermutation.XYZ -> inputPoints.map { hilbert3dEncode10Bit(it.x.toUInt(), it.y.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.XZY -> inputPoints.map { hilbert3dEncode10Bit(it.x.toUInt(), it.z.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.YXZ -> inputPoints.map { hilbert3dEncode10Bit(it.y.toUInt(), it.x.toUInt(), it.z.toUInt()) }
+            Axis3DPermutation.YZX -> inputPoints.map { hilbert3dEncode10Bit(it.y.toUInt(), it.z.toUInt(), it.x.toUInt()) }
+            Axis3DPermutation.ZXY -> inputPoints.map { hilbert3dEncode10Bit(it.z.toUInt(), it.x.toUInt(), it.y.toUInt()) }
+            Axis3DPermutation.ZYX -> inputPoints.map { hilbert3dEncode10Bit(it.z.toUInt(), it.y.toUInt(), it.x.toUInt()) }
+        }
+        else -> error("Only 5 and 10 bit modes are supported.")
+    }
+    return (this zip hilbertCodes).sortedBy { it.second }.map { it.first }
+}

orx-shapes/src/commonMain/kotlin/ordering/Morton2d.kt

@@ -0,0 +1,88 @@
+@file:OptIn(ExperimentalUnsignedTypes::class)
+
+package org.openrndr.extra.shapes.ordering
+
+fun morton2dEncode5Bit(index1: UInt, index2: UInt): UInt {
+    // pack 2 5-bit indices into a 10-bit Morton code
+    var index1: UInt = index1
+    var index2: UInt = index2
+    index1 = index1 and 0x0000001fu
+    index2 = index2 and 0x0000001fu
+    index1 *= 0x01041041u
+    index2 *= 0x01041041u
+    index1 = index1 and 0x10204081u
+    index2 = index2 and 0x10204081u
+    index1 *= 0x00108421u
+    index2 *= 0x00108421u
+    index1 = index1 and 0x15500000u
+    index2 = index2 and 0x15500000u
+    return ((index1 shr 20) or (index2 shr 19))
+}
+
+fun morton2dDecode5Bit(morton: UInt): UIntArray { // unpack 2 5-bit indices from a 10-bit Morton code
+    var value1 = morton;
+    var value2 = (value1 shr 1)
+    value1 = value1 and 0x00000155u
+    value2 = value2 and 0x00000155u
+    value1 = value1 or (value1 shr 1)
+    value2 = value2 or (value2 shr 1)
+    value1 = value1 and 0x00000133u
+    value2 = value2 and 0x00000133u
+    value1 = value1 or (value1 shr 2)
+    value2 = value2 or (value2 shr 2)
+    value1 = value1 and 0x0000010fu
+    value2 = value2 and 0x0000010fu
+    value1 = value1 or (value1 shr 4)
+    value2 = value2 or (value2 shr 4)
+    value1 = value1 and 0x0000001fu
+    value2 = value2 and 0x0000001fu
+    return uintArrayOf(value1, value2)
+}
+
+fun morton2dEncode16Bit(index1: UInt, index2: UInt): UInt { // pack 2 16-bit indices into a 32-bit Morton code
+    var index1: UInt = index1
+    var index2: UInt = index2
+    index1 = index1 and 0x0000ffffu
+    index2 = index2 and 0x0000ffffu
+    index1 = index1 or (index1 shl 8)
+    index2 = index2 or (index2 shl 8)
+    index1 = index1 and 0x00ff00ffu
+    index2 = index2 and 0x00ff00ffu
+    index1 = index1 or (index1 shl 4)
+    index2 = index2 or (index2 shl 4)
+    index1 = index1 and 0x0f0f0f0fu
+    index2 = index2 and 0x0f0f0f0fu
+    index1 = index1 or (index1 shl 2)
+    index2 = index2 or (index2 shl 2)
+    index1 = index1 and 0x33333333u
+    index2 = index2 and 0x33333333u
+    index1 = index1 or (index1 shl 1)
+    index2 = index2 or (index2 shl 1)
+    index1 = index1 and 0x55555555u
+    index2 = index2 and 0x55555555u
+    return (index1 or (index2 shl 1))
+}
+
+fun morton2dDecode16Bit(morton: UInt): UIntArray { // unpack 2 16-bit indices from a 32-bit Morton code
+    var value1 = morton;
+    var value2 = (value1 shr 1);
+    value1 = value1 and 0x55555555u
+    value2 = value2 and 0x55555555u
+    value1 = value1 or (value1 shr 1)
+    value2 = value2 or (value2 shr 1)
+    value1 = value1 and 0x33333333u
+    value2 = value2 and 0x33333333u
+    value1 = value1 or (value1 shr 2)
+    value2 = value2 or (value2 shr 2)
+    value1 = value1 and 0x0f0f0f0fu
+    value2 = value2 and 0x0f0f0f0fu
+    value1 = value1 or (value1 shr 4)
+    value2 = value2 or (value2 shr 4)
+    value1 = value1 and 0x00ff00ffu
+    value2 = value2 and 0x00ff00ffu
+    value1 = value1 or (value1 shr 8)
+    value2 = value2 or (value2 shr 8)
+    value1 = value1 and 0x0000ffffu
+    value2 = value2 and 0x0000ffffu
+    return uintArrayOf(value1, value2)
+}

orx-shapes/src/commonMain/kotlin/ordering/Morton3d.kt

@@ -0,0 +1,154 @@
+package org.openrndr.extra.shapes.ordering
+/*
+https://and-what-happened.blogspot.com/2011/08/fast-2d-and-3d-hilbert-curves-and.html
+ */
+
+
+fun morton3dEncode5Bit(
+    index1: UInt,
+    index2: UInt,
+    index3: UInt
+): UInt { // pack 3 5-bit indices into a 15-bit Morton code
+    var index1 = index1
+    var index2 = index2
+    var index3 = index3
+
+
+    index1 = index1 and 0x0000001fU
+    index2 = index2 and 0x0000001fU
+    index3 = index3 and 0x0000001fU
+    index1 *= 0x01041041U
+    index2 *= 0x01041041U
+    index3 *= 0x01041041U
+    index1 = index1 and 0x10204081U
+    index2 = index2 and 0x10204081U
+    index3 = index3 and 0x10204081U
+    index1 *= 0x00011111U
+    index2 *= 0x00011111U
+    index3 *= 0x00011111U
+    index1 = index1 and 0x12490000U
+    index2 = index2 and 0x12490000U
+    index3 = index3 and 0x12490000U
+    return ((index1 shr 16) or (index2 shr 15) or (index3 shr 14));
+}
+
+fun morton3dDecode5Bit(morton: UInt): UIntArray { // unpack 3 5-bit indices from a 15-bit Morton code
+    var value1 = morton;
+    var value2 = (value1 shr 1);
+    var value3 = (value1 shr 2);
+    value1 = value1 and 0x00001249U
+    value2 = value2 and 0x00001249U
+    value3 = value3 and 0x00001249U
+    value1 = value1 or (value1 shr 2);
+    value2 = value2 or (value2 shr 2);
+    value3 = value3 or (value3 shr 2);
+    value1 = value1 and 0x000010c3U
+    value2 = value2 and 0x000010c3U
+    value3 = value3 and 0x000010c3U
+    value1 = value1 or (value1 shr 4);
+    value2 = value2 or (value2 shr 4);
+    value3 = value3 or (value3 shr 4);
+    value1 = value1 and 0x0000100fU
+    value2 = value2 and 0x0000100fU
+    value3 = value3 and 0x0000100fU
+    value1 = value1 or (value1 shr 8);
+    value2 = value2 or (value2 shr 8);
+    value3 = value3 or (value3 shr 8);
+    value1 = value1 and 0x0000001fU
+    value2 = value2 and 0x0000001fU
+    value3 = value3 and 0x0000001fU
+    return uintArrayOf(value1, value2, value3)
+}
+
+
+fun morton3dEncode10Bit(
+    index1: UInt,
+    index2: UInt,
+    index3: UInt
+): UInt { // pack 3 10-bit indices into a 30-bit Morton code
+
+    var index1 = index1
+    var index2 = index2
+    var index3 = index3
+
+    index1 = index1 and 0x000003ffU
+    index2 = index2 and 0x000003ffU
+    index3 = index3 and 0x000003ffU
+    index1 = index1 or (index1 shl 16)
+    index2 = index2 or (index2 shl 16)
+    index3 = index3 or (index3 shl 16)
+    index1 = index1 and 0x030000ffU
+    index2 = index2 and 0x030000ffU
+    index3 = index3 and 0x030000ffU
+    index1 = index1 or (index1 shl 8)
+    index2 = index2 or (index2 shl 8)
+    index3 = index3 or (index3 shl 8)
+    index1 = index1 and 0x0300f00fU
+    index2 = index2 and 0x0300f00fU
+    index3 = index3 and 0x0300f00fU
+    index1 = index1 or (index1 shl 4)
+    index2 = index2 or (index2 shl 4)
+    index3 = index3 or (index3 shl 4)
+    index1 = index1 and 0x030c30c3U
+    index2 = index2 and 0x030c30c3U
+    index3 = index3 and 0x030c30c3U
+    index1 = index1 or (index1 shl 2)
+    index2 = index2 or (index2 shl 2)
+    index3 = index3 or (index3 shl 2)
+    index1 = index1 and 0x09249249U
+    index2 = index2 and 0x09249249U
+    index3 = index3 and 0x09249249U
+    return (index1 or (index2 shl 1) or (index3 shl 2))
+}
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun morton3dDecode10Bit(morton: UInt): UIntArray { // unpack 3 10-bit indices from a 30-bit Morton code
+    var value1 = morton
+    var value2 = (value1 shr 1)
+    var value3 = (value1 shr 2)
+    value1 = value1 and 0x09249249U
+    value2 = value2 and 0x09249249U
+    value3 = value3 and 0x09249249U
+    value1 = value1 or (value1 shr 2)
+    value2 = value2 or (value2 shr 2)
+    value3 = value3 or (value3 shr 2)
+    value1 = value1 and 0x030c30c3U
+    value2 = value2 and 0x030c30c3U
+    value3 = value3 and 0x030c30c3U
+    value1 = value1 or (value1 shr 4)
+    value2 = value2 or (value2 shr 4)
+    value3 = value3 or (value3 shr 4)
+    value1 = value1 and 0x0300f00fU
+    value2 = value2 and 0x0300f00fU
+    value3 = value3 and 0x0300f00fU
+    value1 = value1 or (value1 shr 8)
+    value2 = value2 or (value2 shr 8)
+    value3 = value3 or (value3 shr 8)
+    value1 = value1 and 0x030000ffU
+    value2 = value2 and 0x030000ffU
+    value3 = value3 and 0x030000ffU
+    value1 = value1 or (value1 shr 16)
+    value2 = value2 or (value2 shr 16)
+    value3 = value3 or (value3 shr 16)
+    value1 = value1 and 0x000003ffU;
+    value2 = value2 and 0x000003ffU;
+    value3 = value3 and 0x000003ffU;
+    return uintArrayOf(value1, value2, value3)
+}
+
+
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun main() {
+    val decoded = morton3dDecode10Bit(300U)
+    val encoded = morton3dEncode10Bit(decoded[0], decoded[1], decoded[2])
+    println(decoded)
+    println(encoded)
+
+    run {
+        val decoded = hilbert3dDecode10Bit(300U)
+        val encoded = hilbert3dEncode10Bit(decoded[0], decoded[1], decoded[2])
+        println(decoded)
+        println(encoded)
+    }
+}

orx-shapes/src/commonMain/kotlin/ordering/MortonToHilbert2d.kt

@@ -0,0 +1,29 @@
+package org.openrndr.extra.shapes.ordering
+
+fun mortonToHilbert2d(morton: UInt, bits: Int): UInt {
+    var hilbert = 0u
+    var remap = 0xb4u
+    var block = (bits shl 1)
+    while (block != 0) {
+        block -= 2
+        var mcode = ((morton shr block) and 3u)
+        var hcode = ((remap shr (mcode shl 1).toInt()) and 3u)
+        remap = remap xor (0x82000028u shr (hcode shl 3).toInt())
+        hilbert = ((hilbert shl 2) + hcode)
+    }
+    return (hilbert);
+}
+
+fun hilbertToMorton2d(hilbert: UInt, bits: Int): UInt {
+    var morton = 0u
+    var remap = 0xb4u
+    var block = (bits shl 1)
+    while (block != 0) {
+        block -= 2;
+        var hcode = ((hilbert shr block) and 3u)
+        var mcode = ((remap shr (hcode shl 1).toInt()) and 3u)
+        remap = remap xor (0x330000ccu shr (hcode shl 3).toInt())
+        morton = ((morton shl 2) + mcode)
+    }
+    return (morton)
+}

orx-shapes/src/commonMain/kotlin/ordering/MortonToHilbert3d.kt

@@ -0,0 +1,57 @@
+package org.openrndr.extra.shapes.ordering
+
+/* https://threadlocalmutex.com/?p=149 */
+@OptIn(ExperimentalUnsignedTypes::class)
+private val mortonToHilbertTable = uintArrayOf(
+    48u, 33u, 27u, 34u, 47u, 78u, 28u, 77u,
+    66u, 29u, 51u, 52u, 65u, 30u, 72u, 63u,
+    76u, 95u, 75u, 24u, 53u, 54u, 82u, 81u,
+    18u, 3u, 17u, 80u, 61u, 4u, 62u, 15u,
+    0u, 59u, 71u, 60u, 49u, 50u, 86u, 85u,
+    84u, 83u, 5u, 90u, 79u, 56u, 6u, 89u,
+    32u, 23u, 1u, 94u, 11u, 12u, 2u, 93u,
+    42u, 41u, 13u, 14u, 35u, 88u, 36u, 31u,
+    92u, 37u, 87u, 38u, 91u, 74u, 8u, 73u,
+    46u, 45u, 9u, 10u, 7u, 20u, 64u, 19u,
+    70u, 25u, 39u, 16u, 69u, 26u, 44u, 43u,
+    22u, 55u, 21u, 68u, 57u, 40u, 58u, 67u
+)
+
+@OptIn(ExperimentalUnsignedTypes::class)
+private val hilbertToMortonTable = uintArrayOf(
+    48u, 33u, 35u, 26u, 30u, 79u, 77u, 44u,
+    78u, 68u, 64u, 50u, 51u, 25u, 29u, 63u,
+    27u, 87u, 86u, 74u, 72u, 52u, 53u, 89u,
+    83u, 18u, 16u, 1u, 5u, 60u, 62u, 15u,
+    0u, 52u, 53u, 57u, 59u, 87u, 86u, 66u,
+    61u, 95u, 91u, 81u, 80u, 2u, 6u, 76u,
+    32u, 2u, 6u, 12u, 13u, 95u, 91u, 17u,
+    93u, 41u, 40u, 36u, 38u, 10u, 11u, 31u,
+    14u, 79u, 77u, 92u, 88u, 33u, 35u, 82u,
+    70u, 10u, 11u, 23u, 21u, 41u, 40u, 4u,
+    19u, 25u, 29u, 47u, 46u, 68u, 64u, 34u,
+    45u, 60u, 62u, 71u, 67u, 18u, 16u, 49u
+)
+
+@OptIn(ExperimentalUnsignedTypes::class)
+private fun transformCurve(input: UInt, bits: Int, lookupTable: UIntArray): UInt {
+    var transform = 0u
+    var out = 0u
+    for (i in 3 * (bits - 1) downTo 0 step 3) {
+        transform = lookupTable[(transform or ((input shr i) and 7U)).toInt()];
+        out = (out shl 3) or (transform and 7U)
+        transform = transform and (7U).inv()
+    }
+    return out;
+}
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun mortonToHilbert3d(mortonIndex: UInt, bits: Int = 10): UInt {
+    return transformCurve(mortonIndex, bits, mortonToHilbertTable)
+}
+
+
+@OptIn(ExperimentalUnsignedTypes::class)
+fun hilbertToMorton3d(hilbertIndex: UInt, bits: Int = 10): UInt {
+    return transformCurve(hilbertIndex, bits, hilbertToMortonTable)
+}