Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add vtkm/Geometry.h and test it.

Browse Source 
This commit adds several geometric constructs to vtk-m
in the `vtkm/Geometry.h` header. They may be used from
both the execution and control environments.

We also add methods to perform projection and Gram-Schmidt
orthonormalization to `vtkm/VectorAnalysis.h`.

See `docs/changelog/geometry.md` included in this commit
for more information.
This commit is contained in:
David Thompson 2018-06-08 00:27:30 -04:00
parent 7c25841053
commit 880d8a989e
15 changed files with 1899 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

75
docs/changelog/geometry.md Normal file
View File

@ -0,0 +1,75 @@
# 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.

14
vtkm/CMakeLists.txt
View File

@ -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

425
vtkm/Geometry.h Normal file
View File

@ -0,0 +1,425 @@
//============================================================================
// 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

647
vtkm/Geometry.hxx Normal file
View File

@ -0,0 +1,647 @@
//============================================================================
// 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.
//============================================================================
namespace vtkm
{
// -----------------------------------------------------------------------------
// Ray
template <typename CoordType, int Dim, bool IsTwoSided>
template <int Dim_, typename std::enable_if<Dim_ == 2, int>::type>
Ray<CoordType, Dim, IsTwoSided>::Ray()
: Origin{ 0.f }
, Direction{ 1.f, 0.f }
{
}
template <typename CoordType, int Dim, bool IsTwoSided>
template <int Dim_, typename std::enable_if<Dim_ == 3, int>::type>
Ray<CoordType, Dim, IsTwoSided>::Ray()
: Origin{ 0.f }
, Direction{ 1.f, 0.f, 0.f }
{
}
template <typename CoordType, int Dim, bool IsTwoSided>
Ray<CoordType, Dim, IsTwoSided>::Ray(const LineSegment<CoordType, Dim>& segment)
: Origin(segment.Endpoints[0])
, Direction(vtkm::Normal(segment.Direction()))
{
}
template <typename CoordType, int Dim, bool IsTwoSided>
Ray<CoordType, Dim, IsTwoSided>::Ray(const Vector& point, const Vector& direction)
: Origin(point)
, Direction(vtkm::Normal(direction))
{
}
template <typename CoordType, int Dim, bool IsTwoSided>
typename Ray<CoordType, Dim, IsTwoSided>::Vector Ray<CoordType, Dim, IsTwoSided>::Evaluate(
CoordType param) const
{
auto pointOnLine = this->Origin + this->Direction * param;
return pointOnLine;
}
template <typename CoordType, int Dim, bool IsTwoSided>
bool Ray<CoordType, Dim, IsTwoSided>::IsValid() const
{
// At least on Ubuntu 17.10, cuda 9.1 will fail with an internal
// compiler error when calling vtkm::IsInf() here. But the fix
// below works. The fix should be removed as soon as our dashboards
// allow it.
#if __CUDACC_VER_MAJOR__ == 9 && __CUDACC_VER_MINOR__ == 1
return !isinf(this->Direction[0]);
#else
return !vtkm::IsInf(this->Direction[0]);
#endif
}
template <typename CoordType, int Dim, bool IsTwoSided>
CoordType Ray<CoordType, Dim, IsTwoSided>::DistanceTo(const Vector& point) const
{
Vector closest;
CoordType param;
return this->DistanceTo(point, param, closest);
}
template <typename CoordType, int Dim, bool IsTwoSided>
CoordType Ray<CoordType, Dim, IsTwoSided>::DistanceTo(const Vector& point,
CoordType& param,
Vector& projectedPoint) const
{
const auto& dir = this->Direction;
auto mag2 = vtkm::MagnitudeSquared(dir);
if (mag2 <= static_cast<CoordType>(0.0f))
{
// We have a point, not a line segment.
projectedPoint = this->Origin;
param = static_cast<CoordType>(0.0f);
return vtkm::Magnitude(point - this->Origin);
}
// Find the closest point on the line, then clamp to the ray if the
// parameter value is negative.
param = vtkm::Dot(point - this->Origin, dir) / mag2;
if (!TwoSided)
{
param = vtkm::Max(param, static_cast<CoordType>(0.0f));
}
// Compute the distance between the closest point and the input point.
projectedPoint = this->Evaluate(param);
auto dist = vtkm::Magnitude(point - projectedPoint);
return dist;
}
template <typename CoordType, int Dim, bool IsTwoSided>
template <bool OtherTwoSided, int Dim_, typename std::enable_if<Dim_ == 2, int>::type>
bool Ray<CoordType, Dim, IsTwoSided>::Intersect(const Ray<CoordType, Dim, OtherTwoSided>& other,
Vector& point,
CoordType tol)
{
auto d1 = this->Direction;
auto d2 = other.Direction;
auto denom = d1[0] * d2[1] - d1[1] * d2[0];
if (vtkm::Abs(denom) < tol)
{ // The lines are coincident or at least parallel.
return false;
}
const auto& a = this->Origin;
const auto& b = other.Origin;
CoordType numerU = a[1] * d2[0] + d2[1] * b[0] - b[1] * d2[0] - d2[1] * a[0];
CoordType uParam = numerU / denom;
point = a + uParam * d1;
if (TwoSided && OtherTwoSided)
{
return true;
}
else
{
CoordType numerV = d1[0] * (a[1] - b[1]) - d1[1] * (a[0] - b[0]);
CoordType vParam = numerV / denom;
return (TwoSided || (uParam + tol) > 0) && (OtherTwoSided || (vParam + tol) > 0);
}
}
// -----------------------------------------------------------------------------
// LineSegment
template <typename CoordType, int Dim>
template <int Dim_, typename std::enable_if<Dim_ == 2, int>::type>
LineSegment<CoordType, Dim>::LineSegment()
: Endpoints{ { 0.f }, { 1.f, 0.f } }
{
}
template <typename CoordType, int Dim>
template <int Dim_, typename std::enable_if<Dim_ == 3, int>::type>
LineSegment<CoordType, Dim>::LineSegment()
: Endpoints{ { 0.f }, { 1.f, 0.f, 0.f } }
{
}
template <typename CoordType, int Dim>
LineSegment<CoordType, Dim>::LineSegment(const Vector& p0, const Vector& p1)
: Endpoints{ p0, p1 }
{
}
template <typename CoordType, int Dim>
bool LineSegment<CoordType, Dim>::IsSingular(CoordType tol2) const
{
return vtkm::MagnitudeSquared(this->Direction()) < tol2;
}
template <typename CoordType, int Dim>
template <int Dim_, typename std::enable_if<Dim_ == 2, int>::type>
Ray<CoordType, Dim, true> LineSegment<CoordType, Dim>::PerpendicularBisector() const
{
const Vector dir = this->Direction();
const Vector perp(-dir[1], dir[0]);
const Vector mid = this->Center();
return Ray<CoordType, Dim, true>(mid, perp);
}
template <typename CoordType, int Dim>
template <int Dim_, typename std::enable_if<Dim_ == 3, int>::type>
Plane<CoordType> LineSegment<CoordType, Dim>::PerpendicularBisector() const
{
return Plane<CoordType>(this->Center(), this->Direction());
}
template <typename CoordType, int Dim>
typename LineSegment<CoordType, Dim>::Vector LineSegment<CoordType, Dim>::Evaluate(
CoordType param) const
{
auto pointOnLine = this->Endpoints[0] * (1.0f - param) + this->Endpoints[1] * param;
return pointOnLine;
}
template <typename CoordType, int Dim>
CoordType LineSegment<CoordType, Dim>::DistanceTo(const Vector& point) const
{
Vector closest;
CoordType param;
return this->DistanceTo(point, param, closest);
}
template <typename CoordType, int Dim>
CoordType LineSegment<CoordType, Dim>::DistanceTo(const Vector& point,
CoordType& param,
Vector& projectedPoint) const
{
auto dir = this->Endpoints[1] - this->Endpoints[0];
auto mag2 = vtkm::MagnitudeSquared(dir);
if (mag2 <= static_cast<CoordType>(0.0f))
{
// We have a point, not a line segment.
projectedPoint = this->Endpoints[0];
param = static_cast<CoordType>(0.0f);
return vtkm::Magnitude(point - this->Endpoints[0]);
}
// Find the closest point on the line, then clamp to the line segment
param = vtkm::Clamp(vtkm::Dot(point - this->Endpoints[0], dir) / mag2,
static_cast<CoordType>(0.0f),
static_cast<CoordType>(1.0f));
// Compute the distance between the closest point and the input point.
projectedPoint = this->Evaluate(param);
auto dist = vtkm::Magnitude(point - projectedPoint);
return dist;
}
template <typename CoordType, int Dim>
template <int Dim_, typename std::enable_if<Dim_ == 2, int>::type>
bool LineSegment<CoordType, Dim>::IntersectInfinite(const LineSegment<CoordType, Dim>& other,
Vector& point,
CoordType tol)
{
auto d1 = this->Direction();
auto d2 = other.Direction();
auto denom = d1[0] * d2[1] - d1[1] * d2[0];
if (vtkm::Abs(denom) < tol)
{ // The lines are coincident or at least parallel.
return false;
}
const auto& a = this->Endpoints;
const auto& b = other.Endpoints;
CoordType numerX = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * -d2[0] -
(b[0][0] * b[1][1] - b[0][1] * b[1][0]) * -d1[0];
CoordType numerY = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) * -d2[1] -
(b[0][0] * b[1][1] - b[0][1] * b[1][0]) * -d1[1];
point = Vector(numerX / denom, numerY / denom);
return true;
}
// -----------------------------------------------------------------------------
// Plane
template <typename CoordType>
Plane<CoordType>::Plane()
: Origin{ 0.f, 0.f, 0.f }
, Normal{ 0.f, 0.f, 1.f }
{
}
template <typename CoordType>
Plane<CoordType>::Plane(const Vector& origin, const Vector& normal, CoordType tol2)
: Origin(origin)
, Normal(vtkm::Normal(normal))
{
if (tol2 > 0.0f && vtkm::MagnitudeSquared(normal) < tol2)
{
auto inf = vtkm::Infinity<CoordType>();
Normal = vtkm::Vec<CoordType, 3>(inf, inf, inf);
}
}
template <typename CoordType>
CoordType Plane<CoordType>::DistanceTo(const Vector& point) const
{
auto dist = vtkm::Dot(point - this->Origin, this->Normal);
return dist;
}
template <typename CoordType>
typename Plane<CoordType>::Vector Plane<CoordType>::ClosestPoint(const Vector& point) const
{
auto vop = vtkm::Project(point - this->Origin, this->Normal);
auto closest = point - vop;
return closest;
}
template <typename CoordType>
template <bool IsTwoSided>
bool Plane<CoordType>::Intersect(const Ray<CoordType, 3, IsTwoSided>& ray,
CoordType& parameter,
Vector& point,
bool& lineInPlane,
CoordType tol) const
{
CoordType d0 = this->DistanceTo(ray.Origin);
CoordType dirDot = vtkm::Dot(this->Normal, ray.Direction);
// If the ray/line lies parallel to the plane, the intersection is degenerate:
if (vtkm::Abs(dirDot) < tol)
{
lineInPlane = (vtkm::Abs(d0) < tol);
return false;
}
lineInPlane = false;
parameter = -d0 / dirDot;
// If we have a ray (not a line) and it points away from the
// side of the plane where its origin lies, then there is no
// intersection.
if (!IsTwoSided)
{
if (parameter < 0.0f)
{
return false;
}
}
// Check whether an endpoint lies in the plane:
if (vtkm::Abs(d0) < tol)
{
parameter = static_cast<CoordType>(0.0f);
point = ray.Origin;
return true;
}
// The perpendicular distance of the origin to the plane
// forms one side of a triangle whose hypotenuse is the
// parameter value (because ray.Direction has unit length).
// The dot product of the plane normal and ray direction
// is the cosine of the angle between the hypotenuse and
// the shortest path to the plane, so....
point = ray.Origin + parameter * ray.Direction;
return true;
}
template <typename CoordType>
bool Plane<CoordType>::Intersect(const LineSegment<CoordType>& segment,
CoordType& parameter,
bool& lineInPlane) const
{
Vector point;
return this->Intersect(segment, parameter, point, lineInPlane);
}
template <typename CoordType>
bool Plane<CoordType>::Intersect(const LineSegment<CoordType>& segment,
CoordType& parameter,
Vector& point,
bool& lineInPlane) const
{
CoordType d0 = this->DistanceTo(segment.Endpoints[0]);
CoordType d1 = this->DistanceTo(segment.Endpoints[1]);
if (d0 == 0 && d1 == 0)
{
// The entire segment lies in the plane.
lineInPlane = true;
return true;
}
lineInPlane = false;
// Check whether an endpoint lies in the plane:
if (d0 == 0)
{
parameter = static_cast<CoordType>(0.0f);
point = segment.Endpoints[0];
return true;
}
if (d1 == 0)
{
parameter = static_cast<CoordType>(1.0f);
point = segment.Endpoints[1];
return true;
}
// See whether endpoints lie on opposite sides of the plane:
//
// Note that the vtkm::SignBit comments below cause an internal compiler
// error on cuda 9.1, ubuntu 17.10 or they would be used:
bool c0 = d0 < 0; // !!vtkm::SignBit(d0);
bool c1 = d1 < 0; // !!vtkm::SignBit(d1);
CoordType a0 = vtkm::Abs(d0);
CoordType a1 = vtkm::Abs(d1);
if (c0 == c1)
{
// Both endpoints lie to the same side of the plane, so
// there is no intersection. Report the closest endpoint.
parameter = static_cast<CoordType>(a0 < a1 ? 0.0f : 1.0f);
point = segment.Endpoints[a0 < a1 ? 0 : 1];
return false;
}
// Endpoint distances have the opposite sign; there must be
// an intersection. It must occur at distance 0, and distance
// varies linearly from d0 to d1, so...
parameter = a0 / (a0 + a1);
point = segment.Endpoints[0] * (1.0f - parameter) + segment.Endpoints[1] * parameter;
return true;
}
template <typename CoordType>
bool Plane<CoordType>::Intersect(const Plane<CoordType>& other,
Ray<CoordType, 3, true>& ray,
bool& coincident,
CoordType tol2) const
{
auto dir = vtkm::Cross(this->Normal, other.Normal);
auto mag2 = vtkm::MagnitudeSquared(dir);
if (mag2 < tol2)
{ // The planes are parallel.
auto dist = this->DistanceTo(other.Origin);
coincident = dist * dist < tol2;
return false;
}
// The planes intersect. We want to find a point on the new plane
// and we want it to be near the other plane base points to avoid
// precision issues in the future.
// So, project each plane origin to the other plane along a line
// perpendicular to the plane and the output line. Both of these
// points are on the output line. Average the two points. The result
// will still be on the line and will be closer to the two base points.
auto nn = vtkm::Normal(dir);
auto moveDir01 = vtkm::Cross(this->Normal, nn);
auto moveDir02 = vtkm::Cross(other.Normal, nn);
Ray<CoordType, 3, true> bra(this->Origin, moveDir01);
Ray<CoordType, 3, true> brb(other.Origin, moveDir02);
vtkm::Vec<CoordType, 3> p0a;
vtkm::Vec<CoordType, 3> p0b;
CoordType pDummyA, pDummyB;
bool bDummyA, bDummyB;
auto tol = vtkm::Sqrt(tol2);
this->Intersect(brb, pDummyA, p0a, bDummyA, tol);
other.Intersect(bra, pDummyB, p0b, bDummyB, tol);
ray = vtkm::Ray<CoordType, 3, true>((p0a + p0b) * 0.5f, nn);
return true;
}
// -----------------------------------------------------------------------------
// Sphere
template <typename CoordType, int Dim>
Sphere<CoordType, Dim>::Sphere()
: Center{ 0.f }
, Radius(static_cast<CoordType>(1.f))
{
}
template <typename CoordType, int Dim>
Sphere<CoordType, Dim>::Sphere(const Vector& center, CoordType radius)
: Center(center)
, Radius(radius <= 0.f ? static_cast<CoordType>(-1.0f) : radius)
{
}
template <typename CoordType, int Dim>
bool Sphere<CoordType, Dim>::Contains(const Vector& point, CoordType tol2) const
{
return this->Classify(point, tol2) < 0;
}
template <typename CoordType, int Dim>
int Sphere<CoordType, Dim>::Classify(const Vector& point, CoordType tol2) const
{
if (!this->IsValid())
{
return 1; // All points are outside invalid spheres.
}
auto d2 = vtkm::MagnitudeSquared(point - this->Center);
auto r2 = this->Radius * this->Radius;
return d2 < r2 - tol2 ? -1 : (d2 > r2 + tol2 ? +1 : 0);
}
// -----------------------------------------------------------------------------
// Construction techniques
template <typename CoordType, bool IsTwoSided>
vtkm::Plane<CoordType> make_PlaneFromPointAndLine(const vtkm::Vec<CoordType, 3>& point,
const vtkm::Ray<CoordType, 3, IsTwoSided>& ray,
CoordType tol2)
{
auto tmpDir = point - ray.Origin;
return vtkm::Plane<CoordType>(point, vtkm::Cross(ray.Direction, tmpDir), tol2);
}
template <typename CoordType>
vtkm::Plane<CoordType> make_PlaneFromPointAndLineSegment(
const vtkm::Vec<CoordType, 3>& point,
const vtkm::LineSegment3<CoordType>& segment,
CoordType tol2)
{
auto tmpDir = point - segment.Endpoints[0];
return vtkm::Plane<CoordType>(point, vtkm::Cross(segment.Direction(), tmpDir), tol2);
}
template <typename CoordType>
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)
{
constexpr int Dim = 2;
using Vector = typename vtkm::Circle<CoordType>::Vector;
vtkm::LineSegment<CoordType, Dim> l01(p0, p1);
vtkm::LineSegment<CoordType, Dim> l02(p0, p2);
auto pb01 = l01.PerpendicularBisector();
auto pb02 = l02.PerpendicularBisector();
Vector center;
CoordType radius;
if (!pb01.IsValid() || !pb02.IsValid())
{
center = Vector(0.f, 0.f);
radius = -1.f;
return vtkm::Circle<CoordType>(center, radius);
}
auto isValid = pb01.Intersect(pb02, center, tol);
radius = isValid ? vtkm::Magnitude(center - p0) : static_cast<CoordType>(-1.0f);
if (!isValid)
{ // Initialize center to something.
center = Vector(vtkm::Nan<CoordType>(), vtkm::Nan<CoordType>());
}
return vtkm::Circle<CoordType>(center, radius);
}
template <typename CoordType>
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)
{
// Choose p3 such that the min(p3 - p[012]) is larger than any other choice of p3.
// From: http://steve.hollasch.net/cgindex/geometry/sphere4pts.html,
// retrieved 2018-06-12 and paraphrased in terms of local variable names:
//
// If circlePointInPlaneOfP3-p3 is much smaller than
// circlePointInPlaneOfP3-circleCenterWorld, then the
// sphere center will be very close to circleCenterWorld
// and subject to error.
// It's best to choose p3 so that the least of p0-p3, p1-p3, and
// p2-p3 is larger than for any other.
CoordType d0 = vtkm::MagnitudeSquared(a1 - a0);
CoordType d1 = vtkm::MagnitudeSquared(a2 - a0);
CoordType d2 = vtkm::MagnitudeSquared(a3 - a0);
CoordType d3 = vtkm::MagnitudeSquared(a2 - a1);
CoordType d4 = vtkm::MagnitudeSquared(a3 - a1);
CoordType d5 = vtkm::MagnitudeSquared(a3 - a2);
CoordType sel0 = vtkm::Min(d0, vtkm::Min(d1, d2));
CoordType sel1 = vtkm::Min(d0, vtkm::Min(d3, d4));
CoordType sel2 = vtkm::Min(d1, vtkm::Min(d3, d5));
CoordType sel3 = vtkm::Min(d2, vtkm::Min(d4, d5));
CoordType selm = vtkm::Max(vtkm::Max(sel0, sel1), vtkm::Max(sel2, sel3));
vtkm::Vec<CoordType, 3> p0 = a0;
vtkm::Vec<CoordType, 3> p1 = a1;
vtkm::Vec<CoordType, 3> p2 = a2;
vtkm::Vec<CoordType, 3> p3 = a3;
if (sel0 == selm)
{
p3 = a0;
p0 = a3;
}
else if (sel1 == selm)
{
p3 = a1;
p1 = a3;
}
else if (sel2 == selm)
{
p3 = a2;
p2 = a3;
}
else // sel3 == selm
{
// do nothing.
}
vtkm::Vec<vtkm::Vec<CoordType, 3>, 3> axes;
axes[1] = p1 - p0;
axes[2] = p2 - p0;
axes[0] = vtkm::Cross(axes[1], axes[2]);
vtkm::Vec<vtkm::Vec<CoordType, 3>, 3> basis;
int rank = vtkm::Orthonormalize(axes, basis, tol);
if (rank < 3)
{
vtkm::Sphere<CoordType, 3> invalid;
invalid.Radius = -1.0f;
return invalid;
}
// Project points to plane and fit a circle.
auto p0_p = vtkm::Vec<CoordType, 2>{ 0.f }; // This is p0's new coordinate...
auto p1_p = vtkm::Vec<CoordType, 2>(vtkm::ProjectedDistance(axes[1], basis[1]),
vtkm::ProjectedDistance(axes[1], basis[2]));
auto p2_p = vtkm::Vec<CoordType, 2>(vtkm::ProjectedDistance(axes[2], basis[1]),
vtkm::ProjectedDistance(axes[2], basis[2]));
auto circle = make_CircleFrom3Points(p0_p, p1_p, p2_p);
if (!circle.IsValid())
{
vtkm::Sphere<CoordType, 3> invalid;
invalid.Radius = -1.0f;
return invalid;
}
vtkm::Vec<CoordType, 3> circleCenterWorld =
p0 + circle.Center[0] * basis[1] + circle.Center[1] * basis[2];
vtkm::Line3<CoordType> centerRay(circleCenterWorld, basis[0]);
// If our remaining unused point (p3) lines on centerRay,
// use one of the other points to locate the sphere's center:
vtkm::Vec<CoordType, 3> circlePointInPlaneOfP3;
if (vtkm::Abs(centerRay.DistanceTo(p3)) < tol)
{
circlePointInPlaneOfP3 = p0;
}
else
{
Plane<CoordType> pp3(circleCenterWorld, basis[0]);
circlePointInPlaneOfP3 =
circleCenterWorld + vtkm::Normal(pp3.ClosestPoint(p3) - circleCenterWorld) * circle.Radius;
}
auto bisectorPlane =
vtkm::LineSegment3<CoordType>(circlePointInPlaneOfP3, p3).PerpendicularBisector();
vtkm::Vec<CoordType, 3> sphereCenter;
CoordType param;
bool lineInPlane;
if (!bisectorPlane.Intersect(centerRay, param, sphereCenter, lineInPlane, tol))
{
vtkm::Sphere<CoordType, 3> invalid;
invalid.Radius = -1.0f;
return invalid;
}
auto sphereRadius = vtkm::Magnitude(sphereCenter - p3);
vtkm::Sphere<CoordType, 3> result(sphereCenter, sphereRadius);
;
return result;
}
} // namespace vtkm

