EBGeometry_AnalyticDistanceFunctions.hpp

Go to the documentation of this file.

/* EBGeometry
 * Copyright © 2022 Robert Marskar
 * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory.
 */
 
#ifndef EBGeometry_AnalyticDistanceFunctions
#define EBGeometry_AnalyticDistanceFunctions
 
// Std includes
#include <algorithm>
#include <random>
 
// Our includes
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_SignedDistanceFunction.hpp"
#include "EBGeometry_NamespaceHeader.hpp"
 
template <class T>
constexpr const T
clamp(const T& v, const T& lo, const T& hi)
{
  return (v < lo) ? lo : (v > hi) ? hi : v;
}
 
template <class T>
constexpr const Vec3T<T>
clamp(const Vec3T<T>& v, const Vec3T<T>& lo, const Vec3T<T>& hi)
{
  return Vec3T<T>(clamp(v[0], lo[0], hi[0]), clamp(v[1], lo[1], hi[1]), clamp(v[2], lo[2], hi[2]));
}
 
constexpr static std::array<int, 256> s_perlinPermutationTable = {
  151, 160, 137, 91,  90,  15,  131, 13,  201, 95,  96,  53,  194, 233, 7,   225, 140, 36,  103, 30,  69,  142,
  8,   99,  37,  240, 21,  10,  23,  190, 6,   148, 247, 120, 234, 75,  0,   26,  197, 62,  94,  252, 219, 203,
  117, 35,  11,  32,  57,  177, 33,  88,  237, 149, 56,  87,  174, 20,  125, 136, 171, 168, 68,  175, 74,  165,
  71,  134, 139, 48,  27,  166, 77,  146, 158, 231, 83,  111, 229, 122, 60,  211, 133, 230, 220, 105, 92,  41,
  55,  46,  245, 40,  244, 102, 143, 54,  65,  25,  63,  161, 1,   216, 80,  73,  209, 76,  132, 187, 208, 89,
  18,  169, 200, 196, 135, 130, 116, 188, 159, 86,  164, 100, 109, 198, 173, 186, 3,   64,  52,  217, 226, 250,
  124, 123, 5,   202, 38,  147, 118, 126, 255, 82,  85,  212, 207, 206, 59,  227, 47,  16,  58,  17,  182, 189,
  28,  42,  223, 183, 170, 213, 119, 248, 152, 2,   44,  154, 163, 70,  221, 153, 101, 155, 167, 43,  172, 9,
  129, 22,  39,  253, 19,  98,  108, 110, 79,  113, 224, 232, 178, 185, 112, 104, 218, 246, 97,  228, 251, 34,
  242, 193, 238, 210, 144, 12,  191, 179, 162, 241, 81,  51,  145, 235, 249, 14,  239, 107, 49,  192, 214, 31,
  181, 199, 106, 157, 184, 84,  204, 176, 115, 121, 50,  45,  127, 4,   150, 254, 138, 236, 205, 93,  222, 114,
  67,  29,  24,  72,  243, 141, 128, 195, 78,  66,  215, 61,  156, 180};
 
template <class T>
class PlaneSDF : public SignedDistanceFunction<T>
{
public:
  PlaneSDF() = delete;
 
  PlaneSDF(const Vec3T<T>& a_point, const Vec3T<T>& a_normal) noexcept
  {
    m_point  = a_point;
    m_normal = a_normal / a_normal.length();
  }
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    return dot((a_point - m_point), m_normal);
  }
 
protected:
  Vec3T<T> m_point;
 
  Vec3T<T> m_normal;
};
 
template <class T>
class SphereSDF : public SignedDistanceFunction<T>
{
public:
  SphereSDF() = delete;
 
  SphereSDF(const Vec3T<T>& a_center, const T& a_radius) noexcept
  {
    this->m_center = a_center;
    this->m_radius = a_radius;
  }
 
  SphereSDF(const SphereSDF& a_other) noexcept
  {
    this->m_center = a_other.m_center;
    this->m_radius = a_other.m_radius;
  }
 
  virtual ~SphereSDF() noexcept = default;
 
