vtkMath.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-FileCopyrightText: Copyright 2011 Sandia Corporation
// SPDX-License-Identifier: LicenseRef-BSD-3-Clause-Sandia-USGov
#ifndef vtkMath_h
#define vtkMath_h
 
#include "vtkCommonCoreModule.h" // For export macro
#include "vtkMathPrivate.hxx"    // For Matrix meta-class helpers
#include "vtkMatrixUtilities.h"  // For Matrix wrapping / mapping
#include "vtkObject.h"
#include "vtkSmartPointer.h" // For vtkSmartPointer.
#include "vtkTypeTraits.h"   // For type traits
 
#include "vtkMathConfigure.h" // For <cmath> and VTK_HAS_ISNAN etc.
 
#include <algorithm>   // for std::clamp
#include <cassert>     // assert() in inline implementations.
#include <type_traits> // for type_traits
 
#ifndef DBL_MIN
#define VTK_DBL_MIN 2.2250738585072014e-308
#else // DBL_MIN
#define VTK_DBL_MIN DBL_MIN
#endif // DBL_MIN
 
#ifndef DBL_EPSILON
#define VTK_DBL_EPSILON 2.2204460492503131e-16
#else // DBL_EPSILON
#define VTK_DBL_EPSILON DBL_EPSILON
#endif // DBL_EPSILON
 
#ifndef VTK_DBL_EPSILON
#ifndef DBL_EPSILON
#define VTK_DBL_EPSILON 2.2204460492503131e-16
#else // DBL_EPSILON
#define VTK_DBL_EPSILON DBL_EPSILON
#endif // DBL_EPSILON
#endif // VTK_DBL_EPSILON
 
VTK_ABI_NAMESPACE_BEGIN
class vtkDataArray;
class vtkPoints;
class vtkMathInternal;
class vtkMinimalStandardRandomSequence;
class vtkBoxMuellerRandomSequence;
VTK_ABI_NAMESPACE_END
 
namespace vtk_detail
{
VTK_ABI_NAMESPACE_BEGIN
// forward declaration
template <typename OutT>
void RoundDoubleToIntegralIfNecessary(double val, OutT* ret);
VTK_ABI_NAMESPACE_END
} // end namespace vtk_detail
 
VTK_ABI_NAMESPACE_BEGIN
class VTKCOMMONCORE_EXPORT vtkMath : public vtkObject
{
public:
  static vtkMath* New();
  vtkTypeMacro(vtkMath, vtkObject);
  void PrintSelf(ostream& os, vtkIndent indent) override;
 
  static constexpr double Pi() { return 3.141592653589793; }
 
 
  static float RadiansFromDegrees(float degrees);
  static double RadiansFromDegrees(double degrees);
 
 
  static float DegreesFromRadians(float radians);
  static double DegreesFromRadians(double radians);
 
#if 1
  static int Round(float f) { return static_cast<int>(f + (f >= 0.0 ? 0.5 : -0.5)); }
  static int Round(double f) { return static_cast<int>(f + (f >= 0.0 ? 0.5 : -0.5)); }
#endif
 
  template <typename OutT>
  static void RoundDoubleToIntegralIfNecessary(double val, OutT* ret)
  {
    // Can't specialize template methods in a template class, so we move the
    // implementations to a external namespace.
    vtk_detail::RoundDoubleToIntegralIfNecessary(val, ret);
  }
 
  static int Floor(double x);
 
  static int Ceil(double x);
 
  static int CeilLog2(vtkTypeUInt64 x);
 
  template <class T>
  static T Min(const T& a, const T& b);
 
  template <class T>
  static T Max(const T& a, const T& b);
 
  static bool IsPowerOfTwo(vtkTypeUInt64 x);
 
  static int NearestPowerOfTwo(int x);
 
  static vtkTypeInt64 Factorial(int N);
 
  static vtkTypeInt64 Binomial(int m, int n);
 
  static int* BeginCombination(int m, int n);
 
  static int NextCombination(int m, int n, int* combination);
 
  static void FreeCombination(int* combination);
 
  static void RandomSeed(int s);
 
  static int GetSeed();
 
  static double Random();
 
  static double Random(double min, double max);
 
  static double Gaussian();
 
  static double Gaussian(double mean, double std);
 
  template <class VectorT1, class VectorT2>
  static void Assign(const VectorT1& a, VectorT2&& b)
  {
    b[0] = a[0];
    b[1] = a[1];
    b[2] = a[2];
  }
 
  static void Assign(const double a[3], double b[3]) { vtkMath::Assign<>(a, b); }
 
  static void Add(const float a[3], const float b[3], float c[3])
  {
    for (int i = 0; i < 3; ++i)
    {
      c[i] = a[i] + b[i];
    }
  }
 
  static void Add(const double a[3], const double b[3], double c[3])
  {
    for (int i = 0; i < 3; ++i)
    {
      c[i] = a[i] + b[i];
    }
  }
 
  template <class VectorT1, class VectorT2, class VectorT3>
  static void Add(VectorT1&& a, VectorT2&& b, VectorT3& c)
  {
    for (int i = 0; i < 3; ++i)
    {
      c[i] = a[i] + b[i];
    }
  }
 
  static void Subtract(const float a[3], const float b[3], float c[3])
  {
    for (int i = 0; i < 3; ++i)
    {
      c[i] = a[i] - b[i];
    }
  }
 
  static void Subtract(const double a[3], const double b[3], double c[3])
  {
    for (int i = 0; i < 3; ++i)
    {
      c[i] = a[i] - b[i];
    }
  }
 
  template <class VectorT1, class VectorT2, class VectorT3>
  static void Subtract(const VectorT1& a, const VectorT2& b, VectorT3&& c)
  {
    c[0] = a[0] - b[0];
    c[1] = a[1] - b[1];
    c[2] = a[2] - b[2];
  }
 
