//============================================================================= // // 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_VectorAnalysis_h #define vtk_m_VectorAnalysis_h // This header file defines math functions that deal with linear albegra functions #include #include #include #include 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 inline VTKM_EXEC_CONT ValueType Lerp(const ValueType& value0, const ValueType& value1, const WeightType& weight) { using ScalarType = typename detail::FloatingPointReturnType::Type; return static_cast((WeightType(1) - weight) * static_cast(value0) + weight * static_cast(value1)); } template VTKM_EXEC_CONT vtkm::Vec Lerp(const vtkm::Vec& value0, const vtkm::Vec& value1, const WeightType& weight) { return (WeightType(1) - weight) * value0 + weight * value1; } template VTKM_EXEC_CONT vtkm::Vec Lerp(const vtkm::Vec& value0, const vtkm::Vec& value1, const vtkm::Vec& weight) { const vtkm::Vec One(ValueType(1)); return (One - weight) * value0 + weight * value1; } // ---------------------------------------------------------------------------- /// \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 VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type MagnitudeSquared(const T& x) { using U = typename detail::FloatingPointReturnType::Type; return static_cast(vtkm::Dot(x, x)); } // ---------------------------------------------------------------------------- namespace detail { template VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type MagnitudeTemplate( T x, vtkm::TypeTraitsScalarTag) { return static_cast::Type>(vtkm::Abs(x)); } template VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type 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 VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type Magnitude(const T& x) { return detail::MagnitudeTemplate(x, typename vtkm::TypeTraits::DimensionalityTag()); } // ---------------------------------------------------------------------------- namespace detail { template VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type RMagnitudeTemplate( T x, vtkm::TypeTraitsScalarTag) { return T(1) / vtkm::Abs(x); } template VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type 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 VTKM_EXEC_CONT typename detail::FloatingPointReturnType::Type RMagnitude(const T& x) { return detail::RMagnitudeTemplate(x, typename vtkm::TypeTraits::DimensionalityTag()); } // ---------------------------------------------------------------------------- namespace detail { template VTKM_EXEC_CONT T NormalTemplate(T x, vtkm::TypeTraitsScalarTag) { return vtkm::CopySign(T(1), x); } template VTKM_EXEC_CONT 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 VTKM_EXEC_CONT T Normal(const T& x) { return detail::NormalTemplate(x, typename vtkm::TypeTraits::DimensionalityTag()); } // ---------------------------------------------------------------------------- /// \brief Changes a vector to be normal. /// /// The given vector is scaled to be unit length. /// template VTKM_EXEC_CONT void Normalize(T& x) { x = vtkm::Normal(x); } // ---------------------------------------------------------------------------- /// \brief Find the cross product of two vectors. /// template VTKM_EXEC_CONT vtkm::Vec::Type, 3> Cross( const vtkm::Vec& x, const vtkm::Vec& y) { return vtkm::Vec::Type, 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 VTKM_EXEC_CONT vtkm::Vec::Type, 3> TriangleNormal(const vtkm::Vec& a, const vtkm::Vec& b, const vtkm::Vec& c) { 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 VTKM_EXEC_CONT vtkm::Vec Project(const vtkm::Vec& v, const vtkm::Vec& u) { T uu = vtkm::Dot(u, u); T uv = vtkm::Dot(u, v); T factor = uv / uu; vtkm::Vec 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 VTKM_EXEC_CONT T ProjectedDistance(const vtkm::Vec& v, const vtkm::Vec& 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 VTKM_EXEC_CONT int Orthonormalize(const vtkm::Vec, N>& inputs, vtkm::Vec, N>& outputs, T tol = static_cast(1e-6)) { int j = 0; // j is the number of non-zero-length, non-collinear inputs encountered. vtkm::Vec, 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{ 0. }; } return j; } } // namespace vtkm #endif //vtk_m_VectorAnalysis_h