  const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  Vec3T<T>&
  getCenter() noexcept
  {
    return m_center;
  }
 
  const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  T&
  getRadius() noexcept
  {
    return m_radius;
  }
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    return (a_point - m_center).length() - m_radius;
  }
 
protected:
  Vec3T<T> m_center;
 
  T m_radius;
};
 
template <class T>
class BoxSDF : public SignedDistanceFunction<T>
{
public:
  BoxSDF() = delete;
 
  BoxSDF(const Vec3T<T>& a_loCorner, const Vec3T<T>& a_hiCorner) noexcept
  {
    this->m_loCorner = a_loCorner;
    this->m_hiCorner = a_hiCorner;
  }
 
  virtual ~BoxSDF() noexcept
  {}
 
  const Vec3T<T>&
  getLowCorner() const noexcept
  {
    return m_loCorner;
  }
 
  Vec3T<T>&
  getLowCorner() noexcept
  {
    return m_loCorner;
  }
 
  const Vec3T<T>&
  getHighCorner() const noexcept
  {
    return m_hiCorner;
  }
 
  Vec3T<T>&
  getHighCorner() noexcept
  {
    return m_hiCorner;
  }
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    // For each coordinate direction, we have delta[dir] if a_point[dir] falls
    // between xLo and xHi. In this case delta[dir] will be the signed distance
    // to the closest box face in the dir-direction. Otherwise, if a_point[dir]
    // is outside the corner we have delta[dir] > 0.
    const Vec3T<T> delta(std::max(m_loCorner[0] - a_point[0], a_point[0] - m_hiCorner[0]),
                         std::max(m_loCorner[1] - a_point[1], a_point[1] - m_hiCorner[1]),
                         std::max(m_loCorner[2] - a_point[2], a_point[2] - m_hiCorner[2]));
 
    // Note: max is max(Vec3T<T>, Vec3T<T>) and not std::max. It returns a
    // vector with coordinate-wise largest components. Note that the first part
    // std::min(...) is the signed distance on the inside of the box (delta will
    // have negative components). The other part max(Vec3T<T>::zero(), ...) is
    // for outside the box.
    const T d = std::min(T(0.0), delta[delta.maxDir(false)]) + max(Vec3T<T>::zero(), delta).length();
 
    return d;
  }
 
protected:
  Vec3T<T> m_loCorner;
 
  Vec3T<T> m_hiCorner;
};
 
template <class T>
class TorusSDF : public SignedDistanceFunction<T>
{
public:
  TorusSDF() = delete;
 
  TorusSDF(const Vec3T<T>& a_center, const T& a_majorRadius, const T& a_minorRadius) noexcept
  {
    this->m_center      = a_center;
    this->m_majorRadius = a_majorRadius;
    this->m_minorRadius = a_minorRadius;
  }
 
  virtual ~TorusSDF() noexcept
  {}
 
  const Vec3T<T>&
  getCenter() const noexcept
  {
    return m_center;
  }
 
  Vec3T<T>&
  getCenter() noexcept
  {
    return m_center;
  }
 
  const T&
  getMajorRadius() const noexcept
  {
    return m_majorRadius;
  }
 
  T&
  getMajorRadius() noexcept
  {
    return m_majorRadius;
  }
 
  const T&
  getMinorRadius() const noexcept
  {
    return m_minorRadius;
  }
 
  T&
  getMinorRadius() noexcept
  {
    return m_minorRadius;
  }
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    const Vec3T<T> p   = a_point - m_center;
    const T        rho = sqrt(p[0] * p[0] + p[1] * p[1]) - m_majorRadius;
    const T        d   = sqrt(rho * rho + p[2] * p[2]) - m_minorRadius;
 
    return d;
  }
 
protected:
  Vec3T<T> m_center;
 
  T m_majorRadius;
 
  T m_minorRadius;
};
 
template <class T>
class CylinderSDF : public SignedDistanceFunction<T>
{
public:
  CylinderSDF() = delete;
 
