ImplicitFunction.h

Go to the documentation of this file.

//============================================================================
//  The contents of this file are covered by the Viskores license. See
//  LICENSE.txt for details.
//
//  By contributing to this file, all contributors agree to the Developer
//  Certificate of Origin Version 1.1 (DCO 1.1) as stated in DCO.txt.
//============================================================================
 
//============================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt for details.
//
//  This software is distributed WITHOUT ANY WARRANTY; without even
//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
//  PURPOSE.  See the above copyright notice for more information.
//============================================================================
#ifndef viskores_ImplicitFunction_h
#define viskores_ImplicitFunction_h
 
#include <viskores/Bounds.h>
#include <viskores/Math.h>
#include <viskores/VecVariable.h>
#include <viskores/VectorAnalysis.h>
#include <viskores/exec/Variant.h>
 
// For interface class only.
#include <viskores/cont/ExecutionAndControlObjectBase.h>
 
namespace viskores
{
 
//============================================================================
namespace internal
{
 
template <typename Derived>
class ImplicitFunctionBase : public viskores::cont::ExecutionAndControlObjectBase
{
public:
  using Scalar = viskores::FloatDefault;
  using Vector = viskores::Vec<Scalar, 3>;
 
  VISKORES_EXEC_CONT Scalar Value(Scalar x, Scalar y, Scalar z) const
  {
    return reinterpret_cast<const Derived*>(this)->Value(Vector(x, y, z));
  }
 
  VISKORES_EXEC_CONT Vector Gradient(Scalar x, Scalar y, Scalar z) const
  {
    return reinterpret_cast<const Derived*>(this)->Gradient(Vector(x, y, z));
  }
 
  VISKORES_CONT Derived PrepareForExecution(viskores::cont::DeviceAdapterId,
                                            viskores::cont::Token&) const
  {
    return *reinterpret_cast<const Derived*>(this);
  }
 
  VISKORES_CONT Derived PrepareForControl() const
  {
    return *reinterpret_cast<const Derived*>(this);
  }
};
 
} // namespace viskores::internal
 
//============================================================================
template <typename FunctionType>
class ImplicitFunctionValueFunctor
{
public:
  using Scalar = typename FunctionType::Scalar;
  using Vector = typename FunctionType::Vector;
 
  ImplicitFunctionValueFunctor() = default;
 
  VISKORES_EXEC_CONT ImplicitFunctionValueFunctor(
    const viskores::internal::ImplicitFunctionBase<FunctionType>& function)
    : Function(reinterpret_cast<const FunctionType&>(function))
  {
  }
 
  VISKORES_EXEC_CONT ImplicitFunctionValueFunctor(const FunctionType& function)
    : Function(function)
  {
  }
 
  VISKORES_EXEC_CONT Scalar operator()(const Vector& point) const
  {
    return this->Function.Value(point);
  }
 
private:
  FunctionType Function;
};
 
template <typename FunctionType>
class ImplicitFunctionGradientFunctor
{
public:
  using Scalar = typename FunctionType::Scalar;
  using Vector = typename FunctionType::Vector;
 
  ImplicitFunctionGradientFunctor() = default;
 
  VISKORES_EXEC_CONT ImplicitFunctionGradientFunctor(
    const viskores::internal::ImplicitFunctionBase<FunctionType>& function)
    : Function(reinterpret_cast<const FunctionType&>(function))
  {
  }
 
  VISKORES_EXEC_CONT ImplicitFunctionGradientFunctor(const FunctionType& function)
    : Function(function)
  {
  }
 
  VISKORES_EXEC_CONT Vector operator()(const Vector& point) const
  {
    return this->Function->Gradient(point);
  }
 
private:
  FunctionType Function;
};
 
//============================================================================
 
class VISKORES_ALWAYS_EXPORT Box : public internal::ImplicitFunctionBase<Box>
{
public:
  VISKORES_EXEC_CONT Box()
    : MinPoint(Vector(Scalar(-0.5)))
    , MaxPoint(Vector(Scalar(0.5)))
  {
  }
 
