dox/Common/vtkMath.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkMath.h,v $

  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 2007 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@ca.sandia.gov,dcthomp@sandia.gov

=========================================================================*/
#ifndef __vtkMath_h
#define __vtkMath_h

#include "vtkObject.h"

#ifndef DBL_EPSILON
#  define VTK_DBL_EPSILON    2.2204460492503131e-16
#else  // DBL_EPSILON
#  define VTK_DBL_EPSILON    DBL_EPSILON
#endif  // DBL_EPSILON

class vtkDataArray;

class VTK_COMMON_EXPORT vtkMath : public vtkObject
{
public:
  static vtkMath *New();
  vtkTypeRevisionMacro(vtkMath,vtkObject);
  void PrintSelf(ostream& os, vtkIndent indent);


  static float Pi() {return 3.14159265358979f;};
  static float DegreesToRadians() {return 0.017453292f;};
  static float RadiansToDegrees() {return 57.2957795131f;};


  static double DoubleDegreesToRadians() {return 0.017453292519943295;};
  static double DoublePi() {return 3.1415926535897932384626;};
  static double DoubleRadiansToDegrees() {return 57.29577951308232;};


  static int Round(float f) {
    return static_cast<int>(f + (f >= 0 ? 0.5 : -0.5)); }
  static int Round(double f) {
    return static_cast<int>(f + (f >= 0 ? 0.5 : -0.5)); }

  static int Floor(double 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 float Dot(const float x[3], const float y[3]) {
    return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]);};


  static double Dot(const double x[3], const double y[3]) {
    return (x[0]*y[0] + x[1]*y[1] + x[2]*y[2]);};
  

  static void Outer(const float x[3], const float y[3], float A[3][3]) {
    for (int i=0; i < 3; i++)
      for (int j=0; j < 3; j++)
        A[i][j] = x[i]*y[j];
  }
  // Description:
  // Outer product of two 3-vectors (double-precision version).
  static void Outer(const double x[3], const double y[3], double A[3][3]) {
    for (int i=0; i < 3; i++)
      for (int j=0; j < 3; j++)
        A[i][j] = x[i]*y[j];
  }

  static void Cross(const float x[3], const float y[3], float z[3]);

  static void Cross(const double x[3], const double y[3], double z[3]);


  static float Norm(const float* x, int n); 
  static double Norm(const double* x, int n); 


  static float Norm(const float x[3]) {
    return static_cast<float> (sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]));};
  

  static double Norm(const double x[3]) {
    return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);};
  
  static float Normalize(float x[3]);

  static double Normalize(double x[3]);


  static void Perpendiculars(const double x[3], double y[3], double z[3], 
                             double theta);
  static void Perpendiculars(const float x[3], float y[3], float z[3],
                             double theta);

  static float Distance2BetweenPoints(const float x[3], const float y[3]);

  static double Distance2BetweenPoints(const double x[3], const double y[3]);


  static float Dot2D(const float x[3], const float y[3]) {
    return (x[0]*y[0] + x[1]*y[1]);};
  

  static double Dot2D(const double x[3], const double y[3]) {
    return (x[0]*y[0] + x[1]*y[1]);};


  static void Outer2D(const float x[3], const float y[3], float A[3][3]) {
    for (int i=0; i < 2; i++)
      for (int j=0; j < 2; j++)
        A[i][j] = x[i]*y[j];
  }
  // Description:
  // Outer product of two 2-vectors (float version). z-comp is ignored
  static void Outer2D(const double x[3], const double y[3], double A[3][3]) {
    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[3]) {
    return static_cast<float> (sqrt(x[0]*x[0] + x[1]*x[1]));};


  static double Norm2D(const double x[3]) {
    return sqrt(x[0]*x[0] + x[1]*x[1]);};

  static float Normalize2D(float x[3]);

  static double Normalize2D(double x[3]);


  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(const double **A, const 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 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 void RandomSeed(long s);  

  static long GetSeed();
  
  static double Random();  

  static double Random(double min, double max);


  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 double* SolveCubic(double c0, double c1, double c2, double c3);

  static double* SolveQuadratic(double c0, double c1, double c2);

  static double* SolveLinear(double c0, double c1);


  static int SolveCubic(double c0, double c1, double c2, double c3, 
                        double *r1, double *r2, double *r3, int *num_roots);


  static int SolveQuadratic(double c0, double c1, double c2, 
                            double *r1, double *r2, int *num_roots);
  
  static int SolveQuadratic( double* c, double* r, int* m );

  static int SolveLinear(double c0, double c1, double *r1, int *num_roots);


  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]);
  static double* RGBToHSV(double r, double g, double b);
  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]);
  static double* HSVToRGB(double h, double s, double v);
  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]);


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


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


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


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


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


  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 int AreBoundsInitialized(double bounds[6]){
    if (bounds[1]-bounds[0]<0.0)
      {
      return 0;
      }
    return 1;
  }


  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 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 int ExtentIsWithinOtherExtent(int extent1[6], int extent2[6]);
  
  static int BoundsIsWithinOtherBounds(double bounds1[6], double bounds2[6], double delta[3]);
  
  static int PointIsWithinBounds(double point[3], double bounds[6], double delta[3]);
  

  static double Inf();
  static double NegInf();
  static double Nan();