  CylinderSDF(const Vec3T<T>& a_center1, const Vec3T<T>& a_center2, const T& a_radius) noexcept
  {
    this->m_center1 = a_center1;
    this->m_center2 = a_center2;
    this->m_radius  = a_radius;
 
    // Some derived quantities that are needed for SDF computations.
    m_center = (m_center2 + m_center1) * 0.5;
    m_length = (m_center2 - m_center1).length();
    m_axis   = (m_center2 - m_center1) / m_length;
  }
 
  virtual ~CylinderSDF() noexcept
  {}
 
  const Vec3T<T>&
  getCenter1() const noexcept
  {
    return m_center1;
  }
 
  const Vec3T<T>&
  getCenter2() const noexcept
  {
    return m_center2;
  }
 
  const T&
  getRadius() const noexcept
  {
    return m_radius;
  }
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    T d = std::numeric_limits<T>::infinity();
 
    if (m_length > 0.0 && m_radius > 0.0) {
      const Vec3T<T> point = a_point - m_center;
      const T        para  = dot(point, m_axis);
      const Vec3T<T> ortho = point - para * m_axis; // Distance from cylinder axis.
 
      const T w = ortho.length() - m_radius;       // Distance from cylinder wall. < 0
                                                   // on inside and > 0 on outside.
      const T h = std::abs(para) - 0.5 * m_length; // Distance from cylinder top.  < 0 on
                                                   // inside and > 0 on outside.
 
      constexpr T zero = T(0.0);
 
      if (w <= zero && h <= zero) { // Inside cylinder
        d = (std::abs(w) < std::abs(h)) ? w : h;
      }
      else if (w <= zero && h > zero) { // Above one of the endcaps.
        d = h;
      }
      else if (w > zero && h < zero) { // Outside radius but between the
                                       // endcaps.
        d = w;
      }
      else {
        d = sqrt(w * w + h * h);
      }
    }
 
    return d;
  }
 
protected:
  Vec3T<T> m_center1;
 
  Vec3T<T> m_center2;
 
  Vec3T<T> m_center;
 
  Vec3T<T> m_axis;
 
  T m_length;
 
  T m_radius;
};
 
template <class T>
class InfiniteCylinderSDF : public SignedDistanceFunction<T>
{
public:
  InfiniteCylinderSDF() = delete;
 
  InfiniteCylinderSDF(const Vec3T<T>& a_center, const T& a_radius, const size_t a_axis) noexcept
  {
    m_center = a_center;
    m_radius = a_radius;
    m_axis   = a_axis;
  }
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    Vec3T<T> delta = a_point - m_center;
    delta[m_axis]  = 0.0;
 
    const T d = delta.length() - m_radius;
 
    return d;
  }
 
protected:
  Vec3T<T> m_center;
 
  T m_radius;
 
  size_t m_axis;
};
 
template <class T>
class CapsuleSDF : public SignedDistanceFunction<T>
{
public:
  CapsuleSDF() = delete;
 
  CapsuleSDF(const Vec3T<T>& a_tip1, const Vec3T<T> a_tip2, const T& a_radius) noexcept
  {
    const Vec3T<T> axis = (a_tip2 - a_tip1) / length(a_tip2 - a_tip1);
    m_center1           = a_tip1 + a_radius * axis;
    m_center2           = a_tip2 - a_radius * axis;
    m_radius            = a_radius;
  }
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    const Vec3T<T> v1 = a_point - m_center1;
    const Vec3T<T> v2 = m_center2 - m_center1;
 
    const T h = clamp(dot(v1, v2) / dot(v2, v2), T(0.0), T(1.0));
    const T d = length(v1 - h * v2) - m_radius;
 
    return d;
  }
 
protected:
  Vec3T<T> m_center1;
 
  Vec3T<T> m_center2;
 
  T m_radius;
};
 
template <class T>
class InfiniteConeSDF : public SignedDistanceFunction<T>
{
public:
  InfiniteConeSDF() = delete;
 
  InfiniteConeSDF(const Vec3T<T>& a_tip, const T& a_angle) noexcept
  {
    constexpr T pi = 3.14159265358979323846;
 
    m_tip = a_tip;
    m_c.x = std::sin(0.5 * a_angle * pi / 180.0);
    m_c.y = std::cos(0.5 * a_angle * pi / 180.0);
  }
 
