VectorAnalysis.h 10.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
//=============================================================================
//
//  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.
//
11
//  Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
12 13 14
//  Copyright 2015 UT-Battelle, LLC.
//  Copyright 2015 Los Alamos National Security.
//
15
//  Under the terms of Contract DE-NA0003525 with NTESS,
16 17 18 19 20 21 22 23 24
//  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

luz.paz's avatar
luz.paz committed
25
// This header file defines math functions that deal with linear albegra functions
26 27 28

#include <vtkm/Math.h>
#include <vtkm/TypeTraits.h>
29
#include <vtkm/Types.h>
30 31
#include <vtkm/VecTraits.h>

32 33
namespace vtkm
{
34 35 36 37 38 39 40 41 42

// ----------------------------------------------------------------------------
/// \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.
///
43
template <typename ValueType, typename WeightType>
44 45 46
inline VTKM_EXEC_CONT ValueType Lerp(const ValueType& value0,
                                     const ValueType& value1,
                                     const WeightType& weight)
47
{
48 49 50
  using ScalarType = typename detail::FloatingPointReturnType<ValueType>::Type;
  return static_cast<ValueType>((WeightType(1) - weight) * static_cast<ScalarType>(value0) +
                                weight * static_cast<ScalarType>(value1));
51
}
52 53 54 55
template <typename ValueType, vtkm::IdComponent N, typename WeightType>
VTKM_EXEC_CONT vtkm::Vec<ValueType, N> Lerp(const vtkm::Vec<ValueType, N>& value0,
                                            const vtkm::Vec<ValueType, N>& value1,
                                            const WeightType& weight)
56
{
57
  return (WeightType(1) - weight) * value0 + weight * value1;
58
}
59 60 61 62
template <typename ValueType, vtkm::IdComponent N>
VTKM_EXEC_CONT vtkm::Vec<ValueType, N> Lerp(const vtkm::Vec<ValueType, N>& value0,
                                            const vtkm::Vec<ValueType, N>& value1,
                                            const vtkm::Vec<ValueType, N>& weight)
63
{
64
  const vtkm::Vec<ValueType, N> One(ValueType(1));
65
  return (One - weight) * value0 + weight * value1;
66 67 68 69 70 71 72 73 74
}

// ----------------------------------------------------------------------------
/// \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.
///
75
template <typename T>
76
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type MagnitudeSquared(const T& x)
77
{
78
  using U = typename detail::FloatingPointReturnType<T>::Type;
79
  return static_cast<U>(vtkm::Dot(x, x));
80 81 82
}

// ----------------------------------------------------------------------------
83 84 85 86
namespace detail
{
template <typename T>
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type MagnitudeTemplate(
87 88
  T x,
  vtkm::TypeTraitsScalarTag)
89
{
90
  return static_cast<typename detail::FloatingPointReturnType<T>::Type>(vtkm::Abs(x));
91 92
}

93 94
template <typename T>
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type MagnitudeTemplate(
95 96
  const T& x,
  vtkm::TypeTraitsVectorTag)
97 98 99
{
  return vtkm::Sqrt(vtkm::MagnitudeSquared(x));
}
100

101 102 103 104 105 106 107 108 109 110
} // 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.
///
111 112
template <typename T>
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type Magnitude(const T& x)
113
{
114
  return detail::MagnitudeTemplate(x, typename vtkm::TypeTraits<T>::DimensionalityTag());
115 116 117
}

// ----------------------------------------------------------------------------
118
namespace detail
119
{
120 121
template <typename T>
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type RMagnitudeTemplate(
122 123
  T x,
  vtkm::TypeTraitsScalarTag)
124 125
{
  return T(1) / vtkm::Abs(x);
126 127
}

128 129
template <typename T>
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type RMagnitudeTemplate(
130 131
  const T& x,
  vtkm::TypeTraitsVectorTag)
132 133 134 135 136 137 138 139 140 141
{
  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.
///
142 143
template <typename T>
VTKM_EXEC_CONT typename detail::FloatingPointReturnType<T>::Type RMagnitude(const T& x)
144
{
145
  return detail::RMagnitudeTemplate(x, typename vtkm::TypeTraits<T>::DimensionalityTag());
146 147 148
}

// ----------------------------------------------------------------------------
149 150 151 152
namespace detail
{
template <typename T>
VTKM_EXEC_CONT T NormalTemplate(T x, vtkm::TypeTraitsScalarTag)
153 154 155 156
{
  return vtkm::CopySign(T(1), x);
}

157 158
template <typename T>
VTKM_EXEC_CONT T NormalTemplate(const T& x, vtkm::TypeTraitsVectorTag)
159
{
160
  return vtkm::RMagnitude(x) * x;
161 162 163 164 165 166 167
}
} // namespace detail

/// \brief Returns a normalized version of the given vector.
///
/// The resulting vector points in the same direction but has unit length.
///
168 169
template <typename T>
VTKM_EXEC_CONT T Normal(const T& x)
170
{
171
  return detail::NormalTemplate(x, typename vtkm::TypeTraits<T>::DimensionalityTag());
172 173 174 175 176 177 178
}

// ----------------------------------------------------------------------------
/// \brief Changes a vector to be normal.
///
/// The given vector is scaled to be unit length.
///
179 180
template <typename T>
VTKM_EXEC_CONT void Normalize(T& x)
181 182 183 184 185 186 187
{
  x = vtkm::Normal(x);
}

// ----------------------------------------------------------------------------
/// \brief Find the cross product of two vectors.
///
188 189
template <typename T>
VTKM_EXEC_CONT vtkm::Vec<typename detail::FloatingPointReturnType<T>::Type, 3> Cross(
190 191
  const vtkm::Vec<T, 3>& x,
  const vtkm::Vec<T, 3>& y)
192
{
193 194
  return vtkm::Vec<typename detail::FloatingPointReturnType<T>::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]);
195 196 197 198 199 200 201 202 203
}

//-----------------------------------------------------------------------------
/// \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.
///
204
template <typename T>
205 206
VTKM_EXEC_CONT vtkm::Vec<typename detail::FloatingPointReturnType<T>::Type, 3>
TriangleNormal(const vtkm::Vec<T, 3>& a, const vtkm::Vec<T, 3>& b, const vtkm::Vec<T, 3>& c)
207
{
208
  return vtkm::Cross(b - a, c - a);
209 210
}

211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
//-----------------------------------------------------------------------------
/// \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;
}

289 290 291
} // namespace vtkm

#endif //vtk_m_VectorAnalysis_h