Improve spatial_hash function to avoid spatial periodicity in the hash function

beluga/include/beluga/algorithm/spatial_hash.hpp

@@ -34,22 +34,31 @@ namespace beluga {
 
 namespace detail {
 
-/// Returns the floor of a value shifted.
+/// Returns the hashed and rotated floor of a value.
 /**
  * \tparam N Number of bits to be used from the integer result, the least significant will be used.
  * \tparam I Result will be shifted by I*N.
- * \param value Input value to be floored and shifted.
+ * \param value Input value to be hashed.
  * \return The calculated result.
  */
 template <std::size_t N, std::size_t I>
-constexpr std::size_t floor_and_shift(double value) {
-  // Compute the largest integer value not greater than value.
-  auto signed_value = static_cast<std::intmax_t>(std::floor(value));
-  // Compute the smallest unsigned integer equal to the source value modulo 2^n
-  // (where n is the number of bits used to represent the destination type).
-  auto unsigned_value = static_cast<std::uintmax_t>(signed_value);
-  // Create a fixed-size sequence of N bits and perform a binary shift left I * N positions.
-  return std::bitset<N>{unsigned_value}.to_ullong() << (I * N);
+constexpr std::size_t floor_and_fibo_hash(double value) {
+  static_assert(std::is_same_v<std::size_t, std::uint64_t>);
+  constexpr auto kFib64 = 11400714819323198485LLU;  // golden ratio for 64 bits
+  // floor the value and convert to integer
+  const auto signed_value = static_cast<std::int64_t>(std::floor(value));
+  // work with unsigned from now on
+  const auto unsigned_value = static_cast<std::uint64_t>(signed_value);
+  // spread number information all through the 64 bits using the fibonacci hash
+  const auto div_hashed_value = kFib64 * unsigned_value;
+  // rotate bits to avoid aliasing between different values of I
+  if constexpr (N * I != 0) {
+    const auto left_hash = (div_hashed_value << N * I);
+    const auto right_hash = (div_hashed_value >> (64 - N * I));
+    return left_hash | right_hash;
+  } else {
+    return div_hashed_value;
+  }
 }
 
 /// Hashes a tuple or array of scalar types, using a resolution for each element and using the same
@@ -68,7 +77,7 @@ constexpr std::size_t hash_impl(
     const std::array<double, sizeof...(Ids)>& resolution,
     [[maybe_unused]] std::index_sequence<Ids...> index_sequence) {
   constexpr auto kBits = std::numeric_limits<std::size_t>::digits / sizeof...(Ids);
-  return (detail::floor_and_shift<kBits, Ids>(std::get<Ids>(value) / resolution[Ids]) | ...);
+  return (detail::floor_and_fibo_hash<kBits, Ids>(std::get<Ids>(value) / resolution[Ids]) ^ ...);
 }
 
 }  // namespace detail

beluga/test/beluga/test_spatial_hash.cpp

@@ -98,4 +98,51 @@ TEST(SpatialHash, Array) {
   EXPECT_EQ(hash1, hash2);
 }
 
+TEST(SpatialHash, NonPeriodicityCheck1) {
+  using Array = std::array<double, 8>;
+  constexpr std::array kClusteringResolution{1., 1., 1., 1., 1., 1., 1., 1.};
+
+  auto uut = beluga::spatial_hash<Array>{kClusteringResolution};
+
+  auto hash0 = uut({0., 0., 0., 0., 0., 0., 0., 0.});
+  auto hash1 = uut({64., 0., 0., 0., 0., 0., 0., 0.});
+  auto hash2 = uut({128., 0., 0., 0., 0., 0., 0., 0.});
+  auto hash3 = uut({192., 0., 0., 0., 0., 0., 0., 0.});
+  auto hash4 = uut({256., 0., 0., 0., 0., 0., 0., 0.});
+
+  EXPECT_NE(hash0, hash1);
+  EXPECT_NE(hash0, hash2);
+  EXPECT_NE(hash0, hash3);
+  EXPECT_NE(hash0, hash4);
+}
+
+TEST(SpatialHash, NonPeriodicityCheck2) {
+  using Array = std::array<double, 8>;
+  constexpr std::array kClusteringResolution{1., 1., 1., 1., 1., 1., 1., 1.};
+  constexpr double kStep = 2.0;
+
+  auto uut = beluga::spatial_hash<Array>{kClusteringResolution};
+
+  auto ref_hash = uut({0., 0., 0., 0., 0., 0., 0., 0.});
+
+  for (double n = kStep; n < 1024.; n += kStep) {
+    auto hash0 = uut({n, 0., 0., 0., 0., 0., 0., 0.});
+    auto hash1 = uut({0., n, 0., 0., 0., 0., 0., 0.});
+    auto hash2 = uut({0., 0., n, 0., 0., 0., 0., 0.});
+    auto hash3 = uut({0., 0., 0., n, 0., 0., 0., 0.});
+    auto hash4 = uut({0., 0., 0., 0., n, 0., 0., 0.});
+    auto hash5 = uut({0., 0., 0., 0., 0., n, 0., 0.});
+    auto hash6 = uut({0., 0., 0., 0., 0., 0., n, 0.});
+    auto hash7 = uut({0., 0., 0., 0., 0., 0., 0., n});
+    EXPECT_NE(ref_hash, hash0);
+    EXPECT_NE(ref_hash, hash1);
+    EXPECT_NE(ref_hash, hash2);
+    EXPECT_NE(ref_hash, hash3);
+    EXPECT_NE(ref_hash, hash4);
+    EXPECT_NE(ref_hash, hash5);
+    EXPECT_NE(ref_hash, hash6);
+    EXPECT_NE(ref_hash, hash7);
+  }
+}
+
 }  // namespace