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 "vtkObject.h"
#include "vtkTypeTraits.h" // For type traits
#include "vtkSmartPointer.h" // For vtkSmartPointer.

#include "vtkMathConfigure.h" // For <cmath> and VTK_HAS_ISNAN etc.

#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 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);

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

  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 );

  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];
  }

  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] );};

  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 static_cast<float> (sqrt( v[0] * v[0] + v[1] * v[1] + v[2] * v[2] ) );};

  static double Norm(const double v[3]) {
    return sqrt( 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]);

  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 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 static_cast<float> (sqrt( x[0] * x[0] + x[1] * x[1] ) );};

  static double Norm2D(const double x[2]) {
    return 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 in[3],
                          float out[3]);
  static void Multiply3x3(const double A[3][3], const double in[3],
                          double out[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]);

  static void MultiplyMatrix(double **A, double **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(float A[3][3]);
  static double Determinant3x3(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]);


  static void Matrix3x3ToQuaternion(const float A[3][3], float quat[4]);
  static void Matrix3x3ToQuaternion(const double A[3][3], double quat[4]);


  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 int SolveLinearSystem(double **A, double *x, int size);

  static int InvertMatrix(double **A, double **AI, int size);

  static int InvertMatrix(double **A, double **AI, int size,
                          int *tmp1Size, double *tmp2Size);

  static int LUFactorLinearSystem(double **A, int *index, int size);

  static int LUFactorLinearSystem(double **A, int *index, int size,
                                  double *tmpSize);

  static void LUSolveLinearSystem(double **A, int *index,
                                  double *x, int size);

  static double EstimateMatrixCondition(double **A, int size);


  static int Jacobi(float **a, float *w, float **v);
  static int Jacobi(double **a, double *w, double **v);


  static int JacobiN(float **a, int n, float *w, float **v);
  static int JacobiN(double **a, int n, double *w, double **v);

  static int SolveHomogeneousLeastSquares(int numberOfSamples, double **xt, int xOrder,
                                double **mt);

  static int 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 double* RGBToHSV(const double rgb[3]) VTK_SIZEHINT(3);
  static double* RGBToHSV(double r, double g, double b) VTK_SIZEHINT(3);
  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 double* HSVToRGB(const double hsv[3]) VTK_SIZEHINT(3);
  static double* HSVToRGB(double h, double s, double v) VTK_SIZEHINT(3);
  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 double *LabToXYZ(const double lab[3]) VTK_SIZEHINT(3);


  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 double *XYZToLab(const double xyz[3]) VTK_SIZEHINT(3);


  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 double *XYZToRGB(const double xyz[3]) VTK_SIZEHINT(3);


  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 double *RGBToXYZ(const double rgb[3]) VTK_SIZEHINT(3);


  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 double *RGBToLab(const double rgb[3]) VTK_SIZEHINT(3);


  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 double *LabToRGB(const double lab[3]) VTK_SIZEHINT(3);


  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(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 int GetAdjustedScalarRange(
    vtkDataArray *array, int comp, double range[2]);

  static vtkTypeBool ExtentIsWithinOtherExtent(int extent1[6], int extent2[6]);

  static vtkTypeBool BoundsIsWithinOtherBounds(double bounds1[6], double bounds2[6], double delta[3]);

  static vtkTypeBool PointIsWithinBounds(double point[3], double bounds[6], double delta[3]);

  static int PlaneIntersectsAABB(double const bounds[6], double const normal[3],
    double const 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);
protected:
  vtkMath() {}
  ~vtkMath() override {}

  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 = ((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[3])
{
  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[3])
{
  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] ) );
}

//----------------------------------------------------------------------------
// 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(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(float A[3][3])
{
  return vtkDeterminant3x3( A );
}

//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(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 (value < min)
  {
    return min;
  }

  if (value > max)
  {
    return max;
  }

  return value;
}

//----------------------------------------------------------------------------
inline void vtkMath::ClampValue(double *value, const double range[2])
{
  if (value && range)
  {
    assert("pre: valid_range" && range[0]<=range[1]);

    if (*value < range[0])
    {
      *value = range[0];
    }
    else if (*value > range[1])
    {
      *value = 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]);

    if (value < range[0])
    {
      *clamped_value = range[0];
    }
    else if (value > range[1])
    {
      *clamped_value = range[1];
    }
    else
    {
      *clamped_value = value;
    }
  }
}

// ---------------------------------------------------------------------------
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
      if(value<range[0])
      {
          result=range[0];
      }
      else
      {
          if(value>range[1])
          {
              result=range[1];
          }
          else
          {
              result=value;
          }
      }

      // 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(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 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
  val = vtkMath::Max(val, static_cast<double>(vtkTypeTraits<OutT>::Min()));
  val = vtkMath::Min(val, static_cast<double>(vtkTypeTraits<OutT>::Max()));
  *ret = static_cast<OutT>((val >= 0.0) ? (val + 0.5) : (val - 0.5));
}
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
  val = vtkMath::Max(val, static_cast<double>(vtkTypeTraits<float>::Min()));
  val = vtkMath::Min(val, static_cast<double>(vtkTypeTraits<float>::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