  virtual ~InfiniteConeSDF()
  {}
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    const Vec3T<T> delta = a_point - m_tip;
    const Vec2T<T> q(sqrt(delta[0] * delta[0] + delta[1] * delta[1]), -delta[2]);
 
    const T d1 = length(q - m_c * std::max(dot(q, m_c), T(0.0)));
    const T d2 = d1 * ((q.x * m_c.y - q.y * m_c.x < 0.0) ? -1.0 : 1.0);
 
    return d2;
  }
 
protected:
  Vec3T<T> m_tip;
 
  Vec2T<T> m_c;
};
 
template <class T>
class ConeSDF : public SignedDistanceFunction<T>
{
public:
  ConeSDF() = delete;
 
  ConeSDF(const Vec3T<T>& a_tip, const T& a_height, const T& a_angle) noexcept
  {
    constexpr T pi = 3.14159265358979323846;
 
    m_tip    = a_tip;
    m_height = a_height;
    m_c.x    = std::sin(0.5 * a_angle * pi / 180.0);
    m_c.y    = std::cos(0.5 * a_angle * pi / 180.0);
  }
 
  virtual ~ConeSDF()
  {}
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    const Vec3T<T> delta = a_point - m_tip;
    const T        dr    = sqrt(delta[0] * delta[0] + delta[1] * delta[1]);
    const T        dz    = delta[2];
 
    constexpr T zero = T(0.0);
    constexpr T one  = T(1.0);
 
    const Vec2T<T> q = m_height * Vec2T<T>(m_c.x / m_c.y, -1.0);
    const Vec2T<T> w = Vec2T<T>(dr, dz);
    const Vec2T<T> a = w - clamp(dot(w, q) / dot(q, q), zero, one) * q;
    const Vec2T<T> b = w - Vec2T<T>(q.x * clamp(w.x / q.x, zero, one), q.y);
 
    auto sign = [](const T& x) { return (x > zero) - (x < zero); };
 
    const T k = sign(q.y);
    const T d = std::min(dot(a, a), dot(b, b));
    const T s = std::max(k * (w.x * q.y - w.y * q.x), k * (w.y - q.y));
 
    return sqrt(d) * sign(s);
  }
 
protected:
  Vec3T<T> m_tip;
 
  Vec2T<T> m_c;
 
  T m_height;
};
 
template <class T>
class RoundedBoxSDF : public SignedDistanceFunction<T>
{
public:
  RoundedBoxSDF() = delete;
 
  RoundedBoxSDF(const Vec3T<T>& a_dimensions, const T a_curvature) noexcept
  {
    this->m_dimensions = 0.5 * a_dimensions;
 
    m_sphere = std::make_shared<SphereSDF<T>>(Vec3T<T>::zero(), a_curvature);
  }
 
  virtual ~RoundedBoxSDF() noexcept
  {}
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    return m_sphere->value(a_point - clamp(a_point, -m_dimensions, m_dimensions));
  }
 
protected:
  std::shared_ptr<SphereSDF<T>> m_sphere;
 
  Vec3T<T> m_dimensions;
};
 
template <class T>
class PerlinSDF : public SignedDistanceFunction<T>
{
public:
  PerlinSDF(const T            a_noiseAmplitude,
            const Vec3T<T>     a_noiseFrequency,
            const T            a_noisePersistence,
            const unsigned int a_noiseOctaves) noexcept
  {
    m_noiseAmplitude   = a_noiseAmplitude;
    m_noiseFrequency   = a_noiseFrequency;
    m_noisePersistence = std::min(1.0, a_noisePersistence);
    m_noiseOctaves     = std::max(static_cast<unsigned int>(1), a_noiseOctaves);
 
    // By default, use Ken Perlin's original permutation table
    for (int i = 0; i < 256; i++) {
      m_permutationTable[i]       = s_perlinPermutationTable[i];
      m_permutationTable[i + 256] = s_perlinPermutationTable[i];
    }
 
#if 1 // Development code
    std::random_device rd;
    std::mt19937       g(0);
 
    this->shuffle(g);
#endif
  }
 
