]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add inverse_Hilbert_space_filling_curve 7280/head
authorDenis Davydov <davydden@gmail.com>
Fri, 5 Oct 2018 09:53:11 +0000 (11:53 +0200)
committerDenis Davydov <davydden@gmail.com>
Sat, 6 Oct 2018 11:02:57 +0000 (13:02 +0200)
doc/news/changes/minor/20181005DenisDavydov [new file with mode: 0644]
include/deal.II/base/utilities.h
source/base/utilities.cc
tests/base/utilities_13.cc [new file with mode: 0644]
tests/base/utilities_13.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20181005DenisDavydov b/doc/news/changes/minor/20181005DenisDavydov
new file mode 100644 (file)
index 0000000..d8e1aa3
--- /dev/null
@@ -0,0 +1,4 @@
+New: Add inverse_Hilbert_space_filling_curve() to map points in dim-dimensional
+space to line index of the Hilbert space filling curve.
+<br>
+(Denis Davydov 2018/10/05)
index 44d83cd612d58383d9774fca88735d1a76611f33..1174563fa8bee0b388eda91efd041d32d4563932 100644 (file)
@@ -52,6 +52,9 @@
 
 DEAL_II_NAMESPACE_OPEN
 
+// forward declare Point
+template <int dim, typename Number>
+class Point;
 
 /**
  * A namespace for utility functions that are not particularly specific to
@@ -72,6 +75,34 @@ namespace Utilities
   std::string
   dealii_version_string();
 
+  /**
+   * Assign to each point in @p points an index using the Hilbert space filling curve.
+   * To that end, a bounding box for @p points will be determined, based on which their
+   * integer coordinates are calculated.
+   * The linear index is given as a dim-collection of bits, from high to low.
+   * This is done in order to keep the maximum resolution in terms of bit depth
+   * along each axis. Note that this dim-integer index can still be easily used
+   * for sorting and ordering, for example using the lexicographic ordering of
+   * tuples of integers.
+   *
+   * The depth of the Hilbert curve (i.e. the number of bits per dimension) by
+   * default is equal to <code>64</code>.
+   */
+  template <int dim, typename Number>
+  std::vector<std::array<std::uint64_t, dim>>
+  inverse_Hilbert_space_filling_curve(
+    const std::vector<Point<dim, Number>> &points,
+    const int                              bits_per_dim = 64);
+
+  /**
+   * Same as above, but for points in integer coordinates.
+   */
+  template <int dim>
+  std::vector<std::array<std::uint64_t, dim>>
+  inverse_Hilbert_space_filling_curve(
+    const std::vector<std::array<std::uint64_t, dim>> &points,
+    const int                                          bits_per_dim = 64);
+
   /**
    * Convert a number @p value to a string, with as many digits as given to
    * fill with leading zeros.
index d3bae9270b333eca6db39c68f0e9b6014e9ee509..baee0da0027a494e51cef9f02f49b7c71fc568be 100644 (file)
@@ -24,6 +24,7 @@
 
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/mpi.h>
+#include <deal.II/base/point.h>
 #include <deal.II/base/thread_local_storage.h>
 #include <deal.II/base/utilities.h>
 
@@ -64,8 +65,6 @@
 #  include <Teuchos_RCP.hpp>
 #endif
 
-
-
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -91,6 +90,184 @@ namespace Utilities
 
 
 
+  template <int dim, typename Number>
+  std::vector<std::array<std::uint64_t, dim>>
+  inverse_Hilbert_space_filling_curve(
+    const std::vector<Point<dim, Number>> &points,
+    const int                              bits_per_dim)
+  {
+    using Integer = std::uint64_t;
+    // take floating point number hopefully with mantissa >= 64bit
+    using LongDouble = long double;
+
+    // return if there is nothing to do
+    if (points.size() == 0)
+      return std::vector<std::array<std::uint64_t, dim>>();
+
+    // get bounding box:
+    Point<dim, Number> bl = points[0], tr = points[0];
+    for (const auto &p : points)
+      for (unsigned int d = 0; d < dim; ++d)
+        {
+          const double cid = p[d];
+          bl[d]            = std::min(cid, bl[d]);
+          tr[d]            = std::max(cid, tr[d]);
+        }
+
+    std::array<LongDouble, dim> extents;
+    for (unsigned int i = 0; i < dim; ++i)
+      {
+        extents[i] =
+          static_cast<LongDouble>(tr[i]) - static_cast<LongDouble>(bl[i]);
+        Assert(extents[i] > 0., ExcMessage("Bounding box is degenerate."));
+      }
+
+    // make sure our conversion from fractional coordinates to
+    // Integers work as expected, namely our cast (LongDouble)max_int
+    const int min_bits =
+      std::min(bits_per_dim,
+               std::min(std::numeric_limits<Integer>::digits,
+                        std::numeric_limits<LongDouble>::digits));
+
+    // based on that get the maximum integer:
+    const Integer max_int = (min_bits == std::numeric_limits<Integer>::digits ?
+                               std::numeric_limits<Integer>::max() :
+                               ((Integer)1 << min_bits) - 1);
+
+    std::vector<std::array<Integer, dim>> int_points(points.size());
+
+    for (unsigned int i = 0; i < points.size(); ++i)
+      {
+        // convert into integers:
+        for (unsigned int d = 0; d < dim; ++d)
+          {
+            const LongDouble v = (static_cast<LongDouble>(points[i][d]) -
+                                  static_cast<LongDouble>(bl[d])) /
+                                 extents[d];
+            Assert(v >= 0. && v <= 1., ExcInternalError());
+            int_points[i][d] =
+              static_cast<Integer>(v * static_cast<LongDouble>(max_int));
+          }
+      }
+
+    // note that we call this with "min_bits"
+    return inverse_Hilbert_space_filling_curve<dim>(int_points, min_bits);
+  }
+
+
+
+  template <int dim>
+  std::vector<std::array<std::uint64_t, dim>>
+  inverse_Hilbert_space_filling_curve(
+    const std::vector<std::array<std::uint64_t, dim>> &points,
+    const int                                          bits_per_dim)
+  {
+    using Integer = std::uint64_t;
+
+    std::vector<std::array<Integer, dim>> int_points(points);
+
+    std::vector<std::array<Integer, dim>> res(int_points.size());
+
+    // follow
+    // J. Skilling, Programming the Hilbert curve, AIP Conf. Proc. 707, 381
+    // (2004); http://dx.doi.org/10.1063/1.1751381 also see
+    // https://stackoverflow.com/questions/499166/mapping-n-dimensional-value-to-a-point-on-hilbert-curve
+    // https://gitlab.com/octopus-code/octopus/blob/develop/src/grid/hilbert.c
+    // https://github.com/trilinos/Trilinos/blob/master/packages/zoltan/src/hsfc/hsfc_hilbert.c
+    // (Zoltan_HSFC_InvHilbertXd)
+    // https://github.com/aditi137/Hilbert/blob/master/Hilbert/hilbert.cpp
+
+    // now we can map to 1D coordinate stored in Transpose format
+    // adopt AxestoTranspose function from the paper, that
+    // transforms in-place between geometrical axes and Hilbert transpose.
+    // Example:   b=5 bits for each of n=3 coordinates.
+    //            15-bit Hilbert integer = A B C D E F G H I J K L M N O is
+    //            stored as its Transpose
+    //                   X[0] = A D G J M                X[2]|
+    //                   X[1] = B E H K N    <------->       | /X[1]
+    //                   X[2] = C F I L O               axes |/
+    //                          high  low                    0------ X[0]
+
+    // Depth of the Hilbert curve
+    Assert(bits_per_dim <= std::numeric_limits<Integer>::digits,
+           ExcMessage("This integer type can not hold " +
+                      std::to_string(bits_per_dim) + " bits."));
+
+    const Integer M = 1 << (bits_per_dim - 1); // largest bit
+
+    for (unsigned int index = 0; index < int_points.size(); ++index)
+      {
+        auto &X = int_points[index];
+        auto &L = res[index];
+
+        // Inverse undo
+        for (Integer q = M; q > 1; q >>= 1)
+          {
+            const Integer p = q - 1;
+            for (unsigned int i = 0; i < dim; i++)
+              {
+                // invert
+                if (X[i] & q)
+                  {
+                    X[0] ^= p;
+                  }
+                // exchange
+                else
+                  {
+                    const Integer t = (X[0] ^ X[i]) & p;
+                    X[0] ^= t;
+                    X[i] ^= t;
+                  }
+              }
+          }
+
+        // Gray encode (inverse of decode)
+        for (unsigned int i = 1; i < dim; i++)
+          X[i] ^= X[i - 1];
+
+        Integer t = 0;
+        for (Integer q = M; q > 1; q >>= 1)
+          if (X[dim - 1] & q)
+            t ^= q - 1;
+        for (unsigned int i = 0; i < dim; i++)
+          X[i] ^= t;
+
+        // now we need to go from index stored in transpose format to
+        // consecutive format, which is better suited for comparators.
+        // we could interleave into some big unsigned int...
+        // https://www.forceflow.be/2013/10/07/morton-encodingdecoding-through-bit-interleaving-implementations/
+        // https://stackoverflow.com/questions/4431522/given-2-16-bit-ints-can-i-interleave-those-bits-to-form-a-single-32-bit-int
+        // ...but we would loose spatial resolution!
+
+        // interleave using brute force, follow TransposetoLine from
+        // https://github.com/aditi137/Hilbert/blob/master/Hilbert/hilbert.cpp
+        {
+          Integer      p = M;
+          unsigned int j = 0;
+          for (unsigned int i = 0; i < dim; ++i)
+            {
+              L[i] = 0;
+              // go through bits using a mask q
+              for (Integer q = M; q > 0; q >>= 1)
+                {
+                  if (X[j] & p)
+                    L[i] |= q;
+                  if (++j == dim)
+                    {
+                      j = 0;
+                      p >>= 1;
+                    }
+                }
+            }
+        }
+
+      } // end of the loop over points
+
+    return res;
+  }
+
+
+
   std::string
   int_to_string(const unsigned int value, const unsigned int digits)
   {
@@ -900,6 +1077,31 @@ namespace Utilities
   template std::string
   to_string<long double>(long double, unsigned int);
 
+  template std::vector<std::array<std::uint64_t, 1>>
+  inverse_Hilbert_space_filling_curve<1, double>(
+    const std::vector<Point<1, double>> &,
+    const int);
+  template std::vector<std::array<std::uint64_t, 1>>
+  inverse_Hilbert_space_filling_curve<1>(
+    const std::vector<std::array<std::uint64_t, 1>> &,
+    const int);
+  template std::vector<std::array<std::uint64_t, 2>>
+  inverse_Hilbert_space_filling_curve<2, double>(
+    const std::vector<Point<2, double>> &,
+    const int);
+  template std::vector<std::array<std::uint64_t, 2>>
+  inverse_Hilbert_space_filling_curve<2>(
+    const std::vector<std::array<std::uint64_t, 2>> &,
+    const int);
+  template std::vector<std::array<std::uint64_t, 3>>
+  inverse_Hilbert_space_filling_curve<3, double>(
+    const std::vector<Point<3, double>> &,
+    const int);
+  template std::vector<std::array<std::uint64_t, 3>>
+  inverse_Hilbert_space_filling_curve<3>(
+    const std::vector<std::array<std::uint64_t, 3>> &,
+    const int);
+
 } // namespace Utilities
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/utilities_13.cc b/tests/base/utilities_13.cc
new file mode 100644 (file)
index 0000000..10fa8b6
--- /dev/null
@@ -0,0 +1,145 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// test Utilities::inverse_Hilbert_space_filling_curve in 2D
+// on the following nine points with bits=4:
+
+//        |
+//     15 |    @---@   @---@   @---@   @---@   @---@   @---@   @---@   @---@
+//        |    |   |   |   |   |   |   |   |   |   |   |   |   |   |   |   |
+//        |    @   @---@   @   @   @---@   @   @   @---@   @   @   @---@   @
+//        |    |           |   |           |   |           |   |           |
+//        |    @---@   @---@   @---@   @---@   @---@   @---@   @---@   @---@
+//        |        |   |           |   |           |   |           |   |
+//        |    @---@   @---@---@---@   @---@   @---@   @---@---@---@   @---@
+//        |    |                           |   |                           |
+//        |    @   @---@---@   @---@---@   @   @   @---@---@   @---@---@   @
+//        |    |   |       |   |       |   |   |   |       |   |       |   |
+// Axes[1]|    @---@   @---@   @---@   @---@   @---@   @---@   @---@   @---@
+//        |            |           |                   |           |
+//        |    @---@   @---@   @---@   @---@   @---@   @---@   @---@   @---@
+//        |    |   |       |   |       |   |   |   |       |   |       |   |
+//        |    @   @---@---@   @---@---@   @---@   @---@---@   @---@---@   @
+//        |    |                                                           |
+//        |    @---@   @---@---@   @---@---@   @---@---@   @---@---@   @---@
+//        |        |   |       |   |       |   |       |   |       |   |
+//        |    @---@   @---@   @---@   @---@   @---@   @---@   @---@   @---@
+//        |    |           |           |           |           |           |
+//        |    @   @---@   @   @---@   @---@   @---@   @---@   @   @---@   @
+//        |    |   |   |   |   |   |       |   |       |   |   |   |   |   |
+//        |    @---@   @---@   @   @---@---@   @---@---@   @   @---@   @---@
+//        |                    |                           |
+//      3 |    5---6   9---@   @   @---@---@   @---@---@   @   @---@   @---@
+//        |    |   |   |   |   |   |       |   |       |   |   |   |   |   |
+//      2 |    4   7---8   @   @---@   @---@   @---@   @---@   @   @---@   @
+//        |    |           |           |           |           |           |
+//      1 |    3---2   @---@   @---@   @---@   @---@   @---@   @---@   @---@
+//        |        |   |       |   |       |   |       |   |       |   |
+//      0 |    0---1   @---@---@   @---@---@   @---@---@   @---@---@   @--255
+//        |
+//         -------------------------------------------------------------------
+//             0   1   2   3          ---> Axes[0]                         15
+
+
+#include <deal.II/base/utilities.h>
+
+#include "../tests.h"
+
+
+void
+test1()
+{
+  deallog << "test1:" << std::endl;
+  const std::vector<Point<2>> points = {Point<2>(0, 0),
+                                        Point<2>(1, 0),
+                                        Point<2>(1, 1),
+                                        Point<2>(0, 1),
+                                        Point<2>(0, 2),
+                                        Point<2>(0, 3),
+                                        Point<2>(1, 3),
+                                        Point<2>(1, 2),
+                                        Point<2>(2, 2),
+                                        Point<2>(2, 3),
+                                        Point<2>(15, 15),
+                                        Point<2>(15, 0)};
+
+  const auto res = Utilities::inverse_Hilbert_space_filling_curve(points, 4);
+
+  for (const auto &p : res)
+    deallog << p[0] << " " << p[1] << std::endl;
+}
+
+// https://github.com/aditi137/Hilbert/blob/master/Hilbert/test.py
+// test int version in 3D
+void
+test2()
+{
+  deallog << "test2:" << std::endl;
+  std::vector<std::array<std::uint64_t, 3>> points = {{5, 10, 20}};
+
+  const auto res = Utilities::inverse_Hilbert_space_filling_curve<3>(points, 5);
+
+  for (const auto &p : res)
+    deallog << p[0] << " " << p[1] << " " << p[2] << std::endl;
+}
+
+
+// finally, a quick test of fractional coordinates -> int conversion inside
+// inverse_Hilbert_space_filling_curve():
+template <typename Integer, typename Number>
+void
+test3()
+{
+  const int bits_int = std::numeric_limits<Integer>::digits;
+  const int bits_num = std::numeric_limits<Number>::digits;
+
+  const int min_bits = std::min(bits_int, bits_num);
+
+  const Integer max =
+    (min_bits == bits_int ? std::numeric_limits<Integer>::max() :
+                            ((Integer)1 << min_bits) - 1);
+
+  // make sure conversion is within bounds:
+  const std::vector<Number> points = {0.2, 0.7};
+  for (const auto &v : points)
+    {
+      const Integer i = (Integer)(v * (Number)max);
+      Assert(i > 0 && i < max, ExcInternalError());
+    }
+
+  {
+    const Number  v = 0.;
+    const Integer i = (Integer)(v * (Number)max);
+    Assert(i == 0, ExcInternalError());
+  }
+
+  {
+    const Number  v = 1.;
+    const Integer i = (Integer)(v * (Number)max);
+    Assert(i == max, ExcInternalError());
+  }
+}
+
+int
+main()
+{
+  initlog();
+
+  test1();
+  test2();
+
+  test3<std::uint64_t, long double>();
+}
diff --git a/tests/base/utilities_13.output b/tests/base/utilities_13.output
new file mode 100644 (file)
index 0000000..a5c97ad
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL::test1:
+DEAL::0 0
+DEAL::0 1
+DEAL::0 2
+DEAL::0 3
+DEAL::0 4
+DEAL::0 5
+DEAL::0 6
+DEAL::0 7
+DEAL::0 8
+DEAL::0 9
+DEAL::10 10
+DEAL::15 15
+DEAL::test2:
+DEAL::7 21 25

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.