  static void MultiplyScalar(float a[3], float s)
  {
    for (int i = 0; i < 3; ++i)
    {
      a[i] *= s;
    }
  }
 
  static void MultiplyScalar2D(float a[2], float s)
  {
    for (int i = 0; i < 2; ++i)
    {
      a[i] *= s;
    }
  }
 
  static void MultiplyScalar(double a[3], double s)
  {
    for (int i = 0; i < 3; ++i)
    {
      a[i] *= s;
    }
  }
 
  static void MultiplyScalar2D(double a[2], double s)
  {
    for (int i = 0; i < 2; ++i)
    {
      a[i] *= s;
    }
  }
 
  static float Dot(const float a[3], const float b[3])
  {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  }
 
  static double Dot(const double a[3], const double b[3])
  {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  }
 
  template <typename ReturnTypeT = double, typename TupleRangeT1, typename TupleRangeT2,
    typename EnableT = typename std::conditional<!std::is_pointer<TupleRangeT1>::value &&
        !std::is_array<TupleRangeT1>::value,
      TupleRangeT1, TupleRangeT2>::type::value_type>
  static ReturnTypeT Dot(const TupleRangeT1& a, const TupleRangeT2& b)
  {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  }
 
  static void Outer(const float a[3], const float b[3], float c[3][3])
  {
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        c[i][j] = a[i] * b[j];
      }
    }
  }
 
  static void Outer(const double a[3], const double b[3], double c[3][3])
  {
    for (int i = 0; i < 3; ++i)
    {
      for (int j = 0; j < 3; ++j)
      {
        c[i][j] = a[i] * b[j];
      }
    }
  }
 
  template <class VectorT1, class VectorT2, class VectorT3>
  static void Cross(VectorT1&& a, VectorT2&& b, VectorT3& c);
 
  static void Cross(const float a[3], const float b[3], float c[3]);
 
  static void Cross(const double a[3], const double b[3], double c[3]);
 
 
  static float Norm(const float* x, int n);
  static double Norm(const double* x, int n);
 
  static float Norm(const float v[3]) { return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); }
 
  static double Norm(const double v[3])
  {
    return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
  }
 
  template <typename ReturnTypeT = double, typename TupleRangeT>
  static ReturnTypeT SquaredNorm(const TupleRangeT& v)
  {
    return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
  }
 
  static float Normalize(float v[3]);
 
  static double Normalize(double v[3]);
 
 
  static void Perpendiculars(const double v1[3], double v2[3], double v3[3], double theta);
  static void Perpendiculars(const float v1[3], float v2[3], float v3[3], double theta);
 
 
  static bool ProjectVector(const float a[3], const float b[3], float projection[3]);
  static bool ProjectVector(const double a[3], const double b[3], double projection[3]);
 
 
  static bool ProjectVector2D(const float a[2], const float b[2], float projection[2]);
  static bool ProjectVector2D(const double a[2], const double b[2], double projection[2]);
 
  template <typename ReturnTypeT = double, typename TupleRangeT1, typename TupleRangeT2,
    typename EnableT = typename std::conditional<!std::is_pointer<TupleRangeT1>::value &&
        !std::is_array<TupleRangeT1>::value,
      TupleRangeT1, TupleRangeT2>::type::value_type>
  static ReturnTypeT Distance2BetweenPoints(const TupleRangeT1& p1, const TupleRangeT2& p2);
 
  static float Distance2BetweenPoints(const float p1[3], const float p2[3]);
 
  static double Distance2BetweenPoints(const double p1[3], const double p2[3]);
 
  static double AngleBetweenVectors(const double v1[3], const double v2[3]);
 
  static double SignedAngleBetweenVectors(
    const double v1[3], const double v2[3], const double vn[3]);
 
  static double GaussianAmplitude(double variance, double distanceFromMean);
 
  static double GaussianAmplitude(double mean, double variance, double position);
 
  static double GaussianWeight(double variance, double distanceFromMean);
 
  static double GaussianWeight(double mean, double variance, double position);
 
  static float Dot2D(const float x[2], const float y[2]) { return x[0] * y[0] + x[1] * y[1]; }
 
  static double Dot2D(const double x[2], const double y[2]) { return x[0] * y[0] + x[1] * y[1]; }
 
  static void Outer2D(const float x[2], const float y[2], float A[2][2])
  {
    for (int i = 0; i < 2; ++i)
    {
      for (int j = 0; j < 2; ++j)
      {
        A[i][j] = x[i] * y[j];
      }
    }
  }
 
  static void Outer2D(const double x[2], const double y[2], double A[2][2])
  {
    for (int i = 0; i < 2; ++i)
    {
      for (int j = 0; j < 2; ++j)
      {
        A[i][j] = x[i] * y[j];
      }
    }
  }
 
  static float Norm2D(const float x[2]) { return std::sqrt(x[0] * x[0] + x[1] * x[1]); }
 
  static double Norm2D(const double x[2]) { return std::sqrt(x[0] * x[0] + x[1] * x[1]); }
 
  static float Normalize2D(float v[2]);
 
  static double Normalize2D(double v[2]);
 
  static float Determinant2x2(const float c1[2], const float c2[2])
  {
    return c1[0] * c2[1] - c2[0] * c1[1];
  }
 
 
  static double Determinant2x2(double a, double b, double c, double d) { return a * d - b * c; }
  static double Determinant2x2(const double c1[2], const double c2[2])
  {
    return c1[0] * c2[1] - c2[0] * c1[1];
  }
 
 
  static void LUFactor3x3(float A[3][3], int index[3]);
  static void LUFactor3x3(double A[3][3], int index[3]);
 
 
  static void LUSolve3x3(const float A[3][3], const int index[3], float x[3]);
  static void LUSolve3x3(const double A[3][3], const int index[3], double x[3]);
 
 
  static void LinearSolve3x3(const float A[3][3], const float x[3], float y[3]);
  static void LinearSolve3x3(const double A[3][3], const double x[3], double y[3]);
 
 
  static void Multiply3x3(const float A[3][3], const float v[3], float u[3]);
  static void Multiply3x3(const double A[3][3], const double v[3], double u[3]);
 
 
  static void Multiply3x3(const float A[3][3], const float B[3][3], float C[3][3]);
  static void Multiply3x3(const double A[3][3], const double B[3][3], double C[3][3]);
 
  template <int RowsT, int MidDimT, int ColsT,
    class LayoutT1 = vtkMatrixUtilities::Layout::Identity,
    class LayoutT2 = vtkMatrixUtilities::Layout::Identity, class MatrixT1, class MatrixT2,
    class MatrixT3>
  static void MultiplyMatrix(const MatrixT1& M1, const MatrixT2& M2, MatrixT3&& M3)
  {
    vtkMathPrivate::MultiplyMatrix<RowsT, MidDimT, ColsT, LayoutT1, LayoutT2>::Compute(M1, M2, M3);
  }
 
  template <int RowsT, int ColsT, class LayoutT = vtkMatrixUtilities::Layout::Identity,
    class MatrixT, class VectorT1, class VectorT2>
  static void MultiplyMatrixWithVector(const MatrixT& M, const VectorT1& X, VectorT2&& Y)
  {
    vtkMathPrivate::MultiplyMatrix<RowsT, ColsT, 1, LayoutT>::Compute(M, X, Y);
  }
 
  template <class ScalarT, int SizeT, class VectorT1, class VectorT2>
  static ScalarT Dot(const VectorT1& x, const VectorT2& y)
  {
    return vtkMathPrivate::ContractRowWithCol<ScalarT, 1, SizeT, 1, 0, 0,
      vtkMatrixUtilities::Layout::Identity, vtkMatrixUtilities::Layout::Transpose>::Compute(x, y);
  }
 
  template <int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity, class MatrixT>
  static typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Determinant(
    const MatrixT& M)
  {
    return vtkMathPrivate::Determinant<SizeT, LayoutT>::Compute(M);
  }
 
  template <int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity, class MatrixT1,
    class MatrixT2>
  static void InvertMatrix(const MatrixT1& M1, MatrixT2&& M2)
  {
    vtkMathPrivate::InvertMatrix<SizeT, LayoutT>::Compute(M1, M2);
  }
 
  template <int RowsT, int ColsT, class LayoutT = vtkMatrixUtilities::Layout::Identity,
    class MatrixT, class VectorT1, class VectorT2>
  static void LinearSolve(const MatrixT& M, const VectorT1& x, VectorT2& y)
  {
    vtkMathPrivate::LinearSolve<RowsT, ColsT, LayoutT>::Compute(M, x, y);
  }
 
  template <class ScalarT, int SizeT, class LayoutT = vtkMatrixUtilities::Layout::Identity,
    class VectorT1, class MatrixT, class VectorT2>
  static ScalarT Dot(const VectorT1& x, const MatrixT& M, const VectorT2& y)
  {
    ScalarT tmp[SizeT];
    vtkMathPrivate::MultiplyMatrix<SizeT, SizeT, 1, LayoutT>::Compute(M, y, tmp);
    return vtkMathPrivate::ContractRowWithCol<ScalarT, 1, SizeT, 1, 0, 0,
      vtkMatrixUtilities::Layout::Identity, vtkMatrixUtilities::Layout::Transpose>::Compute(x, tmp);
  }
 
  static void MultiplyMatrix(const double* const* A, const double* const* B, unsigned int rowA,
    unsigned int colA, unsigned int rowB, unsigned int colB, double** C);
 
 
  static void Transpose3x3(const float A[3][3], float AT[3][3]);
  static void Transpose3x3(const double A[3][3], double AT[3][3]);
 
 
  static void Invert3x3(const float A[3][3], float AI[3][3]);
  static void Invert3x3(const double A[3][3], double AI[3][3]);
 
 
  static void Identity3x3(float A[3][3]);
  static void Identity3x3(double A[3][3]);
 
 
  static double Determinant3x3(const float A[3][3]);
  static double Determinant3x3(const double A[3][3]);
 
  static float Determinant3x3(const float c1[3], const float c2[3], const float c3[3]);
 
  static double Determinant3x3(const double c1[3], const double c2[3], const double c3[3]);
 
  static double Determinant3x3(double a1, double a2, double a3, double b1, double b2, double b3,
    double c1, double c2, double c3);
 
 
  static void QuaternionToMatrix3x3(const float quat[4], float A[3][3]);
  static void QuaternionToMatrix3x3(const double quat[4], double A[3][3]);
  template <class QuaternionT, class MatrixT,
    class EnableT = typename std::enable_if<!vtkMatrixUtilities::MatrixIs2DArray<MatrixT>()>::type>
  static void QuaternionToMatrix3x3(const QuaternionT& q, MatrixT&& A);
 
 
  static void Matrix3x3ToQuaternion(const float A[3][3], float quat[4]);
  static void Matrix3x3ToQuaternion(const double A[3][3], double quat[4]);
  template <class MatrixT, class QuaternionT,
    class EnableT = typename std::enable_if<!vtkMatrixUtilities::MatrixIs2DArray<MatrixT>()>::type>
  static void Matrix3x3ToQuaternion(const MatrixT& A, QuaternionT&& q);
 
 
  static void MultiplyQuaternion(const float q1[4], const float q2[4], float q[4]);
  static void MultiplyQuaternion(const double q1[4], const double q2[4], double q[4]);
 
 
  static void RotateVectorByNormalizedQuaternion(const float v[3], const float q[4], float r[3]);
  static void RotateVectorByNormalizedQuaternion(const double v[3], const double q[4], double r[3]);
 
 
  static void RotateVectorByWXYZ(const float v[3], const float q[4], float r[3]);
  static void RotateVectorByWXYZ(const double v[3], const double q[4], double r[3]);
 
 
  static void Orthogonalize3x3(const float A[3][3], float B[3][3]);
  static void Orthogonalize3x3(const double A[3][3], double B[3][3]);
 
 
  static void Diagonalize3x3(const float A[3][3], float w[3], float V[3][3]);
  static void Diagonalize3x3(const double A[3][3], double w[3], double V[3][3]);
 
 
  static void SingularValueDecomposition3x3(
    const float A[3][3], float U[3][3], float w[3], float VT[3][3]);
  static void SingularValueDecomposition3x3(
    const double A[3][3], double U[3][3], double w[3], double VT[3][3]);
 
  static vtkTypeBool SolveLinearSystemGEPP2x2(
    double a00, double a01, double a10, double a11, double b0, double b1, double& x0, double& x1);
 
  static vtkTypeBool SolveLinearSystem(double** A, double* x, int size);
 
  static vtkTypeBool InvertMatrix(double** A, double** AI, int size);
 
  static vtkTypeBool InvertMatrix(
    double** A, double** AI, int size, int* tmp1Size, double* tmp2Size);
 
  static vtkTypeBool LUFactorLinearSystem(double** A, int* index, int size);
 
  static vtkTypeBool LUFactorLinearSystem(double** A, int* index, int size, double* tmpSize);
 
  static void LUSolveLinearSystem(double** A, int* index, double* x, int size);
 
  static double EstimateMatrixCondition(const double* const* A, int size);
 
 
  static vtkTypeBool Jacobi(float** a, float* w, float** v);
  static vtkTypeBool Jacobi(double** a, double* w, double** v);
 
 
  static vtkTypeBool JacobiN(float** a, int n, float* w, float** v);
  static vtkTypeBool JacobiN(double** a, int n, double* w, double** v);
 
  static vtkTypeBool SolveHomogeneousLeastSquares(
    int numberOfSamples, double** xt, int xOrder, double** mt);
 
  static vtkTypeBool SolveLeastSquares(int numberOfSamples, double** xt, int xOrder, double** yt,
    int yOrder, double** mt, int checkHomogeneous = 1);
 
 
  static void RGBToHSV(const float rgb[3], float hsv[3])
  {
    RGBToHSV(rgb[0], rgb[1], rgb[2], hsv, hsv + 1, hsv + 2);
  }
  static void RGBToHSV(float r, float g, float b, float* h, float* s, float* v);
  static void RGBToHSV(const double rgb[3], double hsv[3])
  {
    RGBToHSV(rgb[0], rgb[1], rgb[2], hsv, hsv + 1, hsv + 2);
  }
  static void RGBToHSV(double r, double g, double b, double* h, double* s, double* v);
 
 
  static void HSVToRGB(const float hsv[3], float rgb[3])
  {
    HSVToRGB(hsv[0], hsv[1], hsv[2], rgb, rgb + 1, rgb + 2);
  }
  static void HSVToRGB(float h, float s, float v, float* r, float* g, float* b);
  static void HSVToRGB(const double hsv[3], double rgb[3])
  {
    HSVToRGB(hsv[0], hsv[1], hsv[2], rgb, rgb + 1, rgb + 2);
  }
  static void HSVToRGB(double h, double s, double v, double* r, double* g, double* b);
 
 
  static void LabToXYZ(const double lab[3], double xyz[3])
  {
    LabToXYZ(lab[0], lab[1], lab[2], xyz + 0, xyz + 1, xyz + 2);
  }
  static void LabToXYZ(double L, double a, double b, double* x, double* y, double* z);
 
 
  static void XYZToLab(const double xyz[3], double lab[3])
  {
    XYZToLab(xyz[0], xyz[1], xyz[2], lab + 0, lab + 1, lab + 2);
  }
  static void XYZToLab(double x, double y, double z, double* L, double* a, double* b);
 
 
  static void XYZToRGB(const double xyz[3], double rgb[3])
  {
    XYZToRGB(xyz[0], xyz[1], xyz[2], rgb + 0, rgb + 1, rgb + 2);
  }
  static void XYZToRGB(double x, double y, double z, double* r, double* g, double* b);
 
 
  static void RGBToXYZ(const double rgb[3], double xyz[3])
  {
    RGBToXYZ(rgb[0], rgb[1], rgb[2], xyz + 0, xyz + 1, xyz + 2);
  }
  static void RGBToXYZ(double r, double g, double b, double* x, double* y, double* z);
 
 
  static void RGBToLab(const double rgb[3], double lab[3])
  {
    RGBToLab(rgb[0], rgb[1], rgb[2], lab + 0, lab + 1, lab + 2);
  }
  static void RGBToLab(double red, double green, double blue, double* L, double* a, double* b);
 
 
  static void LabToRGB(const double lab[3], double rgb[3])
  {
    LabToRGB(lab[0], lab[1], lab[2], rgb + 0, rgb + 1, rgb + 2);
  }
  static void LabToRGB(double L, double a, double b, double* red, double* green, double* blue);
 
 
  static void UninitializeBounds(double bounds[6])
  {
    bounds[0] = 1.0;
    bounds[1] = -1.0;
    bounds[2] = 1.0;
    bounds[3] = -1.0;
    bounds[4] = 1.0;
    bounds[5] = -1.0;
  }
 
 
  static vtkTypeBool AreBoundsInitialized(const double bounds[6])
  {
    if (bounds[1] - bounds[0] < 0.0)
    {
      return 0;
    }
    return 1;
  }
 
  template <class T>
  static T ClampValue(const T& value, const T& min, const T& max);
 
 
  static void ClampValue(double* value, const double range[2]);
  static void ClampValue(double value, const double range[2], double* clamped_value);
  static void ClampValues(double* values, int nb_values, const double range[2]);
  static void ClampValues(
    const double* values, int nb_values, const double range[2], double* clamped_values);
 
  static double ClampAndNormalizeValue(double value, const double range[2]);
 
  template <class T1, class T2>
  static void TensorFromSymmetricTensor(const T1 symmTensor[6], T2 tensor[9]);
 
  template <class T>
  static void TensorFromSymmetricTensor(T tensor[9]);
 
  static int GetScalarTypeFittingRange(
    double range_min, double range_max, double scale = 1.0, double shift = 0.0);
 
  static vtkTypeBool GetAdjustedScalarRange(vtkDataArray* array, int comp, double range[2]);
 
  static vtkTypeBool ExtentIsWithinOtherExtent(const int extent1[6], const int extent2[6]);
 
  static vtkTypeBool BoundsIsWithinOtherBounds(
    const double bounds1[6], const double bounds2[6], const double delta[3]);
 
  static vtkTypeBool PointIsWithinBounds(
    const double point[3], const double bounds[6], const double delta[3]);
 
  static int PlaneIntersectsAABB(
    const double bounds[6], const double normal[3], const double point[3]);
 
  static double Solve3PointCircle(
    const double p1[3], const double p2[3], const double p3[3], double center[3]);
 
  static double Inf();
 
  static double NegInf();
 
  static double Nan();
 
  static vtkTypeBool IsInf(double x);
 
  static vtkTypeBool IsNan(double x);
 
  static bool IsFinite(double x);
 
  static int QuadraticRoot(double a, double b, double c, double min, double max, double* u);
 
  static vtkIdType ComputeGCD(vtkIdType m, vtkIdType n) { return (n ? ComputeGCD(n, m % n) : m); }
 
  enum class ConvolutionMode
  {
    FULL,
    SAME,
    VALID
  };
 
  template <class Iter1, class Iter2, class Iter3>
  static void Convolve1D(Iter1 beginSample, Iter1 endSample, Iter2 beginKernel, Iter2 endKernel,
    Iter3 beginOut, Iter3 endOut, ConvolutionMode mode = ConvolutionMode::FULL)
  {
    int sampleSize = std::distance(beginSample, endSample);
    int kernelSize = std::distance(beginKernel, endKernel);
    int outSize = std::distance(beginOut, endOut);
 
    if (sampleSize <= 0 || kernelSize <= 0 || outSize <= 0)
    {
      return;
    }
 
    int begin = 0;
    int end = outSize;
 
    switch (mode)
    {
      case ConvolutionMode::SAME:
        begin = static_cast<int>(std::ceil((std::min)(sampleSize, kernelSize) / 2.0)) - 1;
        end = begin + (std::max)(sampleSize, kernelSize);
        break;
      case ConvolutionMode::VALID:
        begin = (std::min)(sampleSize, kernelSize) - 1;
        end = begin + std::abs(sampleSize - kernelSize) + 1;
        break;
      case ConvolutionMode::FULL:
      default:
        break;
    }
 
    for (int i = begin; i < end; i++)
    {
      Iter3 out = beginOut + i - begin;
      *out = 0;
      for (int j = (std::max)(i - sampleSize + 1, 0); j <= (std::min)(i, kernelSize - 1); j++)
      {
        *out += *(beginSample + (i - j)) * *(beginKernel + j);
      }
    }
  }
 
  static void GetPointAlongLine(double result[3], double p1[3], double p2[3], const double offset)
  {
    double directionVector[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
    vtkMath::Normalize(directionVector);
    result[0] = p2[0] + (offset * directionVector[0]);
    result[1] = p2[1] + (offset * directionVector[1]);
    result[2] = p2[2] + (offset * directionVector[2]);
  }
 
protected:
  vtkMath() = default;
  ~vtkMath() override = default;
 
  static vtkSmartPointer<vtkMathInternal> Internal;
 
private:
  vtkMath(const vtkMath&) = delete;
  void operator=(const vtkMath&) = delete;
};
 
