Add basic vector analysis methods.

vtkm/CMakeLists.txt

@@ -32,6 +32,7 @@ set(headers
   TypeListTag.h
   Types.h
   TypeTraits.h
+  VectorAnalysis.h
   VecTraits.h
   )
 

vtkm/VectorAnalysis.h

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

vtkm/testing/CMakeLists.txt

@@ -37,6 +37,7 @@ set(unit_tests
   UnitTestTypeListTag.cxx
   UnitTestTypes.cxx
   UnitTestTypeTraits.cxx
+  UnitTestVectorAnalysis.cxx
   UnitTestVecTraits.cxx
   )
 

vtkm/testing/UnitTestVectorAnalysis.cxx

@@ -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);
+}