  VISKORES_EXEC_CONT Box(const Vector& minPoint, const Vector& maxPoint)
    : MinPoint(minPoint)
    , MaxPoint(maxPoint)
  {
  }
 
  VISKORES_EXEC_CONT Box(Scalar xmin,
                         Scalar xmax,
                         Scalar ymin,
                         Scalar ymax,
                         Scalar zmin,
                         Scalar zmax)
    : MinPoint(xmin, ymin, zmin)
    , MaxPoint(xmax, ymax, zmax)
  {
  }
 
  VISKORES_CONT Box(const viskores::Bounds& bounds) { this->SetBounds(bounds); }
 
  VISKORES_CONT void SetMinPoint(const Vector& point) { this->MinPoint = point; }
 
  VISKORES_CONT void SetMaxPoint(const Vector& point) { this->MaxPoint = point; }
 
  VISKORES_EXEC_CONT const Vector& GetMinPoint() const { return this->MinPoint; }
 
  VISKORES_EXEC_CONT const Vector& GetMaxPoint() const { return this->MaxPoint; }
 
  VISKORES_CONT void SetBounds(const viskores::Bounds& bounds)
  {
    this->SetMinPoint({ Scalar(bounds.X.Min), Scalar(bounds.Y.Min), Scalar(bounds.Z.Min) });
    this->SetMaxPoint({ Scalar(bounds.X.Max), Scalar(bounds.Y.Max), Scalar(bounds.Z.Max) });
  }
 
  VISKORES_EXEC_CONT viskores::Bounds GetBounds() const
  {
    return viskores::Bounds(viskores::Range(this->MinPoint[0], this->MaxPoint[0]),
                            viskores::Range(this->MinPoint[1], this->MaxPoint[1]),
                            viskores::Range(this->MinPoint[2], this->MaxPoint[2]));
  }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    Scalar minDistance = viskores::NegativeInfinity32();
    Scalar diff, t, dist;
    Scalar distance = Scalar(0.0);
    viskores::IdComponent inside = 1;
 
    for (viskores::IdComponent d = 0; d < 3; d++)
    {
      diff = this->MaxPoint[d] - this->MinPoint[d];
      if (diff != Scalar(0.0))
      {
        t = (point[d] - this->MinPoint[d]) / diff;
        // Outside before the box
        if (t < Scalar(0.0))
        {
          inside = 0;
          dist = this->MinPoint[d] - point[d];
        }
        // Outside after the box
        else if (t > Scalar(1.0))
        {
          inside = 0;
          dist = point[d] - this->MaxPoint[d];
        }
        else
        {
          // Inside the box in lower half
          if (t <= Scalar(0.5))
          {
            dist = MinPoint[d] - point[d];
          }
          // Inside the box in upper half
          else
          {
            dist = point[d] - MaxPoint[d];
          }
          if (dist > minDistance)
          {
            minDistance = dist;
          }
        }
      }
      else
      {
        dist = viskores::Abs(point[d] - MinPoint[d]);
        if (dist > Scalar(0.0))
        {
          inside = 0;
        }
      }
      if (dist > Scalar(0.0))
      {
        distance += dist * dist;
      }
    }
 
