home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch7_7 / libgm / gmmat4.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-07  |  14.2 KB  |  522 lines
  1. // gmMatrix4.cc - 4x4 element matrix class
  2. //
  3. // libgm++: Graphics Math Library
  4. // Ferdi Scheepers and Stephen F May
  5. // 15 June 1994
  6.  
  7. #include "gmMat4.h"
  8. #include "gmVec3.h"
  9. #include "gmVec4.h"
  10.  
  11. // private function: RCD
  12. // - dot product of row i of matrix A and row j of matrix B
  13.  
  14. inline double RCD(const gmMatrix4& A, const gmMatrix4& B, int i, int j)
  15. {
  16.   return A[i][0] * B[0][j] + A[i][1] * B[1][j] + A[i][2] * B[2][j] +
  17.          A[i][3] * B[3][j];
  18. }
  19.  
  20. // private function: MINOR
  21. // - MINOR of M[r][c]; r in {0,1,2,3}-{r0,r1,r2}; c in {0,1,2,3}-{c0,c1,c2}
  22.  
  23. inline double MINOR(const gmMatrix4& M,
  24.                     int r0, int r1, int r2, int c0, int c1, int c2)
  25. {
  26.   return M[r0][c0] * (M[r1][c1] * M[r2][c2] - M[r2][c1] * M[r1][c2]) -
  27.      M[r0][c1] * (M[r1][c0] * M[r2][c2] - M[r2][c0] * M[r1][c2]) +
  28.      M[r0][c2] * (M[r1][c0] * M[r2][c1] - M[r2][c0] * M[r1][c1]);
  29. }
  30.  
  31. // CONSTRUCTORS
  32.  
  33. gmMatrix4::gmMatrix4()
  34. {
  35.   assign(0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0);
  36. }
  37.  
  38. gmMatrix4::gmMatrix4(const gmMatrix4& M)
  39. {
  40.   assign(M[0][0], M[0][1], M[0][2], M[0][3],
  41.      M[1][0], M[1][1], M[1][2], M[1][3],
  42.      M[2][0], M[2][1], M[2][2], M[2][3],
  43.      M[3][0], M[3][1], M[3][2], M[3][3]);
  44. }
  45.  
  46. gmMatrix4::gmMatrix4(double a00, double a01, double a02, double a03,
  47.              double a10, double a11, double a12, double a13,
  48.              double a20, double a21, double a22, double a23,
  49.              double a30, double a31, double a32, double a33)
  50. {
  51.   assign(a00, a01, a02, a03,
  52.          a10, a11, a12, a13,
  53.      a20, a21, a22, a23,
  54.      a30, a31, a32, a33);
  55. }
  56.  
  57. // ASSIGNMENT
  58.  
  59. gmMatrix4& gmMatrix4::assign(double a00, double a01, double a02, double a03,
  60.                      double a10, double a11, double a12, double a13,
  61.                      double a20, double a21, double a22, double a23,
  62.                      double a30, double a31, double a32, double a33)
  63. {
  64.   m_[0][0] = a00; m_[0][1] = a01; m_[0][2] = a02; m_[0][3] = a03;
  65.   m_[1][0] = a10; m_[1][1] = a11; m_[1][2] = a12; m_[1][3] = a13;
  66.   m_[2][0] = a20; m_[2][1] = a21; m_[2][2] = a22; m_[2][3] = a23;
  67.   m_[3][0] = a30; m_[3][1] = a31; m_[3][2] = a32; m_[3][3] = a33;
  68.   return *this;
  69. }
  70.  
  71. gmMatrix4& gmMatrix4::operator =(const gmMatrix4& M)
  72. {
  73.   assign(M[0][0], M[0][1], M[0][2], M[0][3],
  74.      M[1][0], M[1][1], M[1][2], M[1][3],
  75.      M[2][0], M[2][1], M[2][2], M[2][3],
  76.      M[3][0], M[3][1], M[3][2], M[3][3]);
  77.   return *this;
  78. }
  79.  
  80. ////////////////////////////////////////////////////////////////////////////
  81. // OPERATORS
  82.  
  83. gmMatrix4& gmMatrix4::operator +=(const gmMatrix4& M)
  84. {
  85.   m_[0][0] += M[0][0]; m_[0][1] += M[0][1];
  86.   m_[0][2] += M[0][2]; m_[0][3] += M[0][3];
  87.   
  88.   m_[1][0] += M[1][0]; m_[1][1] += M[1][1];
  89.   m_[1][2] += M[1][2]; m_[1][3] += M[1][3];
  90.   
  91.   m_[2][0] += M[2][0]; m_[2][1] += M[2][1];
  92.   m_[2][2] += M[2][2]; m_[2][3] += M[2][3];
  93.   
  94.   m_[3][0] += M[3][0]; m_[3][1] += M[3][1];
  95.   m_[3][2] += M[3][2]; m_[3][3] += M[3][3];
  96.   return *this;
  97. }
  98.  
  99. gmMatrix4& gmMatrix4::operator -=(const gmMatrix4& M)
  100. {
  101.   m_[0][0] -= M[0][0]; m_[0][1] -= M[0][1];
  102.   m_[0][2] -= M[0][2]; m_[0][3] -= M[0][3];
  103.   
  104.   m_[1][0] -= M[1][0]; m_[1][1] -= M[1][1];
  105.   m_[1][2] -= M[1][2]; m_[1][3] -= M[1][3];
  106.   
  107.   m_[2][0] -= M[2][0]; m_[2][1] -= M[2][1];
  108.   m_[2][2] -= M[2][2]; m_[2][3] -= M[2][3];
  109.   
  110.   m_[3][0] -= M[3][0]; m_[3][1] -= M[3][1];
  111.   m_[3][2] -= M[3][2]; m_[3][3] -= M[3][3];
  112.   return *this;
  113. }
  114.  
  115. gmMatrix4& gmMatrix4::operator *=(const gmMatrix4& M)
  116. {
  117.   assign(RCD(*this, M, 0, 0), RCD(*this, M, 0, 1), 
  118.      RCD(*this, M, 0, 2), RCD(*this, M, 0, 3),
  119.      RCD(*this, M, 1, 0), RCD(*this, M, 1, 1),
  120.      RCD(*this, M, 1, 2), RCD(*this, M, 1, 3),
  121.      RCD(*this, M, 2, 0), RCD(*this, M, 2, 1),
  122.      RCD(*this, M, 2, 2), RCD(*this, M, 2, 3),
  123.      RCD(*this, M, 3, 0), RCD(*this, M, 3, 1),
  124.      RCD(*this, M, 3, 2), RCD(*this, M, 3, 3));
  125.   return *this;
  126. }
  127.  
  128. gmMatrix4& gmMatrix4::operator *=(double d)
  129. {
  130.   m_[0][0] *= d; m_[0][1] *= d; m_[0][2] *= d; m_[0][3] *= d;
  131.   m_[1][0] *= d; m_[1][1] *= d; m_[1][2] *= d; m_[1][3] *= d;
  132.   m_[2][0] *= d; m_[2][1] *= d; m_[2][2] *= d; m_[2][3] *= d;
  133.   m_[3][0] *= d; m_[3][1] *= d; m_[3][2] *= d; m_[3][3] *= d;
  134.   return *this;
  135. }
  136.  
  137. gmMatrix4& gmMatrix4::operator /=(double d)
  138. {
  139.   double di = 1 / d;
  140.   m_[0][0] *= di; m_[0][1] *= di; m_[0][2] *= di; m_[0][3] *= di;
  141.   m_[1][0] *= di; m_[1][1] *= di; m_[1][2] *= di; m_[1][3] *= di;
  142.   m_[2][0] *= di; m_[2][1] *= di; m_[2][2] *= di; m_[2][3] *= di;
  143.   m_[3][0] *= di; m_[3][1] *= di; m_[3][2] *= di; m_[3][3] *= di;
  144.   return *this;
  145. }
  146.  
  147. gmMatrix4 gmMatrix4::operator +(const gmMatrix4& M) const
  148. {
  149.   return gmMatrix4(m_[0][0] + M[0][0], m_[0][1] + M[0][1],
  150.            m_[0][2] + M[0][2], m_[0][3] + M[0][3],
  151.            m_[1][0] + M[1][0], m_[1][1] + M[1][1],
  152.            m_[1][2] + M[1][2], m_[1][3] + M[1][3],
  153.            m_[2][0] + M[2][0], m_[2][1] + M[2][1],
  154.            m_[2][2] + M[2][2], m_[2][3] + M[2][3],
  155.            m_[3][0] + M[3][0], m_[3][1] + M[3][1],
  156.            m_[3][2] + M[3][2], m_[3][3] + M[3][3]);
  157. }
  158.  
  159. gmMatrix4 gmMatrix4::operator -(const gmMatrix4& M) const
  160. {
  161.   return gmMatrix4(m_[0][0] - M[0][0], m_[0][1] - M[0][1],
  162.            m_[0][2] - M[0][2], m_[0][3] - M[0][3],
  163.            m_[1][0] - M[1][0], m_[1][1] - M[1][1],
  164.            m_[1][2] - M[1][2], m_[1][3] - M[1][3],
  165.            m_[2][0] - M[2][0], m_[2][1] - M[2][1],
  166.            m_[2][2] - M[2][2], m_[2][3] - M[2][3],
  167.            m_[3][0] - M[3][0], m_[3][1] - M[3][1],
  168.            m_[3][2] - M[3][2], m_[3][3] - M[3][3]);
  169. }
  170.  
  171. gmMatrix4 gmMatrix4::operator -() const
  172. {
  173.   return gmMatrix4(-m_[0][0], -m_[0][1], -m_[0][2], -m_[0][3],
  174.            -m_[1][0], -m_[1][1], -m_[1][2], -m_[1][3],
  175.            -m_[2][0], -m_[2][1], -m_[2][2], -m_[2][3],
  176.            -m_[3][0], -m_[3][1], -m_[3][2], -m_[3][3]);
  177. }
  178.  
  179. gmMatrix4 gmMatrix4::operator *(const gmMatrix4& M) const
  180. {
  181.   return gmMatrix4(RCD(*this, M, 0, 0), RCD(*this, M, 0, 1), 
  182.            RCD(*this, M, 0, 2), RCD(*this, M, 0, 3),
  183.            RCD(*this, M, 1, 0), RCD(*this, M, 1, 1),
  184.            RCD(*this, M, 1, 2), RCD(*this, M, 1, 3),
  185.            RCD(*this, M, 2, 0), RCD(*this, M, 2, 1),
  186.            RCD(*this, M, 2, 2), RCD(*this, M, 2, 3),
  187.            RCD(*this, M, 3, 0), RCD(*this, M, 3, 1),
  188.            RCD(*this, M, 3, 2), RCD(*this, M, 3, 3));
  189. }
  190.  
  191. gmMatrix4 gmMatrix4::operator *(double d) const
  192. {
  193.   return gmMatrix4(m_[0][0] * d, m_[0][1] * d, m_[0][2] * d, m_[0][3] * d,
  194.            m_[1][0] * d, m_[1][1] * d, m_[1][2] * d, m_[1][3] * d,
  195.            m_[2][0] * d, m_[2][1] * d, m_[2][2] * d, m_[2][3] * d,
  196.            m_[3][0] * d, m_[3][1] * d, m_[3][2] * d, m_[3][3] * d);
  197. }
  198.  
  199. gmMatrix4 gmMatrix4::operator /(double d) const
  200. {
  201.   assert(!gmIsZero(d));
  202.   double di = 1 / d;
  203.   return gmMatrix4(m_[0][0] * di, m_[0][1] * di, m_[0][2] * di, m_[0][3] * di,
  204.            m_[1][0] * di, m_[1][1] * di, m_[1][2] * di, m_[1][3] * di,
  205.            m_[2][0] * di, m_[2][1] * di, m_[2][2] * di, m_[2][3] * di,
  206.            m_[3][0] * di, m_[3][1] * di, m_[3][2] * di, m_[3][3] * di);
  207. }
  208.  
  209. gmMatrix4 operator *(double d, const gmMatrix4& M)
  210. {
  211.   return gmMatrix4(M[0][0] * d, M[0][1] * d, M[0][2] * d, M[0][3] * d,
  212.            M[1][0] * d, M[1][1] * d, M[1][2] * d, M[1][3] * d,
  213.            M[2][0] * d, M[2][1] * d, M[2][2] * d, M[2][3] * d,
  214.            M[3][0] * d, M[3][1] * d, M[3][2] * d, M[3][3] * d);
  215. }
  216.  
  217. bool gmMatrix4::operator ==(const gmMatrix4& M) const
  218. {
  219.   return (gmFuzEQ(m_[0][0], M[0][0]) && gmFuzEQ(m_[0][1], M[0][1]) &&
  220.       gmFuzEQ(m_[0][2], M[0][2]) && gmFuzEQ(m_[0][3], M[0][3]) &&
  221.      
  222.       gmFuzEQ(m_[1][0], M[1][0]) && gmFuzEQ(m_[1][1], M[1][1]) &&
  223.       gmFuzEQ(m_[1][2], M[1][2]) && gmFuzEQ(m_[1][3], M[1][3]) &&
  224.  
  225.        gmFuzEQ(m_[2][0], M[2][0]) && gmFuzEQ(m_[2][1], M[2][1]) &&
  226.       gmFuzEQ(m_[2][2], M[2][2]) && gmFuzEQ(m_[2][3], M[2][3]) &&
  227.  
  228.       gmFuzEQ(m_[3][0], M[3][0]) && gmFuzEQ(m_[3][1], M[3][1]) &&
  229.       gmFuzEQ(m_[3][2], M[3][2]) && gmFuzEQ(m_[3][3], M[3][3]));
  230. }
  231.  
  232. bool gmMatrix4::operator !=(const gmMatrix4& M) const
  233. {
  234.   return (!(*this == M));
  235. }
  236.  
  237. gmVector4 gmMatrix4::operator *(const gmVector4& v) const
  238. {
  239.   return gmVector4(
  240.     m_[0][0] * v[0] + m_[0][1] * v[1] + m_[0][2] * v[2] + m_[0][3] * v[3],
  241.     m_[1][0] * v[0] + m_[1][1] * v[1] + m_[1][2] * v[2] + m_[1][3] * v[3],
  242.     m_[2][0] * v[0] + m_[2][1] * v[1] + m_[2][2] * v[2] + m_[2][3] * v[3],
  243.     m_[3][0] * v[0] + m_[3][1] * v[1] + m_[3][2] * v[2] + m_[3][3] * v[3]);
  244. }
  245.  
  246. gmVector4 operator *(const gmVector4& v, const gmMatrix4& M)
  247. {
  248.   return gmVector4(
  249.     v[0] * M[0][0] + v[1] * M[1][0] + v[2] * M[2][0] + v[3] * M[3][0],
  250.     v[0] * M[0][1] + v[1] * M[1][1] + v[2] * M[2][1] + v[3] * M[3][1],
  251.     v[0] * M[0][2] + v[1] * M[1][2] + v[2] * M[2][2] + v[3] * M[3][2],
  252.     v[0] * M[0][3] + v[1] * M[1][3] + v[2] * M[2][3] + v[3] * M[3][3]);
  253. }
  254.  
  255. // OPERATIONS
  256.  
  257. gmMatrix4 gmMatrix4::transpose() const
  258. {
  259.   return gmMatrix4(m_[0][0], m_[1][0], m_[2][0], m_[3][0],
  260.            m_[0][1], m_[1][1], m_[2][1], m_[3][1],
  261.            m_[0][2], m_[1][2], m_[2][2], m_[3][2],
  262.            m_[0][3], m_[1][3], m_[2][3], m_[3][3]);
  263. }
  264.  
  265. gmMatrix4 gmMatrix4::inverse() const
  266. {
  267.   assert(!isSingular());
  268.   return adjoint() * gmInv(determinant());
  269. }
  270.  
  271. gmMatrix4 gmMatrix4::adjoint() const
  272. {
  273.   gmMatrix4 A;
  274.   
  275.   A[0][0] =  MINOR(*this, 1, 2, 3, 1, 2, 3);
  276.   A[0][1] = -MINOR(*this, 0, 2, 3, 1, 2, 3);
  277.   A[0][2] =  MINOR(*this, 0, 1, 3, 1, 2, 3);
  278.   A[0][3] = -MINOR(*this, 0, 1, 2, 1, 2, 3);
  279.   A[1][0] = -MINOR(*this, 1, 2, 3, 0, 2, 3);
  280.   A[1][1] =  MINOR(*this, 0, 2, 3, 0, 2, 3);
  281.   A[1][2] = -MINOR(*this, 0, 1, 3, 0, 2, 3);
  282.   A[1][3] =  MINOR(*this, 0, 1, 2, 0, 2, 3);
  283.   A[2][0] =  MINOR(*this, 1, 2, 3, 0, 1, 3);
  284.   A[2][1] = -MINOR(*this, 0, 2, 3, 0, 1, 3);
  285.   A[2][2] =  MINOR(*this, 0, 1, 3, 0, 1, 3);
  286.   A[2][3] = -MINOR(*this, 0, 1, 2, 0, 1, 3);
  287.   A[3][0] = -MINOR(*this, 1, 2, 3, 0, 1, 2);
  288.   A[3][1] =  MINOR(*this, 0, 2, 3, 0, 1, 2);
  289.   A[3][2] = -MINOR(*this, 0, 1, 3, 0, 1, 2);
  290.   A[3][3] =  MINOR(*this, 0, 1, 2, 0, 1, 2);
  291.  
  292.   return A;
  293. }
  294.  
  295. double gmMatrix4::determinant() const
  296. {
  297.   return m_[0][0] * MINOR(*this, 1, 2, 3, 1, 2, 3) -
  298.      m_[0][1] * MINOR(*this, 1, 2, 3, 0, 2, 3) +
  299.      m_[0][2] * MINOR(*this, 1, 2, 3, 0, 1, 3) -
  300.      m_[0][3] * MINOR(*this, 1, 2, 3, 0, 1, 2);
  301. }
  302.  
  303. bool gmMatrix4::isSingular() const
  304. {
  305.   return gmIsZero(determinant());
  306. }
  307.  
  308. gmVector3 gmMatrix4::transform(const gmVector3& v) const
  309. {
  310.   return gmVector3(
  311.     v[0] * m_[0][0] + v[1] * m_[1][0] + v[2] * m_[2][0] + m_[3][0],
  312.     v[0] * m_[0][1] + v[1] * m_[1][1] + v[2] * m_[2][1] + m_[3][1],
  313.     v[0] * m_[0][2] + v[1] * m_[1][2] + v[2] * m_[2][2] + m_[3][2]);
  314. }
  315.  
  316. ////////////////////////////////////////////////////////////////////////////
  317. // TRANSFORMATION MATRICES
  318.  
  319. gmMatrix4 gmMatrix4::identity()
  320. {
  321.   return gmMatrix4(1, 0, 0, 0,
  322.            0, 1, 0, 0,
  323.            0, 0, 1, 0,
  324.            0, 0, 0, 1);
  325. }
  326.  
  327. gmMatrix4 gmMatrix4::rotate(double angle, const gmVector3& axis)
  328. {
  329.   gmMatrix4 m_;
  330.   double length = axis.length();
  331.   double a = axis[0] / length;
  332.   double b = axis[1] / length;
  333.   double c = axis[2] / length;
  334.   double aa = a * a;
  335.   double bb = b * b;
  336.   double cc = c * c;
  337.   double sine = sin(gmRadians(angle));
  338.   double cosine = cos(gmRadians(angle));
  339.   double omcos = 1 - cosine;
  340.  
  341.   m_[0][0] = aa + (1 - aa) * cosine;
  342.   m_[1][1] = bb + (1 - bb) * cosine;
  343.   m_[2][2] = cc + (1 - cc) * cosine;
  344.   m_[0][1] = a * b * omcos + c * sine;
  345.   m_[0][2] = a * c * omcos - b * sine;
  346.   m_[1][0] = a * b * omcos - c * sine;
  347.   m_[1][2] = b * c * omcos + a * sine;
  348.   m_[2][0] = a * c * omcos + b * sine;
  349.   m_[2][1] = b * c * omcos - a * sine;
  350.   m_[0][3] = m_[1][3] = m_[2][3] = m_[3][0] = m_[3][1] = m_[3][2] = 0;
  351.   m_[3][3] = 1;
  352.    
  353.   return m_;
  354. }
  355.  
  356. gmMatrix4 gmMatrix4::scale(double x, double y, double z)
  357. {
  358.   return gmMatrix4(x, 0, 0, 0,
  359.            0, y, 0, 0,
  360.            0, 0, z, 0,
  361.            0, 0, 0, 1);
  362. }
  363.  
  364. gmMatrix4 gmMatrix4::translate(double x, double y, double z)
  365. {
  366.   return gmMatrix4(1, 0, 0, 0,
  367.            0, 1, 0, 0,
  368.            0, 0, 1, 0,
  369.            x, y, z, 1);
  370. }
  371.  
  372. ////////////////////////////////////////////////////////////////////////////
  373. // CUBIC BASIS MATRICES
  374.  
  375. gmMatrix4 gmMatrix4::bezierBasis()
  376. {
  377.   return gmMatrix4(-1,  3, -3,  1,
  378.              3, -6,  3,  0,
  379.            -3,  3,  0,  0,
  380.             1,  0,  0,  0);
  381. }
  382.  
  383. gmMatrix4 gmMatrix4::bsplineBasis()
  384. {  
  385.   return gmMatrix4(-1,  3, -3,  1,
  386.             3, -6,  3,  0,
  387.            -3,  0,  3,  0,
  388.             1,  4,  1,  0) / 6;
  389. }
  390.  
  391. gmMatrix4 gmMatrix4::catmullromBasis()
  392. {
  393.   return gmMatrix4(-1,  3, -3,  1,
  394.             2, -5,  4, -1,
  395.            -1,  0,  1,  0,
  396.             0,  2,  0,  0) / 2;
  397. }
  398.  
  399. gmMatrix4 gmMatrix4::hermiteBasis()
  400. {
  401.   return gmMatrix4( 2,  1, -2,  1,
  402.            -3, -2,  3, -1,
  403.             0,  1,  0,  0,
  404.             1,  0,  0,  0);
  405. }
  406.  
  407. gmMatrix4 gmMatrix4::tensedBSplineBasis(double tension)
  408. {
  409.   gmMatrix4 m;
  410.   double sixth = 1.0 / 6.0;
  411.  
  412.   m[0][0] = sixth * (-tension);
  413.   m[0][1] = sixth * (12.0 - 9.0 * tension);
  414.   m[0][2] = sixth * (9.0 * tension - 12.0);
  415.   m[0][3] = sixth * tension;
  416.  
  417.   m[1][0] = sixth * 3.0 * tension;
  418.   m[1][1] = sixth * (12.0 * tension - 18.0);
  419.   m[1][2] = sixth * (18.0 - 15.0 * tension);
  420.   m[1][3] = 0.0;
  421.  
  422.   m[2][0] = sixth * -3.0 * tension;
  423.   m[2][1] = 0.0;
  424.   m[2][2] = sixth * 3.0 * tension;
  425.   m[2][3] = 0.0;
  426.  
  427.   m[3][0] = sixth * tension;
  428.   m[3][1] = sixth * (6.0 - 2.0 * tension);
  429.   m[3][2] = sixth * tension;
  430.   m[3][3] = 0.0;
  431.  
  432.   return m;
  433. }
  434.  
  435. gmMatrix4 gmMatrix4::cardinalBasis(double tension)
  436. {
  437.   gmMatrix4 m;
  438.   
  439.   m[0][0] = -tension;
  440.   m[0][1] = 2.0 - tension;
  441.   m[0][2] = tension - 2.0;
  442.   m[0][3] = tension;
  443.  
  444.   m[1][0] = 2.0 * tension;
  445.   m[1][1] = tension - 3.0;
  446.   m[1][2] = 3 - 2.0 * tension;
  447.   m[1][3] = -tension;
  448.  
  449.   m[2][0] = -tension;
  450.   m[2][1] = 0.0;
  451.   m[2][2] = tension;
  452.   m[2][3] = 0.0;
  453.  
  454.   m[3][0] = 0.0;
  455.   m[3][1] = 1.0;
  456.   m[3][2] = 0.0;
  457.   m[3][3] = 0.0;
  458.   
  459.   return m;
  460. }
  461.  
  462. gmMatrix4 gmMatrix4::tauBasis(double bias, double tension)
  463. {
  464.   gmMatrix4 m;
  465.   
  466.   m[0][0] = tension * (bias - 1.0);
  467.   m[0][1] = 2.0 - tension * bias;
  468.   m[0][2] = tension * (1.0 - bias) - 2.0;
  469.   m[0][3] = tension * bias;
  470.  
  471.   m[1][0] = tension * 2.0 * (1.0 - bias);
  472.   m[1][1] = tension * (3.0 * bias - 1.0) - 3.0;
  473.   m[1][2] = 3.0 - tension;
  474.   m[1][3] = -tension * bias;
  475.  
  476.   m[2][0] = tension * (bias - 1.0);
  477.   m[2][1] = tension * (1.0 - 2.0 * bias);
  478.   m[2][2] = tension * bias;
  479.   m[2][3] = 0.0;
  480.  
  481.   m[3][0] = 0.0;
  482.   m[3][1] = 1.0;
  483.   m[3][2] = 0.0;
  484.   m[3][3] = 0.0;
  485.   
  486.   return m;
  487. }
  488.  
  489. gmMatrix4 gmMatrix4::betaSplineBasis(double bias, double tension)
  490. {
  491.   gmMatrix4 m;
  492.   double bias2, bias3, d;
  493.  
  494.   bias2 = bias * bias;
  495.   bias3 = bias * bias2;
  496.   d = 1.0 / (tension + 2.0 * bias3 + 4.0 * bias2 + 4.0 * bias + 2.0);
  497.  
  498.   m[0][0] = -d * 2.0 * bias3;
  499.   m[0][1] = d * 2.0 * (tension + bias3 + bias2 + bias);
  500.   m[0][2] = -d * 2.0 * (tension + bias2 + bias + 1.0);
  501.   m[0][3] = d * 2.0;
  502.  
  503.   m[1][0] = d * 6.0 * bias3;
  504.   m[1][1] = -d * 3.0 * (tension + 2.0 * bias3 + 2.0 * bias2);
  505.   m[1][2] = d * 3.0 * (tension + 2 * bias2);
  506.   m[1][3] = 0.0;
  507.   
  508.   m[2][0] = -d * 6.0 * bias3;
  509.   m[2][1] = d * 6.0 * (bias3 - bias);
  510.   m[2][2] = d * 6.0 * bias;
  511.   m[2][3] = 0.0;
  512.  
  513.   m[3][0] = d * 2.0 * bias3;
  514.   m[3][1] = d * (tension + 4 * (bias2 + bias));
  515.   m[3][2] = d * 2.0;
  516.   m[3][3] = 0.0;
  517.  
  518.   return m;
  519. }
  520.  
  521.  
  522.