mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
Merge branch 'vector-analysis' into 'master'
Vector Analysis Add basic vector analysis methods. See merge request !73
This commit is contained in:
commit
79960e8383
@ -32,6 +32,7 @@ set(headers
|
||||
TypeListTag.h
|
||||
Types.h
|
||||
TypeTraits.h
|
||||
VectorAnalysis.h
|
||||
VecTraits.h
|
||||
)
|
||||
|
||||
|
@ -154,7 +154,7 @@ vtkm::Pair<T1,T2> make_Pair(const T1 &firstSrc, const T2 &secondSrc)
|
||||
template<typename T, typename U>
|
||||
struct TypeTraits<vtkm::Pair<T,U> >
|
||||
{
|
||||
typedef TypeTraitsUnkownTag NumericTag;
|
||||
typedef TypeTraitsUnknownTag NumericTag;
|
||||
typedef TypeTraitsScalarTag DimensionalityTag;
|
||||
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
|
@ -26,7 +26,7 @@ namespace vtkm {
|
||||
|
||||
/// Tag used to identify types that aren't Real, Integer, Scalar or Vector.
|
||||
///
|
||||
struct TypeTraitsUnkownTag {};
|
||||
struct TypeTraitsUnknownTag {};
|
||||
|
||||
/// Tag used to identify types that store real (floating-point) numbers. A
|
||||
/// TypeTraits class will typedef this class to NumericTag if it stores real
|
||||
@ -64,13 +64,13 @@ public:
|
||||
/// \brief A tag to determing whether the type is integer or real.
|
||||
///
|
||||
/// This tag is either TypeTraitsRealTag or TypeTraitsIntegerTag.
|
||||
typedef TypeTraitsUnkownTag NumericTag;
|
||||
typedef TypeTraitsUnknownTag NumericTag;
|
||||
|
||||
/// \brief A tag to determine whether the type has multiple components.
|
||||
///
|
||||
/// This tag is either TypeTraitsScalarTag or TypeTraitsVectorTag. Scalars can
|
||||
/// also be treated as vectors.
|
||||
typedef TypeTraitsUnkownTag DimensionalityTag;
|
||||
typedef TypeTraitsUnknownTag DimensionalityTag;
|
||||
|
||||
VTKM_EXEC_CONT_EXPORT static T ZeroInitialization() { return T(); }
|
||||
};
|
||||
|
210
vtkm/VectorAnalysis.h
Normal file
210
vtkm/VectorAnalysis.h
Normal file
@ -0,0 +1,210 @@
|
||||
//=============================================================================
|
||||
//
|
||||
// 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 Sandia Corporation.
|
||||
// Copyright 2015 UT-Battelle, LLC.
|
||||
// Copyright 2015 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
|
||||
// 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_VectorAnalysis_h
|
||||
#define vtk_m_VectorAnalysis_h
|
||||
|
||||
// This header file defines math functions that deal with linear albegra funcitons
|
||||
|
||||
#include <vtkm/Math.h>
|
||||
#include <vtkm/Types.h>
|
||||
#include <vtkm/TypeTraits.h>
|
||||
#include <vtkm/VecTraits.h>
|
||||
|
||||
namespace vtkm {
|
||||
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
/// \brief Returns the linear interpolation of two values based on weight
|
||||
///
|
||||
/// \c Lerp interpolates return the linerar interpolation of v0 and v1 based on w. v0
|
||||
/// and v1 are scalars or vectors of same length. w can either be a scalar or a
|
||||
/// vector of the same length as x and y. If w is outside [0,1] then lerp
|
||||
/// extrapolates. If w=0 => v0 is returned if w=1 => v1 is returned.
|
||||
///
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
T Lerp(const T &v0, const T & v1, const T &w)
|
||||
{
|
||||
return v0 + w * (v1-v0);
|
||||
}
|
||||
template<typename T, vtkm::IdComponent N>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
vtkm::Vec<T,N> Lerp(const vtkm::Vec<T,N> &v0,
|
||||
const vtkm::Vec<T,N> &v1,
|
||||
const T &w)
|
||||
{
|
||||
return v0 + w * (v1-v0);
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
/// \brief Returns the square of the magnitude of a vector.
|
||||
///
|
||||
/// It is usually much faster to compute the square of the magnitude than the
|
||||
/// square, so you should use this function in place of Magnitude or RMagnitude
|
||||
/// when possible.
|
||||
///
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
typename vtkm::VecTraits<T>::ComponentType
|
||||
MagnitudeSquared(const T &x)
|
||||
{
|
||||
return vtkm::dot(x,x);
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
namespace detail {
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
T MagnitudeTemplate(T x, vtkm::TypeTraitsScalarTag)
|
||||
{
|
||||
return vtkm::Abs(x);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
typename vtkm::VecTraits<T>::ComponentType
|
||||
MagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
|
||||
{
|
||||
return vtkm::Sqrt(vtkm::MagnitudeSquared(x));
|
||||
}
|
||||
} // namespace detail
|
||||
|
||||
/// \brief Returns the magnitude of a vector.
|
||||
///
|
||||
/// It is usually much faster to compute MagnitudeSquared, so that should be
|
||||
/// substituted when possible (unless you are just going to take the square
|
||||
/// root, which would be besides the point). On some hardware it is also faster
|
||||
/// to find the reciprocal magnitude, so RMagnitude should be used if you
|
||||
/// actually plan to divide by the magnitude.
|
||||
///
|
||||
template<typename T>
|
||||
typename vtkm::VecTraits<T>::ComponentType
|
||||
Magnitude(const T &x)
|
||||
{
|
||||
return detail::MagnitudeTemplate(
|
||||
x, typename vtkm::TypeTraits<T>::DimensionalityTag());
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
namespace detail {
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
T RMagnitudeTemplate(T x, vtkm::TypeTraitsScalarTag)
|
||||
{
|
||||
return T(1)/vtkm::Abs(x);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
typename vtkm::VecTraits<T>::ComponentType
|
||||
RMagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
|
||||
{
|
||||
return vtkm::RSqrt(vtkm::MagnitudeSquared(x));
|
||||
}
|
||||
} // namespace detail
|
||||
|
||||
/// \brief Returns the reciprocal magnitude of a vector.
|
||||
///
|
||||
/// On some hardware RMagnitude is faster than Magnitude, but neither is
|
||||
/// as fast as MagnitudeSquared.
|
||||
///
|
||||
template<typename T>
|
||||
typename vtkm::VecTraits<T>::ComponentType
|
||||
RMagnitude(const T &x)
|
||||
{
|
||||
return detail::RMagnitudeTemplate(
|
||||
x, typename vtkm::TypeTraits<T>::DimensionalityTag());
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
namespace detail {
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
T NormalTemplate(T x, vtkm::TypeTraitsScalarTag)
|
||||
{
|
||||
return vtkm::CopySign(T(1), x);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
T NormalTemplate(const T &x, vtkm::TypeTraitsVectorTag)
|
||||
{
|
||||
return vtkm::RMagnitude(x)*x;
|
||||
}
|
||||
} // namespace detail
|
||||
|
||||
/// \brief Returns a normalized version of the given vector.
|
||||
///
|
||||
/// The resulting vector points in the same direction but has unit length.
|
||||
///
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
T Normal(const T &x)
|
||||
{
|
||||
return detail::NormalTemplate(
|
||||
x, typename vtkm::TypeTraits<T>::DimensionalityTag());
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
/// \brief Changes a vector to be normal.
|
||||
///
|
||||
/// The given vector is scaled to be unit length.
|
||||
///
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
void Normalize(T &x)
|
||||
{
|
||||
x = vtkm::Normal(x);
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------------------
|
||||
/// \brief Find the cross product of two vectors.
|
||||
///
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
vtkm::Vec<T,3> Cross(const vtkm::Vec<T,3> &x, const vtkm::Vec<T,3> &y)
|
||||
{
|
||||
return vtkm::Vec<T,3>(x[1]*y[2] - x[2]*y[1],
|
||||
x[2]*y[0] - x[0]*y[2],
|
||||
x[0]*y[1] - x[1]*y[0]);
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
/// \brief Find the normal of a triangle.
|
||||
///
|
||||
/// Given three coordinates in space, which, unless degenerate, uniquely define
|
||||
/// a triangle and the plane the triangle is on, returns a vector perpendicular
|
||||
/// to that triangle/plane.
|
||||
///
|
||||
template<typename T>
|
||||
VTKM_EXEC_CONT_EXPORT
|
||||
vtkm::Vec<T,3> TriangleNormal(const vtkm::Vec<T,3> &a,
|
||||
const vtkm::Vec<T,3> &b,
|
||||
const vtkm::Vec<T,3> &c)
|
||||
{
|
||||
return vtkm::Cross(b-a, c-a);
|
||||
}
|
||||
|
||||
|
||||
} // namespace vtkm
|
||||
|
||||
#endif //vtk_m_VectorAnalysis_h
|
@ -37,6 +37,7 @@ set(unit_tests
|
||||
UnitTestTypeListTag.cxx
|
||||
UnitTestTypes.cxx
|
||||
UnitTestTypeTraits.cxx
|
||||
UnitTestVectorAnalysis.cxx
|
||||
UnitTestVecTraits.cxx
|
||||
)
|
||||
|
||||
|
224
vtkm/testing/UnitTestVectorAnalysis.cxx
Normal file
224
vtkm/testing/UnitTestVectorAnalysis.cxx
Normal file
@ -0,0 +1,224 @@
|
||||
//=============================================================================
|
||||
//
|
||||
// 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 Sandia Corporation.
|
||||
// Copyright 2015 UT-Battelle, LLC.
|
||||
// Copyright 2015 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
|
||||
// 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.
|
||||
//
|
||||
//=============================================================================
|
||||
|
||||
#include <vtkm/VectorAnalysis.h>
|
||||
|
||||
#include <vtkm/Types.h>
|
||||
#include <vtkm/VecTraits.h>
|
||||
|
||||
#include <vtkm/testing/Testing.h>
|
||||
|
||||
#include <math.h>
|
||||
|
||||
namespace {
|
||||
|
||||
namespace internal{
|
||||
|
||||
template<typename VectorType>
|
||||
typename vtkm::VecTraits<VectorType>::ComponentType
|
||||
MyMag(const VectorType &vt)
|
||||
{
|
||||
typedef vtkm::VecTraits<VectorType> Traits;
|
||||
double total = 0.0;
|
||||
for (vtkm::IdComponent index = 0; index < Traits::NUM_COMPONENTS; ++index)
|
||||
{
|
||||
total += Traits::GetComponent(vt,index) * Traits::GetComponent(vt,index);
|
||||
}
|
||||
return static_cast<typename Traits::ComponentType>(sqrt(total));
|
||||
}
|
||||
|
||||
template<typename VectorType>
|
||||
VectorType MyNormal(const VectorType& vt)
|
||||
{
|
||||
typedef vtkm::VecTraits<VectorType> Traits;
|
||||
double total = 0.0;
|
||||
for (vtkm::IdComponent index = 0; index < Traits::NUM_COMPONENTS; ++index)
|
||||
{
|
||||
total += Traits::GetComponent(vt,index) * Traits::GetComponent(vt,index);
|
||||
}
|
||||
typename Traits::ComponentType mag = internal::MyMag(vt);
|
||||
VectorType temp = vt;
|
||||
for (vtkm::IdComponent index = 0; index < Traits::NUM_COMPONENTS; ++index)
|
||||
{
|
||||
Traits::SetComponent(temp, index, Traits::GetComponent(vt,index)/mag);
|
||||
}
|
||||
return temp;
|
||||
}
|
||||
|
||||
template<typename T,typename W>
|
||||
T MyLerp (const T &a, const T &b, const W &w)
|
||||
{
|
||||
return (W(1)-w)*a + w*b;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
template<typename VectorType>
|
||||
void TestVector(const VectorType& vector)
|
||||
{
|
||||
typedef typename vtkm::VecTraits<VectorType>::ComponentType ComponentType;
|
||||
|
||||
std::cout << "Testing " << vector << std::endl;
|
||||
|
||||
//to do have to implement a norm and normalized call to verify the math ones
|
||||
//against
|
||||
ComponentType magnitude = vtkm::Magnitude(vector);
|
||||
ComponentType magnitudeCompare = internal::MyMag(vector);
|
||||
VTKM_TEST_ASSERT(test_equal(magnitude, magnitudeCompare),
|
||||
"Magnitude failed test.");
|
||||
|
||||
ComponentType magnitudeSquared = vtkm::MagnitudeSquared(vector);
|
||||
VTKM_TEST_ASSERT(test_equal(magnitude*magnitude, magnitudeSquared),
|
||||
"Magnitude squared test failed.");
|
||||
|
||||
if (magnitudeSquared > 0)
|
||||
{
|
||||
ComponentType rmagnitude = vtkm::RMagnitude(vector);
|
||||
VTKM_TEST_ASSERT(test_equal(1/magnitude, rmagnitude),
|
||||
"Reciprical magnitude failed.");
|
||||
|
||||
VTKM_TEST_ASSERT(test_equal(vtkm::Normal(vector),
|
||||
internal::MyNormal(vector)),
|
||||
"Normalized vector failed test.");
|
||||
|
||||
VectorType normalizedVector = vector;
|
||||
vtkm::Normalize(normalizedVector);
|
||||
VTKM_TEST_ASSERT(test_equal(normalizedVector, internal::MyNormal(vector)),
|
||||
"Inplace Normalized vector failed test.");
|
||||
}
|
||||
}
|
||||
|
||||
template <typename VectorType>
|
||||
void TestLerp(const VectorType& a,
|
||||
const VectorType& b,
|
||||
const VectorType& w,
|
||||
const typename vtkm::VecTraits<VectorType>::ComponentType& wS)
|
||||
{
|
||||
VectorType lhs = internal::MyLerp(a,b,w);
|
||||
VectorType rhs = vtkm::Lerp(a,b,w);
|
||||
VTKM_TEST_ASSERT(test_equal(lhs, rhs),
|
||||
"Vectors with Vector weight do not lerp() correctly");
|
||||
|
||||
VectorType lhsS = internal::MyLerp(a,b,wS);
|
||||
VectorType rhsS = vtkm::Lerp(a,b,wS);
|
||||
VTKM_TEST_ASSERT(test_equal(lhsS, rhsS),
|
||||
"Vectors with Scalar weight do not lerp() correctly");
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
void TestCross(const vtkm::Vec<T,3> &x, const vtkm::Vec<T,3> &y)
|
||||
{
|
||||
typedef vtkm::Vec<T,3> Vec3;
|
||||
Vec3 cross = vtkm::Cross(x, y);
|
||||
|
||||
std::cout << "Testing " << x << " x " << y << " = " << cross << std::endl;
|
||||
|
||||
// The cross product result should be perpendicular to input vectors.
|
||||
VTKM_TEST_ASSERT(test_equal(vtkm::dot(cross,x), T(0.0)),
|
||||
"Cross product not perpendicular.");
|
||||
VTKM_TEST_ASSERT(test_equal(vtkm::dot(cross,y), T(0.0)),
|
||||
"Cross product not perpendicular.");
|
||||
|
||||
// The length of cross product should be the lengths of the input vectors
|
||||
// times the sin of the angle between them.
|
||||
T sinAngle =
|
||||
vtkm::Magnitude(cross) * vtkm::RMagnitude(x) * vtkm::RMagnitude(y);
|
||||
|
||||
// The dot product is likewise the lengths of the input vectors times the
|
||||
// cos of the angle between them.
|
||||
T cosAngle = vtkm::dot(x,y) * vtkm::RMagnitude(x) * vtkm::RMagnitude(y);
|
||||
|
||||
// Test that these are the actual sin and cos of the same angle with a
|
||||
// basic trigonometric identity.
|
||||
VTKM_TEST_ASSERT(test_equal(sinAngle*sinAngle + cosAngle*cosAngle, T(1.0)),
|
||||
"Bad cross product length.");
|
||||
|
||||
// Test finding the normal to a triangle (similar to cross product).
|
||||
Vec3 normal = vtkm::TriangleNormal(x, y, Vec3(0, 0, 0));
|
||||
VTKM_TEST_ASSERT(test_equal(vtkm::dot(normal, x-y), T(0.0)),
|
||||
"Triangle normal is not really normal.");
|
||||
}
|
||||
|
||||
struct TestLinearFunctor
|
||||
{
|
||||
template <typename T>
|
||||
void operator()(const T&) const
|
||||
{
|
||||
typedef vtkm::VecTraits<T> Traits;
|
||||
const vtkm::IdComponent NUM_COMPONENTS = Traits::NUM_COMPONENTS;
|
||||
typedef typename Traits::ComponentType ComponentType;
|
||||
|
||||
T zeroVector = T(ComponentType(0));
|
||||
T normalizedVector = T(vtkm::RSqrt(ComponentType(NUM_COMPONENTS)));
|
||||
T posVec = TestValue(1, T());
|
||||
T negVec = -TestValue(2, T());
|
||||
|
||||
TestVector(zeroVector);
|
||||
TestVector(normalizedVector);
|
||||
TestVector(posVec);
|
||||
TestVector(negVec);
|
||||
|
||||
T weight(ComponentType(0.5));
|
||||
ComponentType weightS(0.5);
|
||||
TestLerp(zeroVector,normalizedVector,weight,weightS);
|
||||
TestLerp(zeroVector,posVec,weight,weightS);
|
||||
TestLerp(zeroVector,negVec,weight,weightS);
|
||||
|
||||
TestLerp(normalizedVector,zeroVector,weight,weightS);
|
||||
TestLerp(normalizedVector,posVec,weight,weightS);
|
||||
TestLerp(normalizedVector,negVec,weight,weightS);
|
||||
|
||||
TestLerp(posVec,zeroVector,weight,weightS);
|
||||
TestLerp(posVec,normalizedVector,weight,weightS);
|
||||
TestLerp(posVec,negVec,weight,weightS);
|
||||
|
||||
TestLerp(negVec,zeroVector,weight,weightS);
|
||||
TestLerp(negVec,normalizedVector,weight,weightS);
|
||||
TestLerp(negVec,posVec,weight,weightS);
|
||||
}
|
||||
};
|
||||
|
||||
struct TestCrossFunctor
|
||||
{
|
||||
template<typename VectorType>
|
||||
void operator()(const VectorType &) const
|
||||
{
|
||||
TestCross(VectorType(1.0f,0.0f,0.0f), VectorType(0.0f,1.0f,0.0f));
|
||||
TestCross(VectorType(1.0f,2.0f,3.0f), VectorType(-3.0f,-1.0f,1.0f));
|
||||
TestCross(VectorType(0.0f,0.0f,1.0f), VectorType(0.001f,0.01f,2.0f));
|
||||
}
|
||||
};
|
||||
|
||||
void TestVectorAnalysis()
|
||||
{
|
||||
vtkm::testing::Testing::TryTypes(TestLinearFunctor(),
|
||||
vtkm::TypeListTagField());
|
||||
vtkm::testing::Testing::TryTypes(TestCrossFunctor(),
|
||||
vtkm::TypeListTagFieldVec3());
|
||||
}
|
||||
|
||||
} // anonymous namespace
|
||||
|
||||
int UnitTestVectorAnalysis(int, char *[])
|
||||
{
|
||||
return vtkm::testing::Testing::Run(TestVectorAnalysis);
|
||||
}
|
Loading…
Reference in New Issue
Block a user