    distance = viskores::Sqrt(distance);
    if (inside)
    {
      return minDistance;
    }
    else
    {
      return distance;
    }
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector& point) const
  {
    viskores::IdComponent minAxis = 0;
    Scalar dist = 0.0;
    Scalar minDist = viskores::Infinity32();
    viskores::IdComponent3 location;
    Vector normal(Scalar(0));
    Vector inside(Scalar(0));
    Vector outside(Scalar(0));
    Vector center((this->MaxPoint + this->MinPoint) * Scalar(0.5));
 
    // Compute the location of the point with respect to the box
    // Point will lie in one of 27 separate regions around or within the box
    // Gradient vector is computed differently in each of the regions.
    for (viskores::IdComponent d = 0; d < 3; d++)
    {
      if (point[d] < this->MinPoint[d])
      {
        // Outside the box low end
        location[d] = 0;
        outside[d] = -1.0;
      }
      else if (point[d] > this->MaxPoint[d])
      {
        // Outside the box high end
        location[d] = 2;
        outside[d] = 1.0;
      }
      else
      {
        location[d] = 1;
        if (point[d] <= center[d])
        {
          // Inside the box low end
          dist = point[d] - this->MinPoint[d];
          inside[d] = -1.0;
        }
        else
        {
          // Inside the box high end
          dist = this->MaxPoint[d] - point[d];
          inside[d] = 1.0;
        }
        if (dist < minDist) // dist is negative
        {
          minDist = dist;
          minAxis = d;
        }
      }
    }
 
    viskores::Id indx = location[0] + 3 * location[1] + 9 * location[2];
    switch (indx)
    {
      // verts - gradient points away from center point
      case 0:
      case 2:
      case 6:
      case 8:
      case 18:
      case 20:
      case 24:
      case 26:
        for (viskores::IdComponent d = 0; d < 3; d++)
        {
          normal[d] = point[d] - center[d];
        }
        viskores::Normalize(normal);
        break;
 
      // edges - gradient points out from axis of cube
      case 1:
      case 3:
      case 5:
      case 7:
      case 9:
      case 11:
      case 15:
      case 17:
      case 19:
      case 21:
      case 23:
      case 25:
        for (viskores::IdComponent d = 0; d < 3; d++)
        {
          if (outside[d] != 0.0)
          {
            normal[d] = point[d] - center[d];
          }
          else
          {
            normal[d] = 0.0;
          }
        }
        viskores::Normalize(normal);
        break;
 
      // faces - gradient points perpendicular to face
      case 4:
      case 10:
      case 12:
      case 14:
      case 16:
      case 22:
        for (viskores::IdComponent d = 0; d < 3; d++)
        {
          normal[d] = outside[d];
        }
        break;
 
      // interior - gradient is perpendicular to closest face
      case 13:
        normal[0] = normal[1] = normal[2] = 0.0;
        normal[minAxis] = inside[minAxis];
        break;
      default:
        VISKORES_ASSERT(false);
        break;
    }
    return normal;
  }
 
private:
  Vector MinPoint;
  Vector MaxPoint;
};
 
//============================================================================
class VISKORES_ALWAYS_EXPORT Cylinder : public viskores::internal::ImplicitFunctionBase<Cylinder>
{
public:
  VISKORES_EXEC_CONT Cylinder()
    : Center(Scalar(0))
    , Axis(Scalar(0), Scalar(1), Scalar(0))
    , Radius(Scalar(0.5))
  {
  }
 
  VISKORES_EXEC_CONT Cylinder(const Vector& axis, Scalar radius)
    : Center(Scalar(0))
    , Axis(axis)
    , Radius(radius)
  {
  }
 
  VISKORES_EXEC_CONT Cylinder(const Vector& center, const Vector& axis, Scalar radius)
    : Center(center)
    , Axis(viskores::Normal(axis))
    , Radius(radius)
  {
  }
 
  VISKORES_CONT void SetCenter(const Vector& center) { this->Center = center; }
 
  VISKORES_CONT void SetAxis(const Vector& axis) { this->Axis = viskores::Normal(axis); }
 
  VISKORES_CONT void SetRadius(Scalar radius) { this->Radius = radius; }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    Vector x2c = point - this->Center;
    FloatDefault proj = viskores::Dot(this->Axis, x2c);
    return viskores::Dot(x2c, x2c) - (proj * proj) - (this->Radius * this->Radius);
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector& point) const
  {
    Vector x2c = point - this->Center;
    FloatDefault t = this->Axis[0] * x2c[0] + this->Axis[1] * x2c[1] + this->Axis[2] * x2c[2];
    viskores::Vec<FloatDefault, 3> closestPoint = this->Center + (this->Axis * t);
    return (point - closestPoint) * FloatDefault(2);
  }
 
