Transform.cpp
上传用户:sycq158
上传日期:2008-10-22
资源大小:15361k
文件大小:12k
源码类别:

游戏

开发平台:

Visual C++

  1. /*________________________________________________
  2.  |                                          
  3.  |  跟踪球算法
  4.  |________________________________________________*/
  5. #include "stdafx.h"
  6. #include "Transform.h"
  7. #include <math.h>
  8. #include <GL/gl.h>
  9. #include <GL/glu.h>
  10. //---------------- vector & matrix operation -----------------------
  11. // V = ( x, y, z )
  12. void VectorCopy(VECTOR & V, float x, float y, float z)
  13. {
  14.    V.x = x;
  15.    V.y = y;
  16.    V.z = z;
  17. }
  18. // V = A 
  19. void VectorCopy(VECTOR & V, VECTOR A)
  20. {
  21.    V.x = A.x;
  22.    V.y = A.y;
  23.    V.z = A.z;
  24. }
  25. // V = A + B
  26. void VectorAdd(VECTOR & V, VECTOR A, VECTOR B)
  27. {
  28.   V.x = A.x + B.x;
  29.   V.y = A.y + B.y;
  30.   V.z = A.z + B.z;
  31. }
  32. // V = A - B
  33. void VectorSub(VECTOR & V, VECTOR A, VECTOR B)
  34. {
  35.   V.x = A.x - B.x;
  36.   V.y = A.y - B.y;
  37.   V.z = A.z - B.z;
  38. }
  39. // V = s * V
  40. void VectorScale(VECTOR & V, float s)
  41. {
  42.    V.x = s * V.x;
  43.    V.y = s * V.y;
  44.    V.z = s * V.z;
  45. }
  46. // V = A x B
  47. void VectorCross(VECTOR & V, VECTOR A, VECTOR B)
  48. {
  49.   V.x = A.y*B.z - A.z*B.y;
  50.   V.y = A.z*B.x - A.x*B.z;
  51.   V.z = A.x*B.y - A.y*B.x;
  52. }
  53. // dot = A . B
  54. void VectorDot(float & dot, VECTOR A, VECTOR B)
  55. {
  56.   dot = A.x*B.x + A.y*B.y + A.z*B.z;
  57. }
  58. // len = |V|
  59. void VectorLength(float & len, VECTOR V)
  60. {
  61.    len = sqrt(V.x*V.x + V.y*V.y + V.z*V.z);
  62. }
  63. // V = V/(|V|)
  64. void VectorNormalize(VECTOR & V)
  65. {
  66.   float len;
  67.   len = sqrt(V.x*V.x + V.y*V.y + V.z*V.z);
  68.   if(len<=0.0000001f)
  69.     {
  70.   V.x = 0.0f;
  71.   V.y = 0.0f;
  72.   V.z = 0.0f;
  73.     }
  74.   else
  75.     {
  76.   V.x /= len;
  77.   V.y /= len;
  78.   V.z /= len;
  79.     }
  80. }
  81. // M = M1 * M2
  82. void MultMatrixf( float *M,  float *M1, float *M2,
  83.  int m1Row, int m1Col, int m2Row, int m2Col)
  84. {
  85.   if ( m1Col != m2Row )  return ;
  86.   int size = m1Row * m2Col;
  87.   
  88.   for( int i=0, i1=0; i<size; i+= m2Col, i1+= m1Col)
  89.     for( int j = 0; j < m2Col; j++)
  90.     {
  91.       M[i+j] = 0.0f;
  92.       for( int k=0, k1=0; k<m1Col; k++, k1+=m2Col)
  93.   M[i+j] = M[i+j] + M1[i1+k] * M2[k1+j];
  94.     }                         
  95. }
  96. //------------------  inverse  -------------------------
  97. static int LUDcmp( float *data, int *index, int n)
  98. {
  99.   const float TINY = 1.0e-20f;
  100.   float *vv= new float [n];
  101.   float temp, big;
  102.   int i, ii, j, jj, k, kk, imax = 0;
  103.   for (i = 0, ii= 0; i < n; i++, ii += n)
  104.   {
  105.     big = 0.0f;
  106.     for(j = 0; j < n; j++)
  107.      if((temp = fabs( data[ii+j])) > big)  big = temp;
  108.        
  109. if( big == 0.0f )
  110.     {
  111.       delete vv;
  112.   return 0;
  113.     }
  114.      
  115. vv[i] = 1.0f /big;
  116.   }
  117.   for( j = 0, jj = 0; j < n; j++, jj += n)
  118.   {
  119.      for( i = 0, ii = 0; i <= j; i++, ii += n)
  120.      {
  121.        temp = data[ii+j];
  122.        for( int k = 0, kk = 0; k < i; k++, kk += n)
  123.        temp -= data[ii+k] * data[kk+j];
  124.        data[ii+j] = temp;
  125.      }
  126.      big = 0.0f;
  127.      imax = j;
  128.      for ( i = j+1, ii = i * n; i < n; i++, ii += n)
  129.      {
  130.    temp = data[ii+j];
  131.    for( k = 0, kk = 0; k < j; k++, kk += n)
  132.       temp -= data[ii+k] * data[kk+j];
  133.    data[ii+j] = temp;
  134.    if(( temp = vv[i] * fabs(temp)) > big)
  135.    {
  136.      big = temp;
  137.      imax = i;
  138.         }
  139.      }
  140.      if( j != imax)
  141.      {
  142.    ii = imax * n;
  143.    for ( k = 0; k < n; k++)
  144.    {
  145.    temp = data[ii+k];
  146.    data[ii+k] = data[jj+k];
  147.    data[jj+k] = temp;
  148.    }
  149.    vv[imax] = vv[j];
  150.      }
  151.      index[j] = imax;
  152.      if( data[jj+j] == 0.0f ) data[jj+j] = TINY;
  153.      if( j < n)
  154.      {
  155.    temp = 1.0f/data[jj+j];
  156.    for ( i = j+1, ii = i * n; i < n; i++, ii += n)
  157.    data[ii+j] *= temp;
  158.      }
  159.    }
  160.    delete vv;
  161.    return 1;
  162. }
  163. void LUBksb( float *data, int *index, float *b, int n)
  164. {
  165.   int ii, i, j, k = -1, ip;
  166.   float sum;
  167.   for ( i = 0, ii = 0; i < n; i++, ii += n)
  168.   {
  169.     ip = index[i];
  170.     sum = b[ip];
  171.     b[ip] = b[i];
  172.     if ( k != -1)
  173.       for ( j = k; j < i; j++)    sum -= data[ii+j]*b[j];
  174.     else if ( sum)
  175.       k = i;
  176.     b[i] = sum;
  177.   }
  178.   for( i = n-1, ii = i * n; i >=0; i--, ii -= n)
  179.   {
  180.     sum = b[i];
  181.     for( j = i + 1; j < n; j++)   sum -= data[ii+j] * b[j];
  182.     b[i] = sum / data[ii+i];
  183.   }
  184. }
  185. int InvMatrixf(float *data, int n)
  186. {
  187.   int *index = new int [n];
  188.   if ( LUDcmp( data, index, n) == 0)
  189.   {
  190.      delete index;
  191.      return 0;
  192.   }
  193.   float *b = new float [n];
  194.   float *inva = new float [n*n];
  195.   for ( int i = 0; i < n; i++)
  196.   {
  197.     for( int j = 0; j < n; j++)   b[j] = 0.0f;
  198.     b[i] = 1.0f;
  199.     LUBksb( data, index, b, n);
  200.     int jj;
  201.     for ( j = 0, jj = 0; j < n; j++, jj += n)   inva[jj+i] = b[j];
  202.   }
  203.   for ( i = 0; i < n*n; i++)   data[i] = inva[i];
  204.   delete b;
  205.   delete inva;
  206.   delete index;
  207.   return 1;
  208. }
  209. // transpose matrix
  210. void TransMatrix( float *M, int m, int n)
  211. {
  212.    float *temp = new float [m*n];
  213.    for( int i = 0, ii = 0; i < m; i++, ii += m)
  214.    {
  215.      for( int j = 0, jj = 0; j < n; j++, jj += n)
  216.         temp[ii+j] = M[jj+i];
  217.    }
  218.    for( i = 0; i < m*n; i++)  M[i] = temp[i];
  219.    delete temp;
  220.    return;
  221. }
  222. //-------------------------------------------------------
  223. // 跟踪球算法实现
  224. //-------------------------------------------------------
  225. // 参考 《科学计算可视化算法与系统》 P232 - P233
  226. //-------------------------------------------------------
  227. // 周昆,1998,3,23
  228. void Trackball(int   event,   float radius ,
  229.   float centerX, float centerY, float mouseX, float mouseY,
  230.   float &axisX , float &axisY , float &axisZ, float &angle )
  231. {
  232. //------ TB trackball ----
  233. static TrackballStruct TB;
  234.   
  235. if(event==0)  // TB_START
  236. {
  237. TB.centerX = centerX; 
  238. TB.centerY = centerY; 
  239. TB.radius  = radius ;
  240. TB.pointX1 = mouseX ;   // !
  241. TB.pointY1 = mouseY ;   // !
  242. return;
  243. }
  244. if(event>0)  // TB_MOVE , TB_END
  245. {
  246. float sinTao, cosTao, sinOmiga, cosOmiga, sinSita, cosSita;
  247. TB.centerX = centerX; 
  248. TB.centerY = centerY; 
  249. TB.radius  = radius ;
  250. TB.pointX2 = mouseX ;   // !
  251. TB.pointY2 = mouseY ;   // !
  252. VECTOR  V1, V2, V3; 
  253. float len1, len2, len3;
  254. float dot;
  255. float dx,dy;
  256. // 计算 sinTao, cosTao
  257. dx = TB.pointX2 - TB.pointX1;
  258. dy = TB.pointY2 - TB.pointY1;
  259. VectorCopy(V1, dx, dy, 0.0f);
  260.   
  261. dx = TB.pointX1 - TB.centerX;
  262. dy = TB.pointY1 - TB.centerY;
  263. VectorCopy(V2, dx, dy, 0.0f);
  264. VectorCross(V3, V2, V1);
  265. VectorLength(len1, V1);
  266. VectorLength(len2, V2);
  267. VectorLength(len3, V3);
  268.   
  269. if( len1<=0.00001f  || len2<=0.00001f )
  270. sinTao = 0.0f;
  271. cosTao = 1.0f;
  272. else
  273. {
  274. if( V3.z>=0.0f ) sinTao =  len3/(len1*len2);
  275. else             sinTao = -len3/(len1*len2); 
  276. VectorDot(dot, V1, V2);
  277. cosTao = dot/(len1*len2);
  278. }
  279. // 计算sinOmiga,cosOmiga
  280. sinOmiga  = sqrt(dx*dx + dy*dy);
  281. sinOmiga  = sinOmiga/TB.radius;
  282. if(sinOmiga>1.0f)
  283. sinOmiga=1.0;
  284. //99/12/28
  285. sinOmiga = 0.0f;
  286. cosOmiga=sqrt(1.0-sinOmiga*sinOmiga);
  287. //计算sinSita,cosSita
  288. len1 = sqrt(dx*dx + dy*dy);
  289. if(len1<=0.00001f)
  290. {  
  291. sinSita = 0.0f;
  292. cosSita = 1.0f;
  293. }
  294. else
  295. {
  296. sinSita = dy/len1;
  297. cosSita = dx/len1;
  298. }
  299. // 计算旋转角度,用《角度》表示
  300. dx = TB.pointX2 - TB.pointX1;
  301. dy = TB.pointY2 - TB.pointY1;
  302. len1= sqrt(dx*dx + dy*dy);
  303. angle = (len1/TB.radius)*45.0f;
  304. TB.pointX1 = TB.pointX2;
  305. TB.pointY1 = TB.pointY2;
  306. float M [1][3] = { 0.0f,     0.0f,    0.0f      };
  307. float M0[1][3] = { 0.0f,     0.0f,    0.0f      };
  308. float M1[1][3] = {-sinTao,   cosTao,  0.0f      };
  309. float M2[3][3] = { cosOmiga, 0.0f,   -sinOmiga,
  310. 0.0f,     1.0f,   0.0f,
  311. sinOmiga, 0.0f,   cosOmiga  };
  312. float M3[3][3] = { cosSita,  sinSita, 0.0f,
  313. -sinSita,  cosSita, 0.0f,
  314. 0.0f,     0.0f,    1.0f      };
  315. MultMatrixf( &M0[0][0], &M1[0][0], &M2[0][0], 1, 3, 3, 3);
  316. MultMatrixf( &M [0][0], &M0[0][0], &M3[0][0], 1, 3, 3, 3);
  317. //----- output -----------
  318. axisX = M[0][0];
  319. axisY = M[0][1];
  320. axisZ = M[0][2];
  321. } // end of TB_MOVE 
  322. }
  323. // transform point P(x,y,z) from object coordinate to world coordinate
  324. void ObjectToWorld(float point[3], float result[3])
  325. {
  326. // (1) get current modelview matrix
  327. float temp[4][4];
  328. glGetFloatv(GL_MODELVIEW_MATRIX, &temp[0][0]);
  329. // (2) transpose the modelview matrix
  330. TransMatrix( &temp[0][0], 4, 4);
  331. float a[4][1];
  332. float b[4][1];
  333. a[0][0]=point[0];
  334. a[1][0]=point[1];
  335. a[2][0]=point[2];
  336. a[3][0]=1.0f ;
  337.     
  338. // (3) transform (x,y,z) from object coordinate to world coordinate    
  339.     MultMatrixf( &b[0][0], &temp[0][0], &a[0][0], 4, 4, 4, 1);
  340.     result[0]=b[0][0];
  341.     result[1]=b[1][0];
  342.     result[2]=b[2][0];
  343. }
  344. // transform point P(x,y,z) from world coordinate to object coordinate
  345. void WorldToObject(float point[3], float result[3])
  346. {
  347.   // (1) get current modelview matrix
  348.   float temp[4][4];
  349.   glGetFloatv(GL_MODELVIEW_MATRIX, &temp[0][0]);
  350.   // (2) transpose the modelview matrix
  351.   TransMatrix( &temp[0][0], 4, 4);
  352.   // (3) inverse the modelview matrix
  353.   if( InvMatrixf(&temp[0][0], 4) )
  354.   {
  355. float a[4][1];
  356. float b[4][1];
  357. a[0][0]=point[0];
  358. a[1][0]=point[1];
  359. a[2][0]=point[2];
  360. a[3][0]=1.0f ;
  361.     
  362. // (4) transform (x,y,z) from world coordinate to object coordinate    
  363.     MultMatrixf( &b[0][0], &temp[0][0], &a[0][0], 4, 4, 4, 1);
  364.     result[0]=b[0][0];
  365.     result[1]=b[1][0];
  366.     result[2]=b[2][0];
  367.   }
  368. }
  369. // rotate in world coordinate
  370. void WorldRotate(float point[3], float normal[3], float angle)
  371. {
  372. float pp[3];
  373. float nn[3];
  374.     float tt[3];
  375. float t[3]={0.0f,0.0f,0.0f};
  376. WorldToObject(point, pp);
  377. WorldToObject(normal, nn);
  378. WorldToObject(t, tt);
  379. nn[0]-=tt[0];
  380. nn[1]-=tt[1];
  381. nn[2]-=tt[2];
  382. glTranslatef(pp[0],pp[1],pp[2]);
  383. glRotatef(angle,nn[0],nn[1],nn[2]);
  384. glTranslatef(-pp[0],-pp[1],-pp[2]);
  385. }
  386. // translate in world coordinate
  387. void WorldTranslate(float p[3])
  388. {
  389.   float pp[3];
  390.   float tt[3];
  391.   float t[3]={0.0f,0.0f,0.0f};
  392.   WorldToObject(p, pp);
  393.   WorldToObject(t, tt);
  394.   pp[0]-=tt[0];
  395.   pp[1]-=tt[1];
  396.   pp[2]-=tt[2];
  397.   glTranslatef(pp[0],pp[1],pp[2]);
  398. }
  399. // calculate eye axes vectors
  400. void SetEyeStruct(float eyeX,    float eyeY,    float eyeZ, 
  401.  float centerX, float centerY, float centerZ, 
  402.  float upX,     float upY,     float upZ,
  403.  EYE  & eye)
  404. {
  405.   eye.origin.x=eyeX;
  406.   eye.origin.y=eyeY;
  407.   eye.origin.z=eyeZ;
  408.   eye.center.x=centerX;
  409.   eye.center.y=centerY;
  410.   eye.center.z=centerZ;
  411.   eye.up.x=upX;
  412.   eye.up.y=upY;
  413.   eye.up.z=upZ;
  414.   // 计算眼睛各坐标轴的矢量
  415.   VECTOR v1,v2,v3;
  416.   VectorSub(v1,eye.center,eye.origin);
  417.   VectorNormalize(v1);
  418.   VectorCopy(v2,eye.up);
  419.   VectorCross(v3,v2,v1);
  420.   VectorNormalize(v3);
  421.   VectorCross(v2,v3,v1);
  422.   VectorCopy(eye.axisX,v3);
  423.   VectorCopy(eye.axisY,v2);
  424.   VectorCopy(eye.axisZ,v1);
  425. }
  426. // transform point P(x,y,z) from world coordinate to eye coordinate
  427. void WorldToEye(float point[3], float result[3], EYE eye)
  428. {
  429. float tmp[3][3];
  430. float a[3][1];
  431. float b[3][1];
  432. tmp[0][0]=(eye.axisX).x;
  433. tmp[1][0]=(eye.axisX).y;
  434. tmp[2][0]=(eye.axisX).z;
  435. tmp[0][1]=(eye.axisY).x;
  436. tmp[1][1]=(eye.axisY).y;
  437. tmp[2][1]=(eye.axisY).z;
  438. tmp[0][2]=(eye.axisZ).x;
  439. tmp[1][2]=(eye.axisZ).y;
  440. tmp[2][2]=(eye.axisZ).z;
  441. InvMatrixf(&tmp[0][0],3);
  442. a[0][0]=point[0];
  443. a[1][0]=point[1];
  444. a[2][0]=point[2];
  445. MultMatrixf( &b[0][0], &tmp[0][0], &a[0][0], 3, 3, 3, 1);
  446. result[0]=b[0][0]-eye.origin.x;
  447. result[1]=b[1][0]-eye.origin.y;
  448. result[2]=b[2][0]+eye.origin.z;
  449. }
  450. // transform point P(x,y,z) from eye coordinate to world coordinate
  451. void EyeToWorld(float point[3], float result[3], EYE eye)
  452. {
  453.   VECTOR v1,v2,v3,v4,v5,v6;
  454.   VectorCopy(v1,eye.axisZ);
  455.   VectorCopy(v2,eye.axisY);
  456.   VectorCopy(v3,eye.axisX);
  457.   VectorScale(v3,point[0]);
  458.   VectorScale(v2,point[1]);
  459.   VectorScale(v1,point[2]);
  460.   VectorAdd(v4,v1,v2);
  461.   VectorAdd(v5,v4,v3);
  462.   VectorAdd(v6,eye.origin,v5);
  463.   result[0]=v6.x;
  464.   result[1]=v6.y;
  465.   result[2]=v6.z;
  466. }
  467. // rotate in eye coordinate
  468. void EyeRotate(float point[3], float normal[3], float angle, 
  469.   EYE eye)
  470. {
  471.   float pp[3];
  472.   float nn[3];
  473.   float tt[3];
  474.   float t[3]={0.0f,0.0f,0.0f};
  475.   EyeToWorld(point, pp, eye);
  476.   EyeToWorld(normal, nn, eye);
  477.   EyeToWorld(t, tt, eye);
  478.   nn[0]-=tt[0];
  479.   nn[1]-=tt[1];
  480.   nn[2]-=tt[2];
  481.   WorldRotate(pp, nn, angle);
  482. }
  483. // translate in eye coordinate
  484. void EyeTranslate(float p[3], EYE eye)
  485. {
  486.   float pp[3];
  487.   float tt[3];
  488.   float t[3]={0.0f,0.0f,0.0f};
  489.   EyeToWorld(p, pp, eye);
  490.   EyeToWorld(t, tt, eye);
  491.   pp[0]-=tt[0];
  492.   pp[1]-=tt[1];
  493.   pp[2]-=tt[2];
  494.   WorldTranslate(pp);
  495. }