protected:
  vtkMath() {};
  ~vtkMath() {};
  
  static long Seed;
private:
  vtkMath(const vtkMath&);  // Not implemented.
  void operator=(const vtkMath&);  // Not implemented.
};

//----------------------------------------------------------------------------
inline vtkTypeInt64 vtkMath::Factorial( int N )
{
  vtkTypeInt64 r = 1;
  while ( N > 1 )
    {
    r *= N--;
    }
  return r;
}

//----------------------------------------------------------------------------
inline int vtkMath::Floor(double x)
{
#if defined i386 || defined _M_IX86
  union { int i[2]; double d; } u;
  // use 52-bit precision of IEEE double to round (x - 0.25) to 
  // the nearest multiple of 0.5, according to prevailing rounding
  // mode which is IEEE round-to-nearest,even
  u.d = (x - 0.25) + 3377699720527872.0; // (2**51)*1.5
  // extract mantissa, use shift to divide by 2 and hence get rid
  // of the bit that gets messed up because the FPU uses
  // round-to-nearest,even mode instead of round-to-nearest,+infinity
  return u.i[0] >> 1;
#else
  return static_cast<int>(floor(x));
#endif
}

//----------------------------------------------------------------------------
inline float vtkMath::Normalize(float x[3])
{
  float den; 
  if ( (den = vtkMath::Norm(x)) != 0.0 )
    {
    for (int i=0; i < 3; i++)
      {
      x[i] /= den;
      }
    }
  return den;
}

//----------------------------------------------------------------------------
inline double vtkMath::Normalize(double x[3])
{
  double den; 
  if ( (den = vtkMath::Norm(x)) != 0.0 )
    {
    for (int i=0; i < 3; i++)
      {
      x[i] /= den;
      }
    }
  return den;
}

//----------------------------------------------------------------------------
inline float vtkMath::Normalize2D(float x[3])
{
  float den; 
  if ( (den = vtkMath::Norm2D(x)) != 0.0 )
    {
    for (int i=0; i < 2; i++)
      {
      x[i] /= den;
      }
    }
  return den;
}

//----------------------------------------------------------------------------
inline double vtkMath::Normalize2D(double x[3])
{
  double den; 
  if ( (den = vtkMath::Norm2D(x)) != 0.0 )
    {
    for (int i=0; i < 2; i++)
      {
      x[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 x[3], 
                                             const float y[3])
{
  return ((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) +
          (x[2]-y[2])*(x[2]-y[2]));
}

//----------------------------------------------------------------------------
inline double vtkMath::Distance2BetweenPoints(const double x[3], 
                                              const double y[3])
{
  return ((x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) +
          (x[2]-y[2])*(x[2]-y[2]));
}

//----------------------------------------------------------------------------
inline double vtkMath::Random(double min, double max)
{
  return (min + vtkMath::Random()*(max-min));
}

//----------------------------------------------------------------------------
// Cross product of two 3-vectors. Result vector in z[3].
inline void vtkMath::Cross(const float x[3], const float y[3], float z[3])
{
  float Zx = x[1]*y[2] - x[2]*y[1]; 
  float Zy = x[2]*y[0] - x[0]*y[2];
  float Zz = x[0]*y[1] - x[1]*y[0];
  z[0] = Zx; z[1] = Zy; z[2] = Zz; 
}

//----------------------------------------------------------------------------
// Cross product of two 3-vectors. Result vector in z[3].
inline void vtkMath::Cross(const double x[3], const double y[3], double z[3])
{
  double Zx = x[1]*y[2] - x[2]*y[1]; 
  double Zy = x[2]*y[0] - x[0]*y[2];
  double Zz = x[0]*y[1] - x[1]*y[0];
  z[0] = Zx; z[1] = Zy; z[2] = Zz; 
}

//BTX
//----------------------------------------------------------------------------
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];
}
//ETX

//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(float A[3][3])
{
  return vtkDeterminant3x3(A);
}

//----------------------------------------------------------------------------
inline double vtkMath::Determinant3x3(double A[3][3])
{
  return vtkDeterminant3x3(A);
}

//----------------------------------------------------------------------------
inline void vtkMath::ClampValue(double *value, const double range[2])
{
  if (value && range)
    {
    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)
    {
    if (value < range[0])
      {
      *clamped_value = range[0];
      }
    else if (value > range[1])
      {
      *clamped_value = range[1];
      }
    else
      {
      *clamped_value = value;
      }
    }
}

#endif

---