private:
  Vector Center;
  Vector Axis;
  Scalar Radius;
};
 
//============================================================================
class VISKORES_ALWAYS_EXPORT Frustum : public viskores::internal::ImplicitFunctionBase<Frustum>
{
public:
  Frustum() = default;
 
  VISKORES_EXEC_CONT Frustum(const Vector points[6], const Vector normals[6])
  {
    this->SetPlanes(points, normals);
  }
 
  VISKORES_EXEC_CONT explicit Frustum(const Vector points[8]) { this->CreateFromPoints(points); }
 
  VISKORES_EXEC void SetPlanes(const Vector points[6], const Vector normals[6])
  {
    for (viskores::Id index : { 0, 1, 2, 3, 4, 5 })
    {
      this->Points[index] = points[index];
    }
    for (viskores::Id index : { 0, 1, 2, 3, 4, 5 })
    {
      this->Normals[index] = normals[index];
    }
  }
 
  VISKORES_EXEC void SetPlane(int idx, const Vector& point, const Vector& normal)
  {
    VISKORES_ASSERT((idx >= 0) && (idx < 6));
    this->Points[idx] = point;
    this->Normals[idx] = normal;
  }
 
  VISKORES_EXEC_CONT void GetPlanes(Vector points[6], Vector normals[6]) const
  {
    for (viskores::Id index : { 0, 1, 2, 3, 4, 5 })
    {
      points[index] = this->Points[index];
    }
    for (viskores::Id index : { 0, 1, 2, 3, 4, 5 })
    {
      normals[index] = this->Normals[index];
    }
  }
 
  VISKORES_EXEC_CONT const Vector* GetPoints() const { return this->Points; }
 
  VISKORES_EXEC_CONT const Vector* GetNormals() const { return this->Normals; }
 
  VISKORES_EXEC_CONT void CreateFromPoints(const Vector points[8])
  {
    // XXX(clang-format-3.9): 3.8 is silly. 3.9 makes it look like this.
    // clang-format off
    int planes[6][3] = {
      { 3, 2, 0 }, { 4, 5, 7 }, { 0, 1, 4 }, { 1, 2, 5 }, { 2, 3, 6 }, { 3, 0, 7 }
    };
    // clang-format on
 
    for (int i = 0; i < 6; ++i)
    {
      const Vector& v0 = points[planes[i][0]];
      const Vector& v1 = points[planes[i][1]];
      const Vector& v2 = points[planes[i][2]];
 
      this->Points[i] = v0;
      this->Normals[i] = viskores::Normal(viskores::TriangleNormal(v0, v1, v2));
    }
  }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    Scalar maxVal = viskores::NegativeInfinity<Scalar>();
    for (viskores::Id index : { 0, 1, 2, 3, 4, 5 })
    {
      const Vector& p = this->Points[index];
      const Vector& n = this->Normals[index];
      const Scalar val = viskores::Dot(point - p, n);
      maxVal = viskores::Max(maxVal, val);
    }
    return maxVal;
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector& point) const
  {
    Scalar maxVal = viskores::NegativeInfinity<Scalar>();
    viskores::Id maxValIdx = 0;
    for (viskores::Id index : { 0, 1, 2, 3, 4, 5 })
    {
      const Vector& p = this->Points[index];
      const Vector& n = this->Normals[index];
      Scalar val = viskores::Dot(point - p, n);
      if (val > maxVal)
      {
        maxVal = val;
        maxValIdx = index;
      }
    }
    return this->Normals[maxValIdx];
  }
 
private:
  Vector Points[6] = { { -0.5f, 0.0f, 0.0f }, { 0.5f, 0.0f, 0.0f },  { 0.0f, -0.5f, 0.0f },
                       { 0.0f, 0.5f, 0.0f },  { 0.0f, 0.0f, -0.5f }, { 0.0f, 0.0f, 0.5f } };
  Vector Normals[6] = { { -1.0f, 0.0f, 0.0f }, { 1.0f, 0.0f, 0.0f },  { 0.0f, -1.0f, 0.0f },
                        { 0.0f, 1.0f, 0.0f },  { 0.0f, 0.0f, -1.0f }, { 0.0f, 0.0f, 1.0f } };
};
 