12
vtkm/Math.h
View File

@ -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

12
vtkm/Math.h.in
View File

@ -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

13
vtkm/TypeListTag.h
View File

@ -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.

78
vtkm/VectorAnalysis.h
View File

@ -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>& v, const vtkm::Vec<T, N>& u)
{
T uu = vtkm::Dot(u, u);
T uv = vtkm::Dot(u, v);
T factor = uv / uu;
return factor;
}
//-----------------------------------------------------------------------------
/// \brief Perform Gram-Schmidt orthonormalization for 3-D vectors.
///
/// See https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process for details.
/// The first output vector will always be parallel to the first input vector.
/// The remaining output vectors will be orthogonal and unit length and have
/// the same handedness as their corresponding input vectors.
///
/// This method is geometric.
/// It does not require a matrix solver.
/// However, unlike the algebraic eigensolver techniques which do use matrix
/// inversion, this method may return zero-length output vectors if some input
/// vectors are collinear. The number of non-zero (to within the specified
/// tolerance, \a tol ) output vectors is the return value.
template <typename T, int N>
VTKM_EXEC_CONT int Orthonormalize(const vtkm::Vec<vtkm::Vec<T, N>, N>& inputs,
vtkm::Vec<vtkm::Vec<T, N>, N>& outputs,
T tol = static_cast<T>(1e-6))
{
int j = 0; // j is the number of non-zero-length, non-collinear inputs encountered.
vtkm::Vec<vtkm::Vec<T, N>, N> u;
for (int i = 0; i < N; ++i)
{
u[j] = inputs[i];
for (int k = 0; k < j; ++k)
{
u[j] -= vtkm::Project(inputs[i], u[k]);
}
T rmag = vtkm::RMagnitude(u[j]);
if (rmag * tol > 1.0)
{
// skip this vector, it is zero-length or collinear with others.
continue;
}
outputs[j] = rmag * u[j];
++j;
}
for (int i = j; i < N; ++i)
{
outputs[j] = Vec<T, N>{ 0. };
}
return j;
}
} // namespace vtkm
#endif //vtk_m_VectorAnalysis_h