//----------------------------------------------------------------------------
inline float vtkMath::RadiansFromDegrees(float x)
{
  return x * 0.017453292f;
}
 
//----------------------------------------------------------------------------
inline double vtkMath::RadiansFromDegrees(double x)
{
  return x * 0.017453292519943295;
}
 
//----------------------------------------------------------------------------
inline float vtkMath::DegreesFromRadians(float x)
{
  return x * 57.2957795131f;
}
 
//----------------------------------------------------------------------------
inline double vtkMath::DegreesFromRadians(double x)
{
  return x * 57.29577951308232;
}
 
//----------------------------------------------------------------------------
inline bool vtkMath::IsPowerOfTwo(vtkTypeUInt64 x)
{
  return ((x != 0) & ((x & (x - 1)) == 0));
}
 
//----------------------------------------------------------------------------
// Credit goes to Peter Hart and William Lewis on comp.lang.python 1997
inline int vtkMath::NearestPowerOfTwo(int x)
{
  unsigned int z = static_cast<unsigned int>(((x > 0) ? x - 1 : 0));
  z |= z >> 1;
  z |= z >> 2;
  z |= z >> 4;
  z |= z >> 8;
  z |= z >> 16;
  return static_cast<int>(z + 1);
}
 
//----------------------------------------------------------------------------
// Modify the trunc() operation provided by static_cast<int>() to get floor(),
// Note that in C++ conditions evaluate to values of 1 or 0 (true or false).
inline int vtkMath::Floor(double x)
{
  int i = static_cast<int>(x);
  return i - (i > x);
}
 