//============================================================================
class VISKORES_ALWAYS_EXPORT Plane : public viskores::internal::ImplicitFunctionBase<Plane>
{
public:
  VISKORES_EXEC_CONT explicit Plane(const Vector& normal = { 0, 0, 1 })
    : Origin(Scalar(0))
    , Normal(normal)
  {
  }
 
  VISKORES_EXEC_CONT Plane(const Vector& origin, const Vector& normal)
    : Origin(origin)
    , Normal(normal)
  {
  }
 
  VISKORES_CONT void SetOrigin(const Vector& origin) { this->Origin = origin; }
 
  VISKORES_CONT void SetNormal(const Vector& normal) { this->Normal = normal; }
 
  VISKORES_EXEC_CONT const Vector& GetOrigin() const { return this->Origin; }
  VISKORES_EXEC_CONT const Vector& GetNormal() const { return this->Normal; }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    return viskores::Dot(point - this->Origin, this->Normal);
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector&) const { return this->Normal; }
 
private:
  Vector Origin;
  Vector Normal;
};
 
//============================================================================
class VISKORES_ALWAYS_EXPORT Sphere : public viskores::internal::ImplicitFunctionBase<Sphere>
{
public:
  VISKORES_EXEC_CONT explicit Sphere(Scalar radius = 0.5)
    : Radius(radius)
    , Center(Scalar(0))
  {
  }
 
  VISKORES_EXEC_CONT Sphere(Vector center, Scalar radius)
    : Radius(radius)
    , Center(center)
  {
  }
 
  VISKORES_CONT void SetRadius(Scalar radius) { this->Radius = radius; }
 
  VISKORES_CONT void SetCenter(const Vector& center) { this->Center = center; }
 
  VISKORES_EXEC_CONT Scalar GetRadius() const { return this->Radius; }
 
  VISKORES_EXEC_CONT const Vector& GetCenter() const { return this->Center; }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    return viskores::MagnitudeSquared(point - this->Center) - (this->Radius * this->Radius);
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector& point) const
  {
    return Scalar(2) * (point - this->Center);
  }
 
private:
  Scalar Radius;
  Vector Center;
};
 
//============================================================================
template <viskores::IdComponent MaxNumPlanes>
class VISKORES_ALWAYS_EXPORT MultiPlane
  : public viskores::internal::ImplicitFunctionBase<MultiPlane<MaxNumPlanes>>
{
public:
  using Scalar = viskores::FloatDefault;
  using Vector = viskores::Vec<Scalar, 3>;
  VISKORES_CONT MultiPlane() = default;
  template <viskores::IdComponent SrcMaxPlanes>
  VISKORES_CONT MultiPlane(const MultiPlane<SrcMaxPlanes>& src)
    : Planes(src.GetPlanes())
  {
  }
  template <viskores::IdComponent SrcMaxPlanes>
  VISKORES_CONT MultiPlane& operator=(const MultiPlane<SrcMaxPlanes>& src)
  {
    this->Planes = viskores::VecVariable<viskores::Plane, MaxNumPlanes>{ src.GetPlanes() };
  }
  VISKORES_CONT void AddPlane(const Vector& origin, const Vector& normal)
  {
    VISKORES_ASSERT(this->Planes.GetNumberOfComponents() < MaxNumPlanes);
    this->Planes.Append(Plane(origin, normal));
  }
  VISKORES_CONT viskores::Plane GetPlane(int idx)
  {
    VISKORES_ASSERT((idx >= 0) && (idx < MaxNumPlanes));
    return this->Planes[idx];
  }
  VISKORES_CONT viskores::VecVariable<viskores::Plane, MaxNumPlanes> GetPlanes() const
  {
    return this->Planes;
  }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    Scalar maxVal = viskores::NegativeInfinity<Scalar>();
    viskores::IdComponent NumPlanes = this->Planes.GetNumberOfComponents();
    for (viskores::IdComponent index = 0; index < NumPlanes; ++index)
    {
      const Vector& p = this->Planes[index].GetOrigin();
      const Vector& n = this->Planes[index].GetNormal();
      const Scalar val = viskores::Dot(point - p, n);
      maxVal = viskores::Max(maxVal, val);
    }
    return maxVal;
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector& point) const
  {
    Scalar maxVal = viskores::NegativeInfinity<Scalar>();
    viskores::IdComponent maxValIdx = 0;
    viskores::IdComponent NumPlanes = Planes.GetNumberOfComponents();
    for (viskores::IdComponent index = 0; index < NumPlanes; ++index)
    {
      const Vector& p = this->Planes[index].GetOrigin();
      const Vector& n = this->Planes[index].GetNormal();
      Scalar val = viskores::Dot(point - p, n);
      if (val > maxVal)
      {
        maxVal = val;
        maxValIdx = index;
      }
    }
    return this->Planes[maxValIdx].GetNormal();
  }
 