1
vtkm/cont/cuda/testing/CMakeLists.txt
View File

@ -28,6 +28,7 @@ set(unit_tests
UnitTestCudaDataSetExplicit.cu
UnitTestCudaDataSetSingleType.cu
UnitTestCudaDeviceAdapter.cu
UnitTestCudaGeometry.cu
UnitTestCudaImplicitFunction.cu
UnitTestCudaMath.cu
UnitTestCudaShareUserProvidedManagedMemory.cu

38
vtkm/cont/cuda/testing/UnitTestCudaGeometry.cu Normal file
View File

@ -0,0 +1,38 @@
//=============================================================================
//
// 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 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 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.
//
//=============================================================================
// Make sure that the tested code is using the device adapter specified. This
// is important in the long run so we don't, for example, use the CUDA device
// for a part of an operation where the TBB device was specified.
#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR
#include <vtkm/cont/RuntimeDeviceTracker.h>
#include <vtkm/cont/cuda/DeviceAdapterCuda.h>
#include <vtkm/testing/TestingGeometry.h>
int UnitTestCudaGeometry(int, char* [])
{
auto tracker = vtkm::cont::GetGlobalRuntimeDeviceTracker();
tracker.ForceDevice(vtkm::cont::DeviceAdapterTagCuda{});
return vtkm::cont::testing::Testing::Run(
UnitTestGeometryNamespace::RunGeometryTests<vtkm::cont::DeviceAdapterTagCuda>);
}

1
vtkm/cont/serial/testing/CMakeLists.txt
View File

@ -28,6 +28,7 @@ set(unit_tests
UnitTestSerialDataSetExplicit.cxx
UnitTestSerialDataSetSingleType.cxx
UnitTestSerialDeviceAdapter.cxx
UnitTestSerialGeometry.cxx
UnitTestSerialImplicitFunction.cxx
UnitTestSerialPointLocatorUniformGrid.cxx
UnitTestSerialVirtualObjectHandle.cxx

38
vtkm/cont/serial/testing/UnitTestSerialGeometry.cxx Normal file
View File

@ -0,0 +1,38 @@
//=============================================================================
//
// 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 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 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.
//
//=============================================================================
// Make sure that the tested code is using the device adapter specified. This
// is important in the long run so we don't, for example, use the CUDA device
// for a part of an operation where the TBB device was specified.
#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_ERROR
#include <vtkm/cont/RuntimeDeviceTracker.h>
#include <vtkm/cont/serial/DeviceAdapterSerial.h>
#include <vtkm/testing/TestingGeometry.h>
int UnitTestSerialGeometry(int, char* [])
{
auto tracker = vtkm::cont::GetGlobalRuntimeDeviceTracker();
tracker.ForceDevice(vtkm::cont::DeviceAdapterTagSerial{});
return vtkm::cont::testing::Testing::Run(
UnitTestGeometryNamespace::RunGeometryTests<vtkm::cont::DeviceAdapterTagSerial>);
}

1
vtkm/testing/CMakeLists.txt
View File

@ -21,6 +21,7 @@
set(headers
Testing.h
TestingMath.h
TestingGeometry.h
OptionParser.h
VecTraitsTests.h
)

498
vtkm/testing/TestingGeometry.h Normal file
View File

@ -0,0 +1,498 @@
//=============================================================================
//
// 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 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 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_testing_TestingGeometry_h
#define vtk_m_testing_TestingGeometry_h
#include <vtkm/Geometry.h>
#include <vtkm/Math.h>
#include <vtkm/TypeListTag.h>
#include <vtkm/VecTraits.h>
#include <vtkm/exec/FunctorBase.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/testing/Testing.h>
#define VTKM_MATH_ASSERT(condition, message) \
if (!(condition)) \
{ \
this->RaiseError(message); \
}
//-----------------------------------------------------------------------------
namespace UnitTestGeometryNamespace
{
class Coords
{
public:
static constexpr vtkm::IdComponent NUM_COORDS = 5;
template <typename T>
VTKM_EXEC_CONT vtkm::Vec<T, 3> EndpointList(vtkm::Int32 i) const
{
vtkm::Float64 coords[NUM_COORDS][3] = {
{ 1.0, 0.0, 0.0 }, { 0.0, 1.0, 0.0 }, { -1.0, 0.0, 0.0 },
{ -2.0, 0.0, 0.0 }, { 0.0, -2.0, 0.0 },
};
vtkm::Int32 j = i % NUM_COORDS;
return vtkm::Vec<T, 3>(
static_cast<T>(coords[j][0]), static_cast<T>(coords[j][1]), static_cast<T>(coords[j][2]));
}
template <typename T>
VTKM_EXEC_CONT vtkm::Vec<T, 3> ClosestToOriginList(vtkm::Int32 i) const
{
vtkm::Float64 coords[NUM_COORDS][3] = {
{ 0.5, 0.5, 0.0 }, { -0.5, 0.5, 0.0 }, { -1.0, 0.0, 0.0 },
{ -1.0, -1.0, 0.0 }, { 0.8, -0.4, 0.0 },
};
vtkm::Int32 j = i % NUM_COORDS;
return vtkm::Vec<T, 3>(
static_cast<T>(coords[j][0]), static_cast<T>(coords[j][1]), static_cast<T>(coords[j][2]));
}
template <typename T>
VTKM_EXEC_CONT T DistanceToOriginList(vtkm::Int32 i) const
{
vtkm::Float64 coords[NUM_COORDS] = {
0.707107, 0.707107, 1.0, 1.41421, 0.894427,
};
vtkm::Int32 j = i % NUM_COORDS;
return static_cast<T>(coords[j]);
}
};
//-----------------------------------------------------------------------------
template <typename T>
struct RayTests : public vtkm::exec::FunctorBase
{
VTKM_EXEC
void operator()(vtkm::Id) const
{
{
using V2 = vtkm::Vec<T, 2>;
using Ray2 = vtkm::Ray2<T>;
Ray2 ray0;
VTKM_MATH_ASSERT(test_equal(ray0.Origin, V2(0., 0.)), "Bad origin for default 2D ray ctor.");
VTKM_MATH_ASSERT(test_equal(ray0.Direction, V2(1., 0.)),
"Bad direction for default 2D ray ctor.");
// Test intersection
Ray2 ray1(V2(-1., 0.), V2(+1., +1.));
Ray2 ray2(V2(+1., 0.), V2(-1., +1.));
V2 point;
bool didIntersect = ray1.Intersect(ray2, point);
VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Ray-pair 1 should intersect.");
VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 1 should intersect at (0,1).");
// Test non intersection
Ray2 ray3(V2(-1., 0.), V2(-1., -1.));
Ray2 ray4(V2(+1., 0.), V2(+1., -1.));
didIntersect = ray1.Intersect(ray4, point);
VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Ray-pair 2 should not intersect.");
VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 2 should intersect at (0,1).");
didIntersect = ray3.Intersect(ray2, point);
VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Ray-pair 3 should not intersect.");
VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 3 should intersect at (0,1).");
didIntersect = ray3.Intersect(ray4, point);
VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Ray-pair 4 should not intersect.");
VTKM_MATH_ASSERT(test_equal(point, V2(0., 1.)), "Ray-pair 4 should intersect at (0,1).");
}
{
using V3 = vtkm::Vec<T, 3>;
vtkm::Ray<T, 3> ray0;
VTKM_MATH_ASSERT(test_equal(ray0.Origin, V3(0., 0., 0.)),
"Bad origin for default 3D ray ctor.");
VTKM_MATH_ASSERT(test_equal(ray0.Direction, V3(1., 0., 0.)),
"Bad direction for default 3D ray ctor.");
}
}
};
template <typename Device>
struct TryRayTests
{
template <typename T>
void operator()(const T&) const
{
vtkm::cont::DeviceAdapterAlgorithm<Device>::Schedule(RayTests<T>(), 1);
}
};
//-----------------------------------------------------------------------------
template <typename T>
struct LineSegmentTests : public vtkm::exec::FunctorBase
{
VTKM_EXEC
void operator()(vtkm::Id) const
{
{
using V2 = vtkm::Vec<T, 2>;
using Line2 = vtkm::Line2<T>;
vtkm::LineSegment<T, 2> seg0;
VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[0], V2(0., 0.)),
"Bad origin for default 2D line segment ctor.");
VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[1], V2(1., 0.)),
"Bad direction for default 2D line segment ctor.");
V2 p0(1., 1.);
V2 p1(3., 3.);
V2 p2(2., 2.);
// V2 p3(static_cast<T>(1.2928932), static_cast<T>(2.7071068));
V2 dir(static_cast<T>(-0.7071068), static_cast<T>(0.7071068));
vtkm::LineSegment<T, 2> seg1(p0, p1);
Line2 ray = seg1.PerpendicularBisector();
VTKM_MATH_ASSERT(test_equal(ray.Origin, p2), "Perpendicular bisector origin failed in 2D.");
VTKM_MATH_ASSERT(test_equal(ray.Direction, dir),
"Perpendicular bisector direction failed in 2D.");
}
{
using V3 = vtkm::Vec<T, 3>;
vtkm::LineSegment<T, 3> seg0;
VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[0], V3(0., 0., 0.)),
"Bad origin for default 3D line segment ctor.");
VTKM_MATH_ASSERT(test_equal(seg0.Endpoints[1], V3(1., 0., 0.)),
"Bad direction for default 3D line segment ctor.");
V3 p0(1., 1., 0.);
V3 p1(3., 3., 0.);
V3 p2(2., 2., 0.);
V3 p3(static_cast<T>(0.70710678), static_cast<T>(0.70710678), 0.);
vtkm::LineSegment<T, 3> seg1(p0, p1);
vtkm::Plane<T> bisector = seg1.PerpendicularBisector();
VTKM_MATH_ASSERT(test_equal(bisector.Origin, p2),
"Perpendicular bisector origin failed in 3D.");
VTKM_MATH_ASSERT(test_equal(bisector.Normal, p3),
"Perpendicular bisector direction failed in 3D.");
}
vtkm::Vec<T, 3> origin(0., 0., 0.);
for (vtkm::IdComponent index = 0; index < Coords::NUM_COORDS; ++index)
{
auto p0 = Coords{}.EndpointList<T>(index);
auto p1 = Coords{}.EndpointList<T>((index + 1) % Coords::NUM_COORDS);
vtkm::LineSegment<T> segment(p0, p1);
vtkm::Vec<T, 3> closest;
T param;
auto dp0 = segment.DistanceTo(p0);
auto dp1 = segment.DistanceTo(p1, param, closest);
VTKM_MATH_ASSERT(test_equal(dp0, 0.0), "Distance to endpoint 0 not zero.");
VTKM_MATH_ASSERT(test_equal(dp1, 0.0), "Distance to endpoint 1 not zero.");
VTKM_MATH_ASSERT(test_equal(param, 1.0), "Parameter value of endpoint 1 not 1.0.");
VTKM_MATH_ASSERT(test_equal(p1, closest), "Closest point not endpoint 1.");
closest = segment.Evaluate(static_cast<T>(0.0));
VTKM_MATH_ASSERT(test_equal(p0, closest), "Evaluated point not endpoint 0.");
auto dpo = segment.DistanceTo(origin, param, closest);
auto clo = Coords{}.ClosestToOriginList<T>(index);
auto dst = Coords{}.DistanceToOriginList<T>(index);
VTKM_MATH_ASSERT(test_equal(closest, clo), "Closest point to origin doesn't match.");
VTKM_MATH_ASSERT(test_equal(dpo, dst), "Distance to origin doesn't match.");
}
}
};
template <typename Device>
struct TryLineSegmentTests
{
template <typename T>
void operator()(const T&) const
{
vtkm::cont::DeviceAdapterAlgorithm<Device>::Schedule(LineSegmentTests<T>(), 1);
}
};
//-----------------------------------------------------------------------------
template <typename T>
struct PlaneTests : public vtkm::exec::FunctorBase
{
VTKM_EXEC
void operator()(vtkm::Id) const
{
vtkm::Vec<T, 3> origin(0., 0., 0.);
vtkm::Vec<T, 3> zvectr(0., 0., 5.); // intentionally not unit length to test normalization.
vtkm::Plane<T> plane;
vtkm::LineSegment<T> segment;
T dist;
bool didIntersect;
bool isLineInPlane;
vtkm::Vec<T, 3> nearest;
vtkm::Vec<T, 3> p0;
vtkm::Vec<T, 3> p1;
T param;
// Test signed plane-point distance
plane = vtkm::Plane<T>(origin, zvectr);
dist = plane.DistanceTo(vtkm::Vec<T, 3>(82., 0.5, 1.25));
VTKM_MATH_ASSERT(test_equal(dist, 1.25), "Bad positive point-plane distance.");
dist = plane.DistanceTo(vtkm::Vec<T, 3>(82., 0.5, -1.25));
VTKM_MATH_ASSERT(test_equal(dist, -1.25), "Bad negative point-plane distance.");
dist = plane.DistanceTo(vtkm::Vec<T, 3>(82., 0.5, 0.0));
VTKM_MATH_ASSERT(test_equal(dist, 0.0), "Bad zero point-plane distance.");
// Test line intersection
{
// Case 1. No intersection
segment = vtkm::LineSegment<T>((p0 = vtkm::Vec<T, 3>(1., 1., 1.)),
(p1 = vtkm::Vec<T, 3>(2., 2., 2.)));
didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane);
VTKM_MATH_ASSERT(test_equal(didIntersect, false), "Plane and line should not intersect (1).");
VTKM_MATH_ASSERT(test_equal(isLineInPlane, false),
"Line improperly reported as in plane (1).");
VTKM_MATH_ASSERT(test_equal(nearest, p0), "Unexpected nearest point (1).");
VTKM_MATH_ASSERT(test_equal(param, 0.0), "Unexpected nearest parameter value (1).");
// Case 2. Degenerate intersection (entire segment lies in plane)
segment = vtkm::LineSegment<T>((p0 = vtkm::Vec<T, 3>(1., 1., 0.)),
(p1 = vtkm::Vec<T, 3>(2., 2., 0.)));
didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane);
VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (2).");
VTKM_MATH_ASSERT(test_equal(isLineInPlane, true),
"Line improperly reported as out of plane (2).");
// Case 3. Endpoint intersection
segment = vtkm::LineSegment<T>((p0 = vtkm::Vec<T, 3>(1., 1., 1.)),
(p1 = vtkm::Vec<T, 3>(2., 2., 0.)));
didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane);
VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (3a).");
VTKM_MATH_ASSERT(test_equal(isLineInPlane, false),
"Line improperly reported as in plane (3a).");
VTKM_MATH_ASSERT(test_equal(param, 1.0), "Invalid parameter for intersection point (3a).");
VTKM_MATH_ASSERT(test_equal(nearest, p1), "Invalid intersection point (3a).");
segment = vtkm::LineSegment<T>((p0 = vtkm::Vec<T, 3>(1., 1., 0.)),
(p1 = vtkm::Vec<T, 3>(2., 2., 1.)));
didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane);
VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (3b).");
VTKM_MATH_ASSERT(test_equal(isLineInPlane, false),
"Line improperly reported as in plane (3b).");
VTKM_MATH_ASSERT(test_equal(param, 0.0), "Invalid parameter for intersection point (3b).");
VTKM_MATH_ASSERT(test_equal(nearest, p0), "Invalid intersection point (3b).");
// Case 4. General-position intersection
segment = vtkm::LineSegment<T>((p0 = vtkm::Vec<T, 3>(-1., -1., -1.)),
(p1 = vtkm::Vec<T, 3>(2., 2., 1.)));
didIntersect = plane.Intersect(segment, param, nearest, isLineInPlane);
VTKM_MATH_ASSERT(test_equal(didIntersect, true), "Plane and line should intersect (4).");
VTKM_MATH_ASSERT(test_equal(isLineInPlane, false),
"Line improperly reported as in plane (4).");
VTKM_MATH_ASSERT(test_equal(param, 0.5), "Invalid parameter for intersection point (4).");
VTKM_MATH_ASSERT(test_equal(nearest, vtkm::Vec<T, 3>(0.5, 0.5, 0)),
"Invalid intersection point (4).");
}
// Test plane-plane intersection
{
using V3 = vtkm::Vec<T, 3>;
using PlaneType = vtkm::Plane<T>;
// Case 1. Coincident planes
p0 = V3(1., 2., 3.);
p1 = V3(5., 7., -6.);
V3 nn = vtkm::Normal(V3(1., 1., 1));
PlaneType pa(p0, nn);
PlaneType pb(p1, nn);
vtkm::Line3<T> ii;
bool coincident;
didIntersect = pa.Intersect(pb, ii, coincident);
VTKM_MATH_ASSERT(test_equal(didIntersect, false),
"Coincident planes should have degenerate intersection.");
VTKM_MATH_ASSERT(test_equal(coincident, true),
"Coincident planes should be marked coincident.");
// Case 2. Offset planes
p1 = V3(5., 6., 7.);
pb = PlaneType(p1, nn);
didIntersect = pa.Intersect(pb, ii, coincident);
VTKM_MATH_ASSERT(test_equal(didIntersect, false),
"Offset planes should have degenerate intersection.");
VTKM_MATH_ASSERT(test_equal(coincident, false),
"Offset planes should not be marked coincident.");
// Case 3. General position
p1 = V3(1., 2., 0.);
V3 n2(0., 0., 1.);
pb = PlaneType(p1, n2);
didIntersect = pa.Intersect(pb, ii, coincident);
VTKM_MATH_ASSERT(test_equal(didIntersect, true),
"Proper planes should have non-degenerate intersection.");
VTKM_MATH_ASSERT(test_equal(coincident, false),
"Proper planes should not be marked coincident.");
VTKM_MATH_ASSERT(test_equal(ii.Origin, V3(2.5, 3.5, 0)),
"Unexpected intersection-line base point.");
VTKM_MATH_ASSERT(test_equal(ii.Direction, vtkm::Normal(V3(1, -1, 0))),
"Unexpected intersection-line direction.");
}
}
};
template <typename Device>
struct TryPlaneTests
{
template <typename T>
void operator()(const T&) const
{
vtkm::cont::DeviceAdapterAlgorithm<Device>::Schedule(PlaneTests<T>(), 1);
}
};
//-----------------------------------------------------------------------------
template <typename T>
struct SphereTests : public vtkm::exec::FunctorBase
{
VTKM_EXEC
void operator()(vtkm::Id) const
{
{
using V2 = vtkm::Vec<T, 2>;
V2 origin(0., 0.);
vtkm::Sphere<T, 2> defaultSphere;
VTKM_MATH_ASSERT(test_equal(defaultSphere.Center, origin), "Default circle not at origin.");
VTKM_MATH_ASSERT(test_equal(defaultSphere.Radius, 1.0), "Default circle not unit radius.");
vtkm::Sphere<T, 2> sphere(origin, -2.);
VTKM_MATH_ASSERT(test_equal(sphere.Radius, -1.0), "Negative radius should be reset to -1.");
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), false),
"Negative radius should leave sphere invalid.");
sphere = vtkm::Circle<T>(origin, 1.0);
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Circle assignment failed.");
VTKM_MATH_ASSERT(test_equal(sphere.Contains(origin), true),
"Circle does not contain its center.");
VTKM_MATH_ASSERT(test_equal(sphere.Classify(V2(1., 0.)), 0), "Circle point not on boundary.");
VTKM_MATH_ASSERT(test_equal(sphere.Classify(V2(0.75, 0.75)), +1),
"Circle contains a point that should be outside.");
V2 p0(static_cast<T>(-0.7071), static_cast<T>(-0.7071));
V2 p1(static_cast<T>(+0.7071), static_cast<T>(-0.7071));
V2 p2(static_cast<T>(0.0), static_cast<T>(1.0));
sphere = make_CircleFrom3Points(p0, p1, p2);
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Could not create 3-point circle.");
V2 p3(1, 1);
V2 p4(3, 4);
V2 p5(5, 12);
sphere = make_CircleFrom3Points(p3, p4, p5);
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Could not create 3-point circle.");
T tol = static_cast<T>(1e-3); // Use a loose tolerance
VTKM_MATH_ASSERT(test_equal(sphere.Center, vtkm::Vec<T, 2>(-12.4f, 12.1f)),
"Invalid circle center.");
VTKM_MATH_ASSERT(test_equal(sphere.Radius, static_cast<T>(17.400291f)),
"Invalid circle radius.");
VTKM_MATH_ASSERT(test_equal(sphere.Classify(p3, tol), 0),
"Generator p3 not on circle boundary.");
VTKM_MATH_ASSERT(test_equal(sphere.Classify(p4, tol), 0),
"Generator p4 not on circle boundary.");
VTKM_MATH_ASSERT(test_equal(sphere.Classify(p5, tol), 0),
"Generator p5 not on circle boundary.");
V2 p6(1, 1);
V2 p7(4, 4);
V2 p8(5, 5);
sphere = make_CircleFrom3Points(p6, p7, p8);
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), false),
"3-point circle construction should fail with points on a line.");
}
{
using V3 = vtkm::Vec<T, 3>;
V3 p0(0., 1., 0.);
V3 p1(1., 0., 0.);
V3 p2(-1., 0., 0.);
V3 p3(0., 0., 1.);
V3 p4 = vtkm::Normal(V3(1., 1., 1.));
V3 origin(0., 0., 0.);
vtkm::Sphere<T, 3> defaultSphere;
VTKM_MATH_ASSERT(test_equal(defaultSphere.Center, origin), "Default sphere not at origin.");
VTKM_MATH_ASSERT(test_equal(defaultSphere.Radius, 1.0), "Default sphere not unit radius.");
vtkm::Sphere<T, 3> sphere = make_SphereFrom4Points(p0, p1, p2, p3, static_cast<T>(1.0e-6));
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Easy sphere 1 not valid.");
VTKM_MATH_ASSERT(test_equal(sphere.Center, origin), "Easy sphere 1 not at origin.");
VTKM_MATH_ASSERT(test_equal(sphere.Radius, 1.0), "Easy sphere 1 not unit radius.");
sphere = make_SphereFrom4Points(p0, p1, p2, p4, static_cast<T>(1.0e-6));
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Easy sphere 2 not valid.");
VTKM_MATH_ASSERT(test_equal(sphere.Center, origin), "Easy sphere 2 not at origin.");
VTKM_MATH_ASSERT(test_equal(sphere.Radius, 1.0), "Easy sphere 2 not unit radius.");
V3 fancyCenter(1, 2, 3);
T fancyRadius(2.5);
V3 fp0 = fancyCenter + fancyRadius * p0;
V3 fp1 = fancyCenter + fancyRadius * p1;
V3 fp2 = fancyCenter + fancyRadius * p2;
V3 fp4 = fancyCenter + fancyRadius * p4;
sphere = make_SphereFrom4Points(fp0, fp1, fp2, fp4, static_cast<T>(1.0e-6));
VTKM_MATH_ASSERT(test_equal(sphere.IsValid(), true), "Medium sphere 1 not valid.");
VTKM_MATH_ASSERT(test_equal(sphere.Center, fancyCenter), "Medium sphere 1 not at (1,2,3).");
VTKM_MATH_ASSERT(test_equal(sphere.Radius, fancyRadius), "Medium sphere 1 not radius 2.5.");
}
}
};
template <typename Device>
struct TrySphereTests
{
template <typename T>
void operator()(const T&) const
{
vtkm::cont::DeviceAdapterAlgorithm<Device>::Schedule(SphereTests<T>(), 1);
}
};
//-----------------------------------------------------------------------------
template <typename Device>
void RunGeometryTests()
{
std::cout << "Tests for rays." << std::endl;
vtkm::testing::Testing::TryTypes(TryRayTests<Device>(), vtkm::TypeListTagFieldScalar());
std::cout << "Tests for line segments." << std::endl;
vtkm::testing::Testing::TryTypes(TryLineSegmentTests<Device>(), vtkm::TypeListTagFieldScalar());
std::cout << "Tests for planes." << std::endl;
vtkm::testing::Testing::TryTypes(TryPlaneTests<Device>(), vtkm::TypeListTagFieldScalar());
std::cout << "Tests for spheres." << std::endl;
vtkm::testing::Testing::TryTypes(TrySphereTests<Device>(), vtkm::TypeListTagFieldScalar());
}
} // namespace UnitTestGeometryNamespace
#endif //vtk_m_testing_TestingGeometry_h