//----------------------------------------------------------------------------
// Modify the trunc() operation provided by static_cast<int>() to get ceil(),
// Note that in C++ conditions evaluate to values of 1 or 0 (true or false).
inline int vtkMath::Ceil(double x)
{
  int i = static_cast<int>(x);
  return i + (i < x);
}
 
//----------------------------------------------------------------------------
template <class T>
inline T vtkMath::Min(const T& a, const T& b)
{
  return (b <= a ? b : a);
}
 
//----------------------------------------------------------------------------
template <class T>
inline T vtkMath::Max(const T& a, const T& b)
{
  return (b > a ? b : a);
}
 
//----------------------------------------------------------------------------
inline float vtkMath::Normalize(float v[3])
{
  float den = vtkMath::Norm(v);
  if (den != 0.0)
  {
    for (int i = 0; i < 3; ++i)
    {
      v[i] /= den;
    }
  }
  return den;
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Normalize(double v[3])
{
  double den = vtkMath::Norm(v);
  if (den != 0.0)
  {
    for (int i = 0; i < 3; ++i)
    {
      v[i] /= den;
    }
  }
  return den;
}
 
//----------------------------------------------------------------------------
inline float vtkMath::Normalize2D(float v[2])
{
  float den = vtkMath::Norm2D(v);
  if (den != 0.0)
  {
    for (int i = 0; i < 2; ++i)
    {
      v[i] /= den;
    }
  }
  return den;
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Normalize2D(double v[2])
{
  double den = vtkMath::Norm2D(v);
  if (den != 0.0)
  {
    for (int i = 0; i < 2; ++i)
    {
      v[i] /= den;
    }
  }
  return den;
}
 
//----------------------------------------------------------------------------
inline float vtkMath::Determinant3x3(const float c1[3], const float c2[3], const float c3[3])
{
  return c1[0] * c2[1] * c3[2] + c2[0] * c3[1] * c1[2] + c3[0] * c1[1] * c2[2] -
    c1[0] * c3[1] * c2[2] - c2[0] * c1[1] * c3[2] - c3[0] * c2[1] * c1[2];
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(const double c1[3], const double c2[3], const double c3[3])
{
  return c1[0] * c2[1] * c3[2] + c2[0] * c3[1] * c1[2] + c3[0] * c1[1] * c2[2] -
    c1[0] * c3[1] * c2[2] - c2[0] * c1[1] * c3[2] - c3[0] * c2[1] * c1[2];
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(
  double a1, double a2, double a3, double b1, double b2, double b3, double c1, double c2, double c3)
{
  return (a1 * vtkMath::Determinant2x2(b2, b3, c2, c3) -
    b1 * vtkMath::Determinant2x2(a2, a3, c2, c3) + c1 * vtkMath::Determinant2x2(a2, a3, b2, b3));
}
 
//----------------------------------------------------------------------------
inline float vtkMath::Distance2BetweenPoints(const float p1[3], const float p2[3])
{
  return ((p1[0] - p2[0]) * (p1[0] - p2[0]) + (p1[1] - p2[1]) * (p1[1] - p2[1]) +
    (p1[2] - p2[2]) * (p1[2] - p2[2]));
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Distance2BetweenPoints(const double p1[3], const double p2[3])
{
  return ((p1[0] - p2[0]) * (p1[0] - p2[0]) + (p1[1] - p2[1]) * (p1[1] - p2[1]) +
    (p1[2] - p2[2]) * (p1[2] - p2[2]));
}
 
//----------------------------------------------------------------------------
template <typename ReturnTypeT, typename TupleRangeT1, typename TupleRangeT2, typename EnableT>
inline ReturnTypeT vtkMath::Distance2BetweenPoints(const TupleRangeT1& p1, const TupleRangeT2& p2)
{
  return ((p1[0] - p2[0]) * (p1[0] - p2[0]) + (p1[1] - p2[1]) * (p1[1] - p2[1]) +
    (p1[2] - p2[2]) * (p1[2] - p2[2]));
}
 
//----------------------------------------------------------------------------
template <class VectorT1, class VectorT2, class VectorT3>
void vtkMath::Cross(VectorT1&& a, VectorT2&& b, VectorT3& c)
{
  c[0] = a[1] * b[2] - a[2] * b[1];
  c[1] = a[2] * b[0] - a[0] * b[2];
  c[2] = a[0] * b[1] - a[1] * b[0];
}
 
//----------------------------------------------------------------------------
// Cross product of two 3-vectors. Result (a x b) is stored in c[3].
inline void vtkMath::Cross(const float a[3], const float b[3], float c[3])
{
  float Cx = a[1] * b[2] - a[2] * b[1];
  float Cy = a[2] * b[0] - a[0] * b[2];
  float Cz = a[0] * b[1] - a[1] * b[0];
  c[0] = Cx;
  c[1] = Cy;
  c[2] = Cz;
}
 
//----------------------------------------------------------------------------
// Cross product of two 3-vectors. Result (a x b) is stored in c[3].
inline void vtkMath::Cross(const double a[3], const double b[3], double c[3])
{
  double Cx = a[1] * b[2] - a[2] * b[1];
  double Cy = a[2] * b[0] - a[0] * b[2];
  double Cz = a[0] * b[1] - a[1] * b[0];
  c[0] = Cx;
  c[1] = Cy;
  c[2] = Cz;
}
 
//----------------------------------------------------------------------------
template <class T>
inline double vtkDeterminant3x3(const T A[3][3])
{
  return A[0][0] * A[1][1] * A[2][2] + A[1][0] * A[2][1] * A[0][2] + A[2][0] * A[0][1] * A[1][2] -
    A[0][0] * A[2][1] * A[1][2] - A[1][0] * A[0][1] * A[2][2] - A[2][0] * A[1][1] * A[0][2];
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(const float A[3][3])
{
  return vtkDeterminant3x3(A);
}
 
//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(const double A[3][3])
{
  return vtkDeterminant3x3(A);
}
 
//----------------------------------------------------------------------------
template <class T>
inline T vtkMath::ClampValue(const T& value, const T& min, const T& max)
{
  assert("pre: valid_range" && min <= max);
 
#if __cplusplus >= 201703L
  return std::clamp(value, min, max);
#else
  // compilers are good at optimizing the ternary operator,
  // use '<' since it is preferred by STL for custom types
  T v = (min < value ? value : min);
  return (v < max ? v : max);
#endif
}
 
//----------------------------------------------------------------------------
inline void vtkMath::ClampValue(double* value, const double range[2])
{
  if (value && range)
  {
    assert("pre: valid_range" && range[0] <= range[1]);
 
    *value = vtkMath::ClampValue(*value, range[0], range[1]);
  }
}
 
//----------------------------------------------------------------------------
inline void vtkMath::ClampValue(double value, const double range[2], double* clamped_value)
{
  if (range && clamped_value)
  {
    assert("pre: valid_range" && range[0] <= range[1]);
 
    *clamped_value = vtkMath::ClampValue(value, range[0], range[1]);
  }
}
 
// ---------------------------------------------------------------------------
inline double vtkMath::ClampAndNormalizeValue(double value, const double range[2])
{
  assert("pre: valid_range" && range[0] <= range[1]);
 
  double result;
  if (range[0] == range[1])
  {
    result = 0.0;
  }
  else
  {
    // clamp
    result = vtkMath::ClampValue(value, range[0], range[1]);
 
    // normalize
    result = (result - range[0]) / (range[1] - range[0]);
  }
 
  assert("post: valid_result" && result >= 0.0 && result <= 1.0);
 
  return result;
}
 
//-----------------------------------------------------------------------------
template <class T1, class T2>
inline void vtkMath::TensorFromSymmetricTensor(const T1 symmTensor[9], T2 tensor[9])
{
  for (int i = 0; i < 3; ++i)
  {
    tensor[4 * i] = symmTensor[i];
  }
  tensor[1] = tensor[3] = symmTensor[3];
  tensor[2] = tensor[6] = symmTensor[5];
  tensor[5] = tensor[7] = symmTensor[4];
}
 
//-----------------------------------------------------------------------------
template <class T>
inline void vtkMath::TensorFromSymmetricTensor(T tensor[9])
{
  tensor[6] = tensor[5]; // XZ
  tensor[7] = tensor[4]; // YZ
  tensor[8] = tensor[2]; // ZZ
  tensor[4] = tensor[1]; // YY
  tensor[5] = tensor[7]; // YZ
  tensor[2] = tensor[6]; // XZ
  tensor[1] = tensor[3]; // XY
}
VTK_ABI_NAMESPACE_END
 
namespace
{
template <class QuaternionT, class MatrixT>
inline void vtkQuaternionToMatrix3x3(const QuaternionT& quat, MatrixT& A)
{
  typedef typename vtkMatrixUtilities::ScalarTypeExtractor<MatrixT>::value_type Scalar;
 
  Scalar ww = quat[0] * quat[0];
  Scalar wx = quat[0] * quat[1];
  Scalar wy = quat[0] * quat[2];
  Scalar wz = quat[0] * quat[3];
 
  Scalar xx = quat[1] * quat[1];
  Scalar yy = quat[2] * quat[2];
  Scalar zz = quat[3] * quat[3];
 
  Scalar xy = quat[1] * quat[2];
  Scalar xz = quat[1] * quat[3];
  Scalar yz = quat[2] * quat[3];
 
  Scalar rr = xx + yy + zz;
  // normalization factor, just in case quaternion was not normalized
  Scalar f = 1 / (ww + rr);
  Scalar s = (ww - rr) * f;
  f *= 2;
 
  typedef vtkMatrixUtilities::Wrapper<3, 3, MatrixT> Wrapper;
 
  Wrapper::template Get<0, 0>(A) = xx * f + s;
  Wrapper::template Get<1, 0>(A) = (xy + wz) * f;
  Wrapper::template Get<2, 0>(A) = (xz - wy) * f;
 
  Wrapper::template Get<0, 1>(A) = (xy - wz) * f;
  Wrapper::template Get<1, 1>(A) = yy * f + s;
  Wrapper::template Get<2, 1>(A) = (yz + wx) * f;
 
  Wrapper::template Get<0, 2>(A) = (xz + wy) * f;
  Wrapper::template Get<1, 2>(A) = (yz - wx) * f;
  Wrapper::template Get<2, 2>(A) = zz * f + s;
}
} // anonymous namespace
 
VTK_ABI_NAMESPACE_BEGIN
//------------------------------------------------------------------------------
inline void vtkMath::QuaternionToMatrix3x3(const float quat[4], float A[3][3])
{
  vtkQuaternionToMatrix3x3(quat, A);
}
 
//------------------------------------------------------------------------------
inline void vtkMath::QuaternionToMatrix3x3(const double quat[4], double A[3][3])
{
  vtkQuaternionToMatrix3x3(quat, A);
}
 
//-----------------------------------------------------------------------------
template <class QuaternionT, class MatrixT, class EnableT>
inline void vtkMath::QuaternionToMatrix3x3(const QuaternionT& q, MatrixT&& A)
{
  vtkQuaternionToMatrix3x3(q, A);
}
VTK_ABI_NAMESPACE_END
 
namespace
{
//------------------------------------------------------------------------------
//  The solution is based on
//  Berthold K. P. Horn (1987),
//  "Closed-form solution of absolute orientation using unit quaternions,"
//  Journal of the Optical Society of America A, 4:629-642
template <class MatrixT, class QuaternionT>
inline void vtkMatrix3x3ToQuaternion(const MatrixT& A, QuaternionT& quat)
{
  typedef typename vtkMatrixUtilities::ScalarTypeExtractor<QuaternionT>::value_type Scalar;
 
  Scalar N[4][4];
 
  typedef vtkMatrixUtilities::Wrapper<3, 3, MatrixT> Wrapper;
 
  // on-diagonal elements
  N[0][0] = Wrapper::template Get<0, 0>(A) + Wrapper::template Get<1, 1>(A) +
    Wrapper::template Get<2, 2>(A);
  N[1][1] = Wrapper::template Get<0, 0>(A) - Wrapper::template Get<1, 1>(A) -
    Wrapper::template Get<2, 2>(A);
  N[2][2] = -Wrapper::template Get<0, 0>(A) + Wrapper::template Get<1, 1>(A) -
    Wrapper::template Get<2, 2>(A);
  N[3][3] = -Wrapper::template Get<0, 0>(A) - Wrapper::template Get<1, 1>(A) +
    Wrapper::template Get<2, 2>(A);
 
  // off-diagonal elements
  N[0][1] = N[1][0] = Wrapper::template Get<2, 1>(A) - Wrapper::template Get<1, 2>(A);
  N[0][2] = N[2][0] = Wrapper::template Get<0, 2>(A) - Wrapper::template Get<2, 0>(A);
  N[0][3] = N[3][0] = Wrapper::template Get<1, 0>(A) - Wrapper::template Get<0, 1>(A);
 
  N[1][2] = N[2][1] = Wrapper::template Get<1, 0>(A) + Wrapper::template Get<0, 1>(A);
  N[1][3] = N[3][1] = Wrapper::template Get<0, 2>(A) + Wrapper::template Get<2, 0>(A);
  N[2][3] = N[3][2] = Wrapper::template Get<2, 1>(A) + Wrapper::template Get<1, 2>(A);
 
  Scalar eigenvectors[4][4], eigenvalues[4];
 
  // convert into format that JacobiN can use,
  // then use Jacobi to find eigenvalues and eigenvectors
  Scalar *NTemp[4], *eigenvectorsTemp[4];
  for (int i = 0; i < 4; ++i)
  {
    NTemp[i] = N[i];
    eigenvectorsTemp[i] = eigenvectors[i];
  }
  vtkMath::JacobiN(NTemp, 4, eigenvalues, eigenvectorsTemp);
 
  // the first eigenvector is the one we want
  quat[0] = eigenvectors[0][0];
  quat[1] = eigenvectors[1][0];
  quat[2] = eigenvectors[2][0];
  quat[3] = eigenvectors[3][0];
}
} // anonymous namespace
 
VTK_ABI_NAMESPACE_BEGIN
//------------------------------------------------------------------------------
inline void vtkMath::Matrix3x3ToQuaternion(const float A[3][3], float quat[4])
{
  vtkMatrix3x3ToQuaternion(A, quat);
}
 
//------------------------------------------------------------------------------
inline void vtkMath::Matrix3x3ToQuaternion(const double A[3][3], double quat[4])
{
  vtkMatrix3x3ToQuaternion(A, quat);
}
 
//-----------------------------------------------------------------------------
template <class MatrixT, class QuaternionT, class EnableT>
inline void vtkMath::Matrix3x3ToQuaternion(const MatrixT& A, QuaternionT&& q)
{
  vtkMatrix3x3ToQuaternion(A, q);
}
VTK_ABI_NAMESPACE_END
 
namespace vtk_detail
{
VTK_ABI_NAMESPACE_BEGIN
// Can't specialize templates inside a template class, so we move the impl here.
template <typename OutT>
void RoundDoubleToIntegralIfNecessary(double val, OutT* ret)
{ // OutT is integral -- clamp and round
  if (!vtkMath::IsNan(val))
  {
    double min = static_cast<double>(vtkTypeTraits<OutT>::Min());
    double max = static_cast<double>(vtkTypeTraits<OutT>::Max());
    val = vtkMath::ClampValue(val, min, max);
    *ret = static_cast<OutT>((val >= 0.0) ? (val + 0.5) : (val - 0.5));
  }
  else
    *ret = 0;
}
template <>
inline void RoundDoubleToIntegralIfNecessary(double val, double* retVal)
{ // OutT is double: passthrough
  *retVal = val;
}
template <>
inline void RoundDoubleToIntegralIfNecessary(double val, float* retVal)
{ // OutT is float -- just clamp (as doubles, then the cast to float is well-defined.)
  if (!vtkMath::IsNan(val))
  {
    double min = static_cast<double>(vtkTypeTraits<float>::Min());
    double max = static_cast<double>(vtkTypeTraits<float>::Max());
    val = vtkMath::ClampValue(val, min, max);
  }
 
  *retVal = static_cast<float>(val);
}
VTK_ABI_NAMESPACE_END
} // end namespace vtk_detail
 
VTK_ABI_NAMESPACE_BEGIN
//-----------------------------------------------------------------------------
#if defined(VTK_HAS_ISINF) || defined(VTK_HAS_STD_ISINF)
#define VTK_MATH_ISINF_IS_INLINE
inline vtkTypeBool vtkMath::IsInf(double x)
{
#if defined(VTK_HAS_STD_ISINF)
  return std::isinf(x);
#else
  return (isinf(x) != 0);    // Force conversion to bool
#endif
}
#endif
 
//-----------------------------------------------------------------------------
#if defined(VTK_HAS_ISNAN) || defined(VTK_HAS_STD_ISNAN)
#define VTK_MATH_ISNAN_IS_INLINE
inline vtkTypeBool vtkMath::IsNan(double x)
{
#if defined(VTK_HAS_STD_ISNAN)
  return std::isnan(x);
#else
  return (isnan(x) != 0);    // Force conversion to bool
#endif
}
#endif
 
//-----------------------------------------------------------------------------
#if defined(VTK_HAS_ISFINITE) || defined(VTK_HAS_STD_ISFINITE) || defined(VTK_HAS_FINITE)
#define VTK_MATH_ISFINITE_IS_INLINE
inline bool vtkMath::IsFinite(double x)
{
#if defined(VTK_HAS_STD_ISFINITE)
  return std::isfinite(x);
#elif defined(VTK_HAS_ISFINITE)
  return (isfinite(x) != 0); // Force conversion to bool
#else
  return (finite(x) != 0); // Force conversion to bool
#endif
}
#endif
 
VTK_ABI_NAMESPACE_END
#endif