private:
  viskores::VecVariable<viskores::Plane, MaxNumPlanes> Planes;
};
 
namespace detail
{
 
struct ImplicitFunctionValueFunctor
{
  template <typename ImplicitFunctionType>
  VISKORES_EXEC_CONT typename ImplicitFunctionType::Scalar operator()(
    const ImplicitFunctionType& function,
    const typename ImplicitFunctionType::Vector& point) const
  {
    return function.Value(point);
  }
};
 
struct ImplicitFunctionGradientFunctor
{
  template <typename ImplicitFunctionType>
  VISKORES_EXEC_CONT typename ImplicitFunctionType::Vector operator()(
    const ImplicitFunctionType& function,
    const typename ImplicitFunctionType::Vector& point) const
  {
    return function.Gradient(point);
  }
};
 
} // namespace detail
 
//============================================================================
template <typename... ImplicitFunctionTypes>
class ImplicitFunctionMultiplexer
  : public viskores::internal::ImplicitFunctionBase<
      ImplicitFunctionMultiplexer<ImplicitFunctionTypes...>>
{
  viskores::exec::Variant<ImplicitFunctionTypes...> Variant;
 
  using Superclass =
    viskores::internal::ImplicitFunctionBase<ImplicitFunctionMultiplexer<ImplicitFunctionTypes...>>;
 
public:
  using Scalar = typename Superclass::Scalar;
  using Vector = typename Superclass::Vector;
 
  ImplicitFunctionMultiplexer() = default;
 
  template <typename FunctionType>
  VISKORES_EXEC_CONT ImplicitFunctionMultiplexer(
    const viskores::internal::ImplicitFunctionBase<FunctionType>& function)
    : Variant(reinterpret_cast<const FunctionType&>(function))
  {
  }
 
  VISKORES_EXEC_CONT Scalar Value(const Vector& point) const
  {
    return this->Variant.CastAndCall(detail::ImplicitFunctionValueFunctor{}, point);
  }
 
  VISKORES_EXEC_CONT Vector Gradient(const Vector& point) const
  {
    return this->Variant.CastAndCall(detail::ImplicitFunctionGradientFunctor{}, point);
  }
};
 
//============================================================================
class ImplicitFunctionGeneral
  : public viskores::ImplicitFunctionMultiplexer<viskores::Box,
                                                 viskores::Cylinder,
                                                 viskores::Frustum,
                                                 viskores::Plane,
                                                 viskores::Sphere,
                                                 viskores::MultiPlane<3>>
{
  using Superclass = viskores::ImplicitFunctionMultiplexer<viskores::Box,
                                                           viskores::Cylinder,
                                                           viskores::Frustum,
                                                           viskores::Plane,
                                                           viskores::Sphere,
                                                           viskores::MultiPlane<3>>;
 
public:
  using Superclass::Superclass;
};
 
} // namespace viskores
 
#endif //viskores_ImplicitFunction_h