vtkMath.h

Go to the documentation of this file.

/*=========================================================================
 
  Program:   Visualization Toolkit
  Module:    vtkMath.h
 
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm 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 2011 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
  license for use of this work by or on behalf of the
  U.S. Government. Redistribution and use in source and binary forms, with
  or without modification, are permitted provided that this Notice and any
  statement of authorship are reproduced on all copies.
 
  Contact: pppebay@sandia.gov,dcthomp@sandia.gov
 
=========================================================================*/
#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.
 
#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
 
class vtkDataArray;
class vtkPoints;
class vtkMathInternal;
class vtkMinimalStandardRandomSequence;
class vtkBoxMuellerRandomSequence;
 
namespace vtk_detail
{
// forward declaration
template <typename OutT>
void RoundDoubleToIntegralIfNecessary(double val, OutT* ret);
} // end namespace vtk_detail
 
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];
    }
  }
 
  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];
      }
    }
  }
 
  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(const double variance, const double distanceFromMean);
 
  static double GaussianAmplitude(const double mean, const double variance, const double position);
 
  static double GaussianWeight(const double variance, const double distanceFromMean);
 
  static double GaussianWeight(const double mean, const double variance, const 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);
 
  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);
      }
    }
  }
 
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]));
}
 
//----------------------------------------------------------------------------
// 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
}
 
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
 
//------------------------------------------------------------------------------
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);
}
 
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
 
//------------------------------------------------------------------------------
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);
}
 
namespace vtk_detail
{
// 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);
}
} // end namespace vtk_detail
 
//-----------------------------------------------------------------------------
#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
 
#endif