  virtual ~PerlinSDF() noexcept
  {}
 
  virtual T
  signedDistance(const Vec3T<T>& a_point) const noexcept override
  {
    T ret = 0.0;
 
    Vec3T<T> curFreq = m_noiseFrequency;
    T        curAmp  = 1.0;
 
    for (unsigned int curOctave = 0; curOctave < m_noiseOctaves; curOctave++) {
      ret += 0.5 * curAmp * (1 + this->noise(a_point * curFreq));
 
      curFreq = curFreq / m_noisePersistence;
      curAmp  = curAmp * m_noisePersistence;
    }
 
    return ret * m_noiseAmplitude;
  };
 
  template <class URNG>
  void
  shuffle(URNG& g) noexcept
  {
 
    for (unsigned int i = 0; i < 256; i++) {
      m_permutationTable[i] = i;
    }
 
    std::shuffle(m_permutationTable.begin(), m_permutationTable.begin() + 256, g);
 
    for (int i = 0; i < 256; i++) {
      m_permutationTable[i + 256] = m_permutationTable[i];
    }
  }
 
  std::array<int, 512>&
  getPermutationTable() noexcept
  {
    return m_permutationTable;
  }
 
protected:
  Vec3T<T> m_noiseFrequency;
 
  T m_noiseAmplitude;
 
  T m_noisePersistence;
 
  std::array<T, 512> m_permutationTable;
 
  unsigned int m_noiseOctaves;
 
  virtual T
  lerp(const T t, const T a, const T b) const noexcept
  {
    return a + t * (b - a);
  };
 
  virtual T
  fade(const T t) const noexcept
  {
    return t * t * t * (t * (t * 6 - 15) + 10);
  };
 
  T
  grad(const int hash, const T x, const T y, const T z) const noexcept
  {
    const int h = hash & 15;
    const T   u = h < 8 ? x : y;
    const T   v = h < 4 ? y : h == 12 || h == 14 ? x : y;
    return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v);
  }
 
  T
  noise(const Vec3T<T>& a_point) const noexcept
  {
    // Lower cube corner
    const int X = static_cast<int>(std::floor(a_point[0])) & 255;
    const int Y = static_cast<int>(std::floor(a_point[1])) & 255;
    const int Z = static_cast<int>(std::floor(a_point[2])) & 255;
 
    // Relative distance wrt lower cube corner
    const double x = a_point[0] - std::floor(a_point[0]);
    const double y = a_point[1] - std::floor(a_point[1]);
    const double z = a_point[2] - std::floor(a_point[2]);
 
    // Fade curves
    const double u = fade(x);
    const double v = fade(y);
    const double w = fade(z);
 
    // Hash coordinates of 8 cube corners
    const int A  = m_permutationTable[X] + Y;
    const int AA = m_permutationTable[A] + Z;
    const int AB = m_permutationTable[A + 1] + Z;
    const int B  = m_permutationTable[X + 1] + Y;
    const int BA = m_permutationTable[B] + Z;
    const int BB = m_permutationTable[B + 1] + Z;
 
    // Add blended results from 8 corners of cube
    return lerp(
      w,
      lerp(v,
           lerp(u, grad(m_permutationTable[AA], x, y, z), grad(m_permutationTable[BA], x - 1, y, z)),
           lerp(u, grad(m_permutationTable[AB], x, y - 1, z), grad(m_permutationTable[BB], x - 1, y - 1, z))),
      lerp(v,
           lerp(u, grad(m_permutationTable[AA + 1], x, y, z - 1), grad(m_permutationTable[BA + 1], x - 1, y, z - 1)),
           lerp(u,
                grad(m_permutationTable[AB + 1], x, y - 1, z - 1),
                grad(m_permutationTable[BB + 1], x - 1, y - 1, z - 1))));
  };
};
 
#include "EBGeometry_NamespaceFooter.hpp"
 
#endif