Commit d8cf1f7b authored by David Thompson's avatar David Thompson Committed by Kitware Robot
Browse files

Merge topic 'geometry-squashed'

880d8a98

 Add `vtkm/Geometry.h` and test it.
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Acked-by: default avatarRobert Maynard <robert.maynard@kitware.com>
Merge-request: !1262
parents 7de4bba7 880d8a98
# New geometry classes and header.
There are now some additional structures available both
the control and execution environments for representing
geometric entities (mostly of dimensions 2 and 3).
These new structures are now in `vtkm/Geometry.h` and
demonstrated/tested in `vtkm/testing/TestingGeometry.h`:
+ `Ray<CoordType, Dimension, IsTwoSided>`.
Instances of this struct represent a semi-infinite line
segment in a 2-D plane or in a 3-D space, depending on the
integer dimension specified as a template parameter.
Its state is the point at the start of the ray (`Origin`)
plus the ray's `Direction`, a unit-length vector.
If the third template parameter (IsTwoSided) is true, then
the ray serves as an infinite line. Otherwise, the ray will
only report intersections in its positive halfspace.
+ `LineSegment<CoordType, Dimension>`.
Instances of this struct represent a finite line segment
in a 2-D plane or in a 3-D space, depending on the integer
dimension specified as a template parameter.
Its state is the coordinates of its `Endpoints`.
+ `Plane<CoordType>`.
Instances of this struct represent a plane in 3-D.
Its state is the coordinates of a base point (`Origin`) and
a unit-length normal vector (`Normal`).
+ `Sphere<CoordType, Dimension>`.
Instances of this struct represent a *d*-dimensional sphere.
Its state is the coordinates of its center plus a radius.
It is also aliased with a `using` statment to `Circle<CoordType>`
for the specific case of 2-D.
These structures provide useful queries and generally
interact with one another.
For instance, it is possible to intersect lines and planes
and compute distances.
For ease of use, there are also several `using` statements
that alias these geometric structures to names that specialize
them for a particular dimension or other template parameter.
As an example, `Ray<CoordType, Dimension, true>` is aliased
to `Line<CoordType, Dimension>` and `Ray<CoordType, 3, true>`
is aliased to `Line3<CoordType>` and `Ray<FloatDefault, 3, true>`
is aliased to `Line3d`.
## Design patterns
If you plan to add a new geometric entity type,
please adopt these conventions:
+ Each geometric entity may be default-constructed.
The default constructor will initialize the state to some
valid unit-length entity, usually with some part of
its state at the origin of the coordinate system.
+ Entities may always be constructed by passing in values
for their internal state.
Alternate construction methods are declared as free functions
such as `make_CircleFrom3Points()`
+ Use template metaprogramming to make methods available
only when the template dimension gives them semantic meaning.
For example, a 2-D line segment's perpendicular bisector
is another line segment, but a 3-D line segment's perpendicular
line segment is a plane.
Note how this is accomplished and apply this pattern to
new geometric entities or new methods on existing entities.
+ Some entities may have invalid state.
If this is possible, the entity will have an `IsValid()` method.
For example, a sphere may be invalid because the user or some
construction technique specified a zero or negative radius.
+ When signed distance is semantically meaningful, provide it
in favor of or in addition to unsigned distance.
+ Accept a tolerance parameter when appropriate,
but provide a sensible default value.
You may want to perform exact arithmetic versions of tests,
but please provide fast, tolerance-based versions as well.
......@@ -35,6 +35,7 @@ set(headers
CellShape.h
CellTraits.h
Flags.h
Geometry.h
Hash.h
ImplicitFunction.h
ListTag.h
......@@ -60,11 +61,20 @@ set(headers
VecVariable.h
VirtualObjectBase.h
UnaryPredicates.h
)
)
set(template_sources
Geometry.hxx
)
vtkm_pyexpander_generated_file(Math.h)
vtkm_declare_headers(${headers})
vtkm_declare_headers(
${headers}
${template_sources}
EXCLUDE_FROM_TESTING
${template_sources}
)
#-----------------------------------------------------------------------------
#first add all the components vtkm that are shared between control and exec
......
//============================================================================
// 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.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_Geometry_h
#define vtk_m_Geometry_h
#include <vtkm/VectorAnalysis.h>
namespace vtkm
{
// Forward declarations of geometric types:
template <typename CoordType, int Dim, bool IsTwoSided>
struct Ray;
template <typename CoordType, int Dim>
struct LineSegment;
template <typename CoordType>
struct Plane;
template <typename CoordType, int Dim>
struct Sphere;
/// Represent an infinite or semi-infinite line segment with a point and a direction
///
/// The \a IsTwoSided template parameter indicates whether the class represents
/// an infinite line extending in both directions from the base point or a
/// semi-infinite ray extending only in the positive direction from the base points.
template <typename CoordType = vtkm::FloatDefault, int Dim = 3, bool IsTwoSided = false>
struct Ray
{
static constexpr int Dimension = Dim;
using Vector = vtkm::Vec<CoordType, Dim>;
static constexpr bool TwoSided = IsTwoSided;
Vector Origin;
Vector Direction; // always stored as a unit length.
/// Construct a default 2-D ray, from (0,0) pointing along the +x axis.
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
VTKM_EXEC_CONT Ray();
/// Construct a default 3-D ray from (0,0,0) pointing along the +x axis.
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 3, int>::type = 0>
VTKM_EXEC_CONT Ray();
/// Construct a ray from a point and direction.
VTKM_EXEC_CONT
Ray(const Vector& point, const Vector& direction);
/// Construct a ray from a line segment.
VTKM_EXEC_CONT
Ray(const LineSegment<CoordType, Dim>& segment);
/// Return whether the ray is valid or not.
///
/// It is possible for an invalid direction (zero length)
/// to be passed to the constructor. When this happens,
/// the constructor divides by zero, leaving Inf in all
/// components.
VTKM_EXEC_CONT
bool IsValid() const;
/// Compute a point along the line. \a param values > 0 lie on the ray.
VTKM_EXEC_CONT
Vector Evaluate(CoordType param) const;
/// Return the minmum distance from \a point to this line/ray.
///
/// Note that when the direction has zero length, this simplifies
/// the distance between \a point and the ray's origin.
/// Otherwise, the distance returned is either the perpendicular
/// distance from \a point to the line or the distance
/// to the ray's origin.
VTKM_EXEC_CONT
CoordType DistanceTo(const Vector& point) const;
/// Return the minimum distance between the ray/line and \a point.
VTKM_EXEC_CONT
CoordType DistanceTo(const Vector& point, CoordType& param, Vector& projectedPoint) const;
/// Compute the non-degenerate point where two 2-D rays intersect, or return false.
///
/// If true is returned, then the rays intersect in a unique point and
/// \a point is set to that location.
///
/// If false is returned, then either
/// (1) the rays are parallel and may be
/// coincident (intersect everywhere) or offset (intersect nowhere); or
/// (2) the lines intersect but not the rays (because the the intersection
/// occurs in the negative parameter space of one or both rays).
/// In the latter case (2), the \a point is still set to the location of the
/// intersection.
///
/// The tolerance \a tol is the minimum acceptable denominator used to
/// compute the intersection point coordinates and thus dictates the
/// maximum distance from the segments at which intersections will be
/// reported as valid.
template <bool OtherTwoSided, int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
VTKM_EXEC_CONT bool Intersect(const Ray<CoordType, Dim, OtherTwoSided>& other,
Vector& point,
CoordType tol = 0.f);
};
/// Represent a finite line segment with a pair of points
template <typename CoordType = vtkm::FloatDefault, int Dim = 3>
struct LineSegment
{
static constexpr int Dimension = Dim;
using Vector = vtkm::Vec<CoordType, Dim>;
Vector Endpoints[2];
/// Construct a default segment from (0,0) to (1,0).
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
VTKM_EXEC_CONT LineSegment();
/// Construct a default segment from (0,0,0) to (1,0,0).
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 3, int>::type = 0>
VTKM_EXEC_CONT LineSegment();
/// Construct a segment spanning points \a p0 and \a p1.
VTKM_EXEC_CONT
LineSegment(const Vector& p0, const Vector& p1);
/// Return whether this line segment has an infinitesimal extent (i.e., whether the endpoints are coincident)
VTKM_EXEC_CONT
bool IsSingular(CoordType tol2 = static_cast<CoordType>(1.0e-6f)) const;
/// Construct a plane bisecting this line segment (only when Dimension is 3).
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 3, int>::type = 0>
VTKM_EXEC_CONT Plane<CoordType> PerpendicularBisector() const;
/// Construct a perpendicular bisector to this line segment (only when Dimension is 2).
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
VTKM_EXEC_CONT Ray<CoordType, Dim, true> PerpendicularBisector() const;
/// Return the midpoint of the line segment
VTKM_EXEC_CONT
Vector Center() const { return this->Evaluate(0.5f); }
/// Return the vector pointing to endpoint 1 from endpoint 0.
/// This vector is not of unit length and in the case of
/// degenerate lines, may have zero length.
/// Call vtkm::Normal() on the return value if you want a normalized result.
VTKM_EXEC_CONT
Vector Direction() const { return this->Endpoints[1] - this->Endpoints[0]; }
/// Compute a point along the line. \a param values in [0,1] lie on the line segment.
VTKM_EXEC_CONT
Vector Evaluate(CoordType param) const;
/// Return the minmum distance from \a point to this line segment.
///
/// Note that when the endpoints are coincident, this simplifies
/// the distance between \a point and either endpoint.
/// Otherwise, the distance returned is either the perpendicular
/// distance from \a point to the line segment or the distance
/// to the nearest endpoint (whichever is smaller).
VTKM_EXEC_CONT
CoordType DistanceTo(const Vector& point) const;
/// Return the minimum distance between the line segment and \a point.
VTKM_EXEC_CONT
CoordType DistanceTo(const Vector& point, CoordType& param, Vector& projectedPoint) const;
/// Compute the non-degenerate point where two (infinite) 2-D line segments intersect, or return false.
///
/// If true is returned, then the lines intersect in a unique point and
/// \a point is set to that location.
///
/// If false is returned, then the lines are parallel and either they are
/// coincident (intersect everywhere) or offset (intersect nowhere).
///
/// The tolerance \a tol is the minimum acceptable denominator used to
/// compute the intersection point coordinates and thus dictates the
/// maximum distance from the segments at which intersections will be
/// reported as valid.
template <int Dim_ = Dim, typename std::enable_if<Dim_ == 2, int>::type = 0>
VTKM_EXEC_CONT bool IntersectInfinite(const LineSegment<CoordType, Dim>& other,
Vector& point,
CoordType tol = 0.f);
};
/// Represent a plane with a base point (origin) and normal vector.
template <typename CoordType = vtkm::FloatDefault>
struct Plane
{
using Vector = vtkm::Vec<CoordType, 3>;
Vector Origin;
Vector Normal;
/// Construct a default plane whose base point is the origin and whose normal is (0,0,1)
VTKM_EXEC_CONT
Plane();
/// Construct a plane with the given \a origin and \a normal.
VTKM_EXEC_CONT
Plane(const Vector& origin, const Vector& normal, CoordType tol2 = static_cast<CoordType>(1e-8f));
/// Return true if the plane's normal is well-defined to within the given tolerance.
VTKM_EXEC_CONT
bool IsValid() const { return !vtkm::IsInf(this->Normal[0]); }
/// Return the **signed** distance from the plane to the point.
VTKM_EXEC_CONT
CoordType DistanceTo(const Vector& point) const;
/// Return the closest point in the plane to the given point.
VTKM_EXEC_CONT
Vector ClosestPoint(const Vector& point) const;
/// Intersect this plane with the ray (or line if the ray is two-sided).
///
/// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection).
/// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the
/// entire ray/line lies in the plane).
/// In the latter case, \a lineInPlane will be true upon exit.
///
/// If this method returns true, then \a parameter will be set to a number indicating
/// where along the ray/line the plane hits and \a point will be set to that location.
/// If the input is a ray, the \a parameter will be non-negative.
template <bool IsTwoSided>
VTKM_EXEC_CONT bool Intersect(const Ray<CoordType, 3, IsTwoSided>& ray,
CoordType& parameter,
Vector& point,
bool& lineInPlane,
CoordType tol = CoordType(1e-6f)) const;
/// Intersect this plane with the line \a segment.
///
/// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection).
/// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the
/// entire line segment lies in the plane).
/// In the latter case, \a lineInPlane will be true upon exit.
///
/// If this method returns true, then \a parameter will be set to a number in [0,1] indicating
/// where along the line segment the plane hits.
VTKM_EXEC_CONT
bool Intersect(const LineSegment<CoordType>& segment,
CoordType& parameter,
bool& lineInPlane) const;
/// Intersect this plane with the line \a segment.
///
/// Returns true if there is a non-degenrate intersection (i.e., an isolated point of intersection).
/// Returns false if there is no intersection *or* if the intersection is degenerate (i.e., the
/// entire line segment lines in the plane).
/// In the latter case, \a lineInPlane will be true upon exit.
///
/// If this method returns true, then \a parameter will be set to a number in [0,1] indicating
/// where along the line segment the plane hits and \a point will be set to that location.
VTKM_EXEC_CONT
bool Intersect(const LineSegment<CoordType>& segment,
CoordType& parameter,
Vector& point,
bool& lineInPlane) const;
/// Intersect this plane with another plane.
///
/// Returns true if there is a non-degenrate intersection (i.e., a line of intersection).
/// Returns false if there is no intersection *or* if the intersection is degenerate
/// (i.e., the planes are coincident).
/// In the latter case, \a coincident will be true upon exit and \a segment will
/// unmodified.
///
/// If this method returns true, then the resulting \a segment will have its
/// base point on the line of intersection and its second point will be a unit
/// length away in the direction of the cross produce of the input plane normals
/// (this plane crossed with the \a other).
///
/// The tolerance \a tol is the minimum squared length of the cross-product
/// of the two plane normals. It is also compared to the squared distance of
/// the base point of \a other away from \a this plane when considering whether
/// the planes are coincident.
VTKM_EXEC_CONT
bool Intersect(const Plane<CoordType>& other,
Ray<CoordType, 3, true>& ray,
bool& coincident,
CoordType tol2 = static_cast<CoordType>(1e-6f)) const;
};
/// Represent a sphere of the given \a Dimension.
/// If a constructor is given an invalid specification, then
/// the Radius of the resulting sphere will be -1.
template <typename CoordType = vtkm::FloatDefault, int Dim = 3>
struct Sphere
{
static constexpr int Dimension = Dim;
using Vector = vtkm::Vec<CoordType, Dim>;
Vector Center;
CoordType Radius;
/// Construct a default sphere (unit radius at the origin).
VTKM_EXEC_CONT Sphere();
/// Construct a sphere from a center point and radius.
VTKM_EXEC_CONT Sphere(const Vector& center, CoordType radius);
/// Return true if the sphere is valid (i.e., has a strictly positive radius).
VTKM_EXEC_CONT
bool IsValid() const { return this->Radius > 0.f; }
/// Return whether the point lies strictly inside the sphere.
VTKM_EXEC_CONT
bool Contains(const Vector& point, CoordType tol2 = 0.f) const;
/// Classify a point as inside (-1), on (0), or outside (+1) of the sphere.
///
/// The tolerance \a tol2 is the maximum allowable difference in squared
/// magnitude between the squared radius and the squared distance between
/// the \a point and Center.
VTKM_EXEC_CONT
int Classify(const Vector& point, CoordType tol2 = 0.f) const;
};
// -----------------------------------------------------------------------------
// Synonyms
//
// These "using" statements aim to make it easier to use the templated
// structs above when working with a particular dimension and/or the
// default floating-point type.
/// Lines are two-sided rays:
template <typename CoordType, int Dim = 3>
using Line = Ray<CoordType, Dim, true>;
// Shortcuts for 2D and 3D rays, lines, and line segments:
template <typename CoordType>
using Ray2 = Ray<CoordType, 2>;
template <typename CoordType>
using Ray3 = Ray<CoordType, 3>;
template <typename CoordType>
using Line2 = Line<CoordType, 2>;
template <typename CoordType>
using Line3 = Line<CoordType, 3>;
template <typename CoordType>
using LineSegment2 = LineSegment<CoordType, 2>;
template <typename CoordType>
using LineSegment3 = LineSegment<CoordType, 3>;
/// Circle is an alias for a 2-Dimensional sphere.
template <typename T>
using Circle = Sphere<T, 2>;
// Aliases for d-dimensional spheres.
template <typename T>
using Sphere2 = Sphere<T, 2>;
template <typename T>
using Sphere3 = Sphere<T, 3>;
// Shortcuts for default floating-point types
using Ray2d = Ray2<vtkm::FloatDefault>;
using Ray3d = Ray3<vtkm::FloatDefault>;
using Line2d = Line2<vtkm::FloatDefault>;
using Line3d = Line3<vtkm::FloatDefault>;
using LineSegment2d = LineSegment2<vtkm::FloatDefault>;
using LineSegment3d = LineSegment3<vtkm::FloatDefault>;
using Plane3d = Plane<vtkm::FloatDefault>;
using Circle2d = Circle<vtkm::FloatDefault>;
using Sphere2d = Sphere2<vtkm::FloatDefault>;
using Sphere3d = Sphere3<vtkm::FloatDefault>;
// -----------------------------------------------------------------------------
// Construction techniques
//
// These are free functions that create instances of geometric structs by taking
// in data that is not identical to the state of the struct and converting it
// into state for the struct.
/// Construct a plane from a point plus one of: a line, a ray, or a line segment.
///
/// The plane returned will contain the point and the line/ray/segment.
/// The plane normal will be the cross product of the line/ray/segment's direction
/// and the vector from the line/ray/segment's origin to the given \a point.
/// If the \a point is collinear with the line/ray/line-segment, an invalid
/// plane will be returned.
template <typename CoordType, bool IsTwoSided>
VTKM_EXEC_CONT vtkm::Plane<CoordType> make_PlaneFromPointAndLine(
const vtkm::Vec<CoordType, 3>& point,
const vtkm::Ray<CoordType, 3, IsTwoSided>& ray,
CoordType tol2 = static_cast<CoordType>(1e-8f));
template <typename CoordType>
VTKM_EXEC_CONT vtkm::Plane<CoordType> make_PlaneFromPointAndLineSegment(
const vtkm::Vec<CoordType, 3>& point,
const vtkm::LineSegment3<CoordType>& segment,
CoordType tol2 = static_cast<CoordType>(1e-8f));
/// Construct a circle from 3 points.
template <typename CoordType>
VTKM_EXEC_CONT vtkm::Circle<CoordType> make_CircleFrom3Points(
const typename vtkm::Vec<CoordType, 2>& p0,
const typename vtkm::Vec<CoordType, 2>& p1,
const typename vtkm::Vec<CoordType, 2>& p2,
CoordType tol = static_cast<CoordType>(1e-6f));
/// Construct a sphere from 4 points.
template <typename CoordType>
VTKM_EXEC_CONT vtkm::Sphere<CoordType, 3> make_SphereFrom4Points(
const vtkm::Vec<CoordType, 3>& a0,
const vtkm::Vec<CoordType, 3>& a1,
const vtkm::Vec<CoordType, 3>& a2,
const vtkm::Vec<CoordType, 3>& a3,
CoordType tol = static_cast<CoordType>(1e-6f));
} // namespace vtkm
#include "vtkm/Geometry.hxx"
#endif // vtk_m_Geometry_h
This diff is collapsed.
......@@ -1815,6 +1815,18 @@ static inline VTKM_EXEC_CONT T Min(const T& x, const T& y)
return detail::Min(x, y, typename vtkm::TypeTraits<T>::DimensionalityTag());
}
/// Clamp \p x to the given range.
///
inline VTKM_EXEC_CONT vtkm::Float32 Clamp(vtkm::Float32 x, vtkm::Float32 lo, vtkm::Float32 hi)
{
return x > lo ? (x < hi ? x : hi) : lo;
}
inline VTKM_EXEC_CONT vtkm::Float64 Clamp(vtkm::Float64 x, vtkm::Float64 lo, vtkm::Float64 hi)
{
return x > lo ? (x < hi ? x : hi) : lo;
}
//-----------------------------------------------------------------------------
//#ifdef VTKM_CUDA
......
......@@ -637,6 +637,18 @@ static inline VTKM_EXEC_CONT T Min(const T& x, const T& y)
return detail::Min(x, y, typename vtkm::TypeTraits<T>::DimensionalityTag());
}
/// Clamp \p x to the given range.
///
inline VTKM_EXEC_CONT vtkm::Float32 Clamp(vtkm::Float32 x, vtkm::Float32 lo, vtkm::Float32 hi)
{
return x > lo ? (x < hi ? x : hi) : lo;
}
inline VTKM_EXEC_CONT vtkm::Float64 Clamp(vtkm::Float64 x, vtkm::Float64 lo, vtkm::Float64 hi)
{
return x > lo ? (x < hi ? x : hi) : lo;
}
//-----------------------------------------------------------------------------
//#ifdef VTKM_CUDA
......
......@@ -92,6 +92,19 @@ struct TypeListTagFieldVec4
{
};
/// A list containing common types for floating-point vectors. Specifically contains
/// floating point vectors of size 2, 3, and 4 with floating point components.
/// Scalars are not included.
///
struct TypeListTagFloatVec : vtkm::ListTagBase<vtkm::Vec<vtkm::Float32, 2>,
vtkm::Vec<vtkm::Float64, 2>,
vtkm::Vec<vtkm::Float32, 3>,
vtkm::Vec<vtkm::Float64, 3>,
vtkm::Vec<vtkm::Float32, 4>,
vtkm::Vec<vtkm::Float64, 4>>
{
};
/// A list containing common types for values in fields. Specifically contains
/// floating point scalars and vectors of size 2, 3, and 4 with floating point
/// components.
......
......@@ -208,6 +208,84 @@ TriangleNormal(const vtkm::Vec<T, 3>& a, const vtkm::Vec<T, 3>& b, const vtkm::V
return vtkm::Cross(b - a, c - a);
}
//-----------------------------------------------------------------------------
/// \brief Project a vector onto another vector.
///
/// This method computes the orthogonal projection of the vector v onto u;
/// that is, it projects its first argument onto its second.
///
/// Note that if the vector \a u has zero length, the output
/// vector will have all its entries equal to NaN.
template <typename T, int N>
VTKM_EXEC_CONT vtkm::Vec<T, N> Project(const vtkm::Vec<T, N>& v, const vtkm::Vec<T, N>& u)
{
T uu = vtkm::Dot(u, u);
T uv = vtkm::Dot(u, v);
T factor = uv / uu;
vtkm::Vec<T, N> result = factor * u;
return result;
}
//-----------------------------------------------------------------------------
/// \brief Project a vector onto another vector, returning only the projected distance.
///
/// This method computes the orthogonal projection of the vector v onto u;
/// that is, it projects its first argument onto its second.
///
/// Note that if the vector \a u has zero length, the output will be NaN.
template <typename T, int N>
VTKM_EXEC_CONT T ProjectedDistance(const vtkm::Vec<T, N