{"id":157,"date":"2018-04-27T20:16:49","date_gmt":"2018-04-27T20:16:49","guid":{"rendered":"http:\/\/threadlocalmutex.com\/?p=157"},"modified":"2018-04-27T20:16:49","modified_gmt":"2018-04-27T20:16:49","slug":"proof-of-concept-3d-hilbert-curves-in-ologn","status":"publish","type":"post","link":"http:\/\/threadlocalmutex.com\/?p=157","title":{"rendered":"Proof of concept: 3D Hilbert curves in O(log(n))"},"content":{"rendered":"

Back in my earlier post<\/a> about 2D Hilbert curves in logarithmic time, I hinted that it would be possible but likely not efficient to extend the same approach to higher dimensions. Just as a proof of concept, I’ve done just that for the Hilbert Index -> XYZ mapping in 3D. Most of the actual code is auto-generated, and I haven’t spent much effort optimizing it in the same way as the 2D version, as it is rather obviously useless in practice. I’m only presenting the 32bit version here, but the code easily extends to arbitrary bit widths with an appropriate large integer library.<\/p>\n

Applying the same would also be possible for the XYZ -> Hilbert Index mapping, however the code would be at least 10x-20x slower due to the exponential blowup of the technique, and probably rather tedious to generate.<\/p>\n

\r\n\r\n\/\/ Shamelessly borrowed from http:\/\/and-what-happened.blogspot.de\/2011\/08\/fast-2d-and-3d-hilbert-curves-and.html\r\nvoid Morton_3D_Decode_10bit(const uint32_t morton,\r\n\tuint32_t& index1, uint32_t& index2, uint32_t& index3)\r\n{ \/\/ unpack 3 10-bit indices from a 30-bit Morton code\r\n\tuint32_t value1 = morton;\r\n\tuint32_t value2 = (value1 >> 1);\r\n\tuint32_t value3 = (value1 >> 2);\r\n\tvalue1 &= 0x09249249;\r\n\tvalue2 &= 0x09249249;\r\n\tvalue3 &= 0x09249249;\r\n\tvalue1 |= (value1 >> 2);\r\n\tvalue2 |= (value2 >> 2);\r\n\tvalue3 |= (value3 >> 2);\r\n\tvalue1 &= 0x030c30c3;\r\n\tvalue2 &= 0x030c30c3;\r\n\tvalue3 &= 0x030c30c3;\r\n\tvalue1 |= (value1 >> 4);\r\n\tvalue2 |= (value2 >> 4);\r\n\tvalue3 |= (value3 >> 4);\r\n\tvalue1 &= 0x0300f00f;\r\n\tvalue2 &= 0x0300f00f;\r\n\tvalue3 &= 0x0300f00f;\r\n\tvalue1 |= (value1 >> 8);\r\n\tvalue2 |= (value2 >> 8);\r\n\tvalue3 |= (value3 >> 8);\r\n\tvalue1 &= 0x030000ff;\r\n\tvalue2 &= 0x030000ff;\r\n\tvalue3 &= 0x030000ff;\r\n\tvalue1 |= (value1 >> 16);\r\n\tvalue2 |= (value2 >> 16);\r\n\tvalue3 |= (value3 >> 16);\r\n\tvalue1 &= 0x000003ff;\r\n\tvalue2 &= 0x000003ff;\r\n\tvalue3 &= 0x000003ff;\r\n\tindex1 = value1;\r\n\tindex2 = value2;\r\n\tindex3 = value3;\r\n}\r\n\r\nvoid transformFromIndex_parallel(uint32_t i0, uint32_t i1, uint32_t i2, uint32_t& t0, uint32_t& t1, uint32_t& t2, uint32_t& t3)\r\n{\r\n\tt0 = i2 | (i0 & i1);\r\n\tt1 = (~i0 & ~i1 & ~i2) | (i0 & i1 & ~i2) | (~i0 & ~i1 & i2);\r\n\tt2 = ~(i0 & i1) ^ i2;\r\n\tt3 = (i0 & ~i1 & i2) | (~i0 & i1 & i2);\r\n}\r\n\r\nvoid multiplyTransform_parallel(\r\n\tuint32_t A0, uint32_t A1, uint32_t A2, uint32_t A3,\r\n\tuint32_t B0, uint32_t B1, uint32_t B2, uint32_t B3,\r\n\tuint32_t& C0, uint32_t& C1, uint32_t& C2, uint32_t& C3)\r\n{\r\n\tC0 = (A0 & ~A1 & ~A2 & B0 & B1 & ~B2 & ~B3) | (A0 & ~A2 & ~A3 & ~B0 & B2) | (A0 & ~A1 & ~A2 & ~A3 & B0 & ~B1 & B3) | (A0 & ~A2 & ~A3 & ~B0 & B1 & B3) | (A0 & A1 & A3 & ~B0 & B1 & ~B2 & ~B3) | (A1 & ~A3 & ~B0 & B1 & ~B2 & ~B3) | (~A0 & A1 & ~A2 & ~A3 & B1 & B3) | (A1 & A2 & B0 & ~B1 & ~B2 & ~B3) | (A0 & ~A1 & A3 & B0 & ~B2 & ~B3) | (~A0 & ~A1 & A3 & B0 & B1 & ~B3) | (A1 & ~A2 & ~A3 & ~B0 & ~B1 & B2) | (~A0 & A1 & A2 & B0 & B2) | (~A0 & A1 & ~A2 & ~B0 & ~B1 & B3) | (A0 & A1 & ~A2 & ~B0 & ~B1 & ~B2) | (~A0 & A1 & ~A3 & B0 & ~B1 & B2) | (~A0 & A2 & ~B0 & ~B1 & B3) | (~A0 & A1 & A3 & ~B0 & B2) | (~A0 & ~A1 & ~A2 & ~A3 & B0) | (A0 & A1 & A3 & ~B1 & ~B3) | (A0 & ~B0 & ~B1 & ~B2 & ~B3) | (A0 & A1 & ~A3 & ~B0 & ~B3) | (~A0 & A2 & B0 & B1 & B3) | (~A1 & A3 & B0 & ~B1 & B3) | (~A0 & A3 & B0 & ~B1 & B3) | (~A0 & ~A1 & A2 & B3) | (A0 & A2 & ~B2 & ~B3) | (~A1 & A2 & B0 & B2) | (~A0 & A3 & B1 & B3) | (~A1 & A3 & B1 & B2);\r\n\tC1 = (~A0 & A3 & B0 & B1 & ~B2 & ~B3) | (A0 & ~A1 & ~A2 & ~A3 & B0 & ~B1 & B3) | (~A0 & ~A2 & ~A3 & B0 & B1 & B2) | (A0 & A1 & A3 & ~B0 & B1 & ~B2 & ~B3) | (A0 & ~A3 & B0 & ~B1 & ~B2 & ~B3) | (~A0 & ~A1 & A3 & B0 & ~B2 & ~B3) | (~A1 & A2 & ~B0 & B1 & ~B2 & ~B3) | (A0 & ~A1 & A3 & B0 & ~B2 & ~B3) | (A0 & ~A1 & A3 & B1 & ~B2 & ~B3) | (~A1 & ~A2 & ~A3 & ~B0 & B1 & B2) | (~A0 & A1 & A2 & B0 & B2) | (A0 & A1 & ~A2 & ~B0 & ~B1 & ~B2) | (~A0 & A1 & ~A3 & B0 & ~B1 & B2) | (~A0 & A2 & B1 & ~B2 & ~B3) | (~A0 & ~A1 & A2 & ~B1 & B2) | (A1 & ~A2 & ~A3 & ~B1 & B3) | (A0 & ~A1 & ~A2 & ~B0 & B2) | (~A0 & ~A3 & ~B0 & B1 & B3) | (~A0 & ~A1 & ~A3 & B1 & ~B2) | (A0 & A1 & A3 & B1 & B2) | (A0 & ~A3 & ~B0 & ~B1 & B3) | (A1 & ~B0 & ~B1 & ~B2 & ~B3) | (A0 & A1 & ~A3 & ~B1 & ~B2) | (A0 & A1 & ~A2 & ~A3 & ~B1) | (~A0 & A3 & B0 & ~B1 & B3) | (A1 & A3 & ~B0 & B1 & B3) | (~A0 & A3 & ~B1 & B2) | (A0 & A2 & B1 & B2) | (~A1 & B0 & B1 & B3);\r\n\tC2 = (A0 & ~A1 & A2 & ~B1 & ~B2 & ~B3) | (~A0 & A1 & ~A2 & B0 & B1 & ~B2 & ~B3) | (A0 & A1 & ~A2 & B0 & ~B1 & ~B2 & ~B3) | (A0 & ~A1 & ~A2 & B0 & B1 & ~B2 & ~B3) | (A0 & A1 & A2 & ~B0 & ~B2 & ~B3) | (~A0 & ~A1 & A3 & B0 & ~B1 & B2) | (~A1 & ~A2 & ~A3 & B0 & ~B1 & B2) | (~A0 & A1 & A2 & ~B0 & B1 & B2) | (~A0 & ~A1 & A3 & ~B0 & B1 & ~B3) | (A1 & A3 & B0 & ~B1 & ~B2 & ~B3) | (~A0 & A1 & A2 & B0 & ~B1 & B2) | (A0 & ~A1 & ~A2 & ~B0 & ~B1 & B3) | (~A0 & ~A1 & ~A3 & ~B0 & ~B1 & B2) | (A1 & ~A2 & ~A3 & B0 & B1 & B2) | (~A0 & ~A1 & ~A3 & B0 & B1 & B2) | (A0 & A1 & A3 & ~B0 & B1 & ~B2 & ~B3) | (A0 & A1 & ~A2 & ~B0 & B1 & B3) | (A0 & A1 & A3 & B0 & ~B1 & ~B2) | (A0 & ~A1 & A3 & B1 & ~B2 & ~B3) | (~A1 & ~A2 & ~A3 & ~B0 & B1 & B2) | (A1 & ~A2 & ~A3 & ~B0 & ~B1 & B2) | (~A0 & A1 & ~A2 & ~B0 & ~B1 & B3) | (A2 & ~B0 & ~B1 & ~B2 & ~B3) | (A0 & ~A1 & A2 & B1 & B3) | (A0 & A1 & A2 & B0 & B3) | (~A0 & A2 & B0 & B1 & B3) | (~A1 & A3 & B0 & ~B1 & B3) | (A1 & A3 & ~B0 & B1 & B3);\r\n\tC3 = (A0 & A1 & ~A2 & ~A3 & ~B0 & B1 & ~B2 & ~B3) | (A0 & ~A1 & ~A2 & ~A3 & B0 & B1 & B2) | (A0 & A1 & ~A2 & ~A3 & ~B0 & ~B1 & B3) | (A0 & A1 & A3 & B0 & B1 & ~B2 & ~B3) | (~A0 & ~A1 & A2 & B0 & ~B1 & ~B3) | (A0 & ~A1 & ~A2 & ~B0 & B1 & B3) | (~A0 & A1 & ~A2 & ~B0 & B1 & B2) | (A0 & ~A1 & ~A2 & ~A3 & B0 & ~B1 & B3) | (~A0 & A1 & ~A3 & B0 & ~B1 & B3) | (A0 & A1 & ~A3 & B0 & ~B1 & B2) | (~A0 & ~A1 & A3 & B0 & ~B2 & ~B3) | (~A1 & A2 & ~B0 & B1 & ~B2 & ~B3) | (~A0 & A1 & ~A2 & ~A3 & B1 & B3) | (A1 & A2 & B0 & ~B1 & ~B2 & ~B3) | (~A0 & ~A1 & A3 & B0 & B1 & ~B3) | (~A0 & ~A1 & ~A2 & ~A3 & B3) | (~A1 & A2 & ~B0 & ~B1 & B3) | (A3 & ~B0 & ~B1 & ~B2 & ~B3) | (A0 & ~A1 & A2 & ~B1 & B2) | (~A2 & ~A3 & B0 & B1 & B3) | (~A0 & A1 & A3 & ~B0 & ~B3) | (A1 & A2 & B0 & B1 & ~B3) | (A0 & A3 & ~B0 & B2);\r\n}\r\n\r\nvoid octantFromIndex_parallel(uint32_t i0, uint32_t i1, uint32_t i2, uint32_t& x, uint32_t& y, uint32_t& z)\r\n{\r\n\tx = i0 ^ i1;\r\n\ty = i1 ^ i2;\r\n\tz = i2;\r\n}\r\n\r\nvoid applyTransform_parallel(uint32_t t0, uint32_t t1, uint32_t t2, uint32_t t3, uint32_t x, uint32_t y, uint32_t z, uint32_t& o0, uint32_t& o1, uint32_t& o2)\r\n{\r\n\to2 = (~t3 | t1 | t0 | ~z)&(t3 | t1 | ~t0 | ~y)&(~t3 | ~t2 | ~z | ~y | x)&(~t3 | t2 | t1 | ~t0 | ~x)&(t3 | ~t2 | t1 | t0 | x)&(t3 | ~t2 | ~t1 | t0 | y)&(~t3 | t2 | ~t1 | t0 | y)&(t3 | ~t2 | ~t1 | ~t0 | ~x)&(~t3 | t2 | ~t1 | ~t0 | ~z)&(t3 | t2 | t1 | t0 | z)&(t3 | t2 | ~t1 | t0 | x)&(t3 | t2 | ~t1 | ~t0 | z)&(~t3 | ~t2 | t1 | ~z | y | x)&(~t3 | ~t2 | t1 | z | y | ~x)&(~t3 | ~t2 | t0 | z | ~y | ~x)&(~t3 | ~t2 | ~t1 | ~t0 | y | x)&(~t3 | t0 | ~z | y | ~x)&(~t2 | t1 | t0 | ~y | x);\r\n\to1 = (~t3 | ~t2 | ~y | ~x)&(~t3 | ~t1 | t0 | ~x)&(~t3 | t1 | t0 | ~y)&(~t3 | ~t2 | t1 | ~z | ~x)&(t3 | ~t2 | z | x)&(~t3 | ~t2 | t1 | z | ~y)&(~t2 | t1 | t0 | z | ~x)&(~t3 | ~t2 | t0 | z | ~y)&(~t2 | ~t1 | ~t0 | z | ~x)&(~t3 | ~t2 | ~t1 | y | x)&(~t3 | t2 | t1 | ~t0 | ~z)&(~t3 | t2 | ~t1 | ~t0 | y)&(t3 | t2 | t1 | t0 | y)&(t3 | t2 | t1 | ~t0 | ~x)&(t2 | ~t1 | t0 | ~z | ~x)&(t3 | ~t1 | t0 | ~z | x)&(t3 | t2 | ~t1 | ~t0 | ~y)&(~t3 | ~t2 | t1 | ~z | y | x)&(t3 | ~t2 | t1 | ~t0 | ~z | x)&(~t3 | ~t2 | ~t0 | z | y | x);\r\n\to0 = (~t3 | ~t2 | ~y | ~x)&(~t2 | ~t1 | t0 | z)&(t2 | t1 | t0 | x)&(t2 | ~t1 | ~t0 | ~x)&(~t3 | t0 | ~z | x)&(~t3 | ~t2 | t1 | ~z | ~x)&(~t3 | ~t2 | t1 | z | ~y)&(~t3 | ~t2 | ~z | ~y | x)&(~t3 | ~t2 | ~t1 | z | ~x)&(t3 | ~t2 | t1 | t0 | y)&(~t3 | t2 | t1 | ~t0 | y)&(~t2 | t1 | ~t0 | ~z | ~x)&(~t3 | t2 | ~t1 | ~z | ~x)&(t3 | ~t2 | ~t1 | ~t0 | ~y)&(t3 | t2 | t1 | ~t0 | z)&(t3 | t2 | ~t1 | t0 | ~y)&(t3 | ~t2 | t1 | ~t0 | ~z | x)&(~t3 | ~t2 | ~t0 | z | y | x)&(~t2 | t0 | z | y | ~x);\r\n}\r\n\r\nvoid hilbertToMorton3D_logarithmic(uint32_t hilbertIndex, uint32_t bits, uint32_t& x, uint32_t& y, uint32_t& z)\r\n{\r\n\tuint32_t i0, i1, i2;\r\n\tMorton_3D_Decode_10bit(hilbertIndex, i0, i1, i2);\r\n\ti0 <<= 32 - bits;\r\n\ti1 <<= 32 - bits;\r\n\ti2 <<= 32 - bits;\r\n\r\n\tuint32_t t0, t1, t2, t3;\r\n\ttransformFromIndex_parallel(i0, i1, i2, t0, t1, t2, t3);\r\n\r\n\t{\r\n\t\tt0 >>= 1;\r\n\t\tt1 >>= 1;\r\n\t\tt2 >>= 1;\r\n\t\tt3 >>= 1;\r\n\t}\r\n\r\n\tfor (uint32_t shift = 1; shift <= 16; shift *= 2)\r\n\t{\r\n\t\tuint32_t s0, s1, s2, s3;\r\n\r\n\t\ts0 = t0 >> shift;\r\n\t\ts1 = t1 >> shift;\r\n\t\ts2 = t2 >> shift;\r\n\t\ts3 = t3 >> shift;\r\n\r\n\t\tuint32_t T0, T1, T2, T3;\r\n\t\tmultiplyTransform_parallel(t0, t1, t2, t3, s0, s1, s2, s3, T0, T1, T2, T3);\r\n\r\n\t\tt0 = T0;\r\n\t\tt1 = T1;\r\n\t\tt2 = T2;\r\n\t\tt3 = T3;\r\n\t}\r\n\r\n\tuint32_t o0, o1, o2;\r\n\toctantFromIndex_parallel(i0, i1, i2, o0, o1, o2);\r\n\tapplyTransform_parallel(t0, t1, t2, t3, o0, o1, o2, x, y, z);\r\n\r\n\tx >>= 32 - bits;\r\n\ty >>= 32 - bits;\r\n\tz >>= 32 - bits;\r\n}\r\n\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"
Back in my earlier post about 2D Hilbert curves in logarithmic time, I hinted that it would be possible but likely not efficient to extend the same approach to higher dimensions. Just as a proof of concept, I’ve done just that for the Hilbert Index -> XYZ mapping in 3D. Most of the actual code […]<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-157","post","type-post","status-publish","format-standard","hentry","category-general"],"_links":{"self":[{"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=\/wp\/v2\/posts\/157","targetHints":{"allow":["GET"]}}],"collection":[{"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=157"}],"version-history":[{"count":3,"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=\/wp\/v2\/posts\/157\/revisions"}],"predecessor-version":[{"id":160,"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=\/wp\/v2\/posts\/157\/revisions\/160"}],"wp:attachment":[{"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=157"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=157"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/threadlocalmutex.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=157"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}