48
vtkm/testing/UnitTestVectorAnalysis.cxx
View File

@ -158,6 +158,21 @@ void TestCross(const vtkm::Vec<T, 3>& x, const vtkm::Vec<T, 3>& y)
"Triangle normal is not really normal.");
}
template <typename VectorBasisType>
void TestOrthonormalize(const VectorBasisType& inputs, int expectedRank)
{
VectorBasisType outputs;
int actualRank = vtkm::Orthonormalize(inputs, outputs);
std::cout << "Testing orthonormalize\n"
<< " Rank " << actualRank << " expected " << expectedRank << "\n"
<< " Basis vectors:\n";
for (int i = 0; i < actualRank; ++i)
{
std::cout << " " << i << " " << outputs[i] << "\n";
}
VTKM_TEST_ASSERT(test_equal(actualRank, expectedRank), "Orthonormalized rank is unexpected.");
}
struct TestLinearFunctor
{
template <typename T>
@ -208,10 +223,43 @@ struct TestCrossFunctor
}
};
struct TestVectorFunctor
{
template <typename VectorType>
void operator()(const VectorType&) const
{
constexpr int NUM_COMPONENTS = VectorType::NUM_COMPONENTS;
using Traits = vtkm::VecTraits<VectorType>;
using ComponentType = typename Traits::ComponentType;
vtkm::Vec<VectorType, NUM_COMPONENTS> basis;
VectorType normalizedVector = VectorType(vtkm::RSqrt(ComponentType(NUM_COMPONENTS)));
VectorType zeroVector = VectorType(ComponentType(0));
// Test with a degenerate set of inputs:
basis[0] = zeroVector;
basis[1] = normalizedVector;
for (int ii = 2; ii < NUM_COMPONENTS; ++ii)
{
basis[ii] = zeroVector;
}
TestOrthonormalize(basis, 1);
// Now test with a valid set of inputs:
for (int ii = 0; ii < NUM_COMPONENTS; ++ii)
{
for (int jj = 0; jj < NUM_COMPONENTS; ++jj)
{
basis[ii][jj] =
ComponentType(jj == ii ? 1.0 : 0.0) + ComponentType(0.05) * ComponentType(jj);
}
}
TestOrthonormalize(basis, NUM_COMPONENTS);
}
};
void TestVectorAnalysis()
{
vtkm::testing::Testing::TryTypes(TestLinearFunctor(), vtkm::TypeListTagField());
vtkm::testing::Testing::TryTypes(TestCrossFunctor(), vtkm::TypeListTagFieldVec3());
vtkm::testing::Testing::TryTypes(TestVectorFunctor(), vtkm::TypeListTagFloatVec());
}
} // anonymous namespace