simgear/scene/sky/clouds3d/mat16fv.cpp

   1 //------------------------------------------------------------------------------
   2 // File : mat16fv.cpp
   3 //------------------------------------------------------------------------------
   4 // GLVU : Copyright 1997 - 2002
   5 //        The University of North Carolina at Chapel Hill
   6 //------------------------------------------------------------------------------
   7 // Permission to use, copy, modify, distribute and sell this software and its
   8 // documentation for any purpose is hereby granted without fee, provided that
   9 // the above copyright notice appear in all copies and that both that copyright
  10 // notice and this permission notice appear in supporting documentation.
  11 // Binaries may be compiled with this software without any royalties or
  12 // restrictions.
  13 //
  14 // The University of North Carolina at Chapel Hill makes no representations
  15 // about the suitability of this software for any purpose. It is provided
  16 // "as is" without express or implied warranty.
  17
  18 //============================================================================
  19 // $Id$
  20 //----------------------------------------------------------------------------
  21 // mat16fv.cpp : opengl-style float[16] matrix routines.
  22 //
  23 // OPENGL STYLE CAMERA VIEWING MATRIX ROUTINES (PREMULT/ROW VECT, 4x4 MATRIX)
  24 // A POINT IN THE WORLD SPACE CAN BE TRANSFORMED TO A PIXEL ON THE CAMERA
  25 // VIEWPLANE BY TRANSFORMING WITH THE COMPOSITE MATRIX: C = V*P*M
  26 // WHERE V IS THE Viewport XFORM, P IS THE Projection XFORM, AND M IS THE
  27 // Modelview XFORM. The associated inverse matrices are also provided.
  28 // Matrices are of the following form:
  29 //   M[16] : 0  4  8 12
  30 //           1  5  9 13
  31 //           2  6 10 14
  32 //           3  7 11 15
  33 //============================================================================
  34
  35 #include <stdio.h>
  36 #include <math.h>
  37 #include "vec3fv.hpp"
  38
  39
  40 float* Copy16fv(float* A, const float* B) // A=B
  41 {
  42   A[0]=B[0];   A[1]=B[1];   A[2]=B[2];   A[3]=B[3];
  43   A[4]=B[4];   A[5]=B[5];   A[6]=B[6];   A[7]=B[7];
  44   A[8]=B[8];   A[9]=B[9];   A[10]=B[10]; A[11]=B[11];
  45   A[12]=B[12]; A[13]=B[13]; A[14]=B[14]; A[15]=B[15];
  46   return A;
  47 }
  48
  49 float* Mult16fv(float* C, const float* A, const float* B) // C=A*B
  50 {
  51   float tC[16];
  52
  53   tC[0]  = A[0]*B[0] + A[4]*B[1] + A[8]*B[2] + A[12]*B[3];
  54   tC[1]  = A[1]*B[0] + A[5]*B[1] + A[9]*B[2] + A[13]*B[3];
  55   tC[2]  = A[2]*B[0] + A[6]*B[1] + A[10]*B[2] + A[14]*B[3];
  56   tC[3]  = A[3]*B[0] + A[7]*B[1] + A[11]*B[2] + A[15]*B[3];
  57
  58   tC[4]  = A[0]*B[4] + A[4]*B[5] + A[8]*B[6] + A[12]*B[7];
  59   tC[5]  = A[1]*B[4] + A[5]*B[5] + A[9]*B[6] + A[13]*B[7];
  60   tC[6]  = A[2]*B[4] + A[6]*B[5] + A[10]*B[6] + A[14]*B[7];
  61   tC[7]  = A[3]*B[4] + A[7]*B[5] + A[11]*B[6] + A[15]*B[7];
  62
  63   tC[8]  = A[0]*B[8] + A[4]*B[9] + A[8]*B[10] + A[12]*B[11];
  64   tC[9]  = A[1]*B[8] + A[5]*B[9] + A[9]*B[10] + A[13]*B[11];
  65   tC[10] = A[2]*B[8] + A[6]*B[9] + A[10]*B[10] + A[14]*B[11];
  66   tC[11] = A[3]*B[8] + A[7]*B[9] + A[11]*B[10] + A[15]*B[11];
  67
  68   tC[12] = A[0]*B[12] + A[4]*B[13] + A[8]*B[14] + A[12]*B[15];
  69   tC[13] = A[1]*B[12] + A[5]*B[13] + A[9]*B[14] + A[13]*B[15];
  70   tC[14] = A[2]*B[12] + A[6]*B[13] + A[10]*B[14] + A[14]*B[15];
  71   tC[15] = A[3]*B[12] + A[7]*B[13] + A[11]*B[14] + A[15]*B[15];
  72
  73   Copy16fv(C,tC);
  74
  75   return(C);
  76 }
  77
  78 float* Mult16fv3fv(float *NewV, const float* M, const float *V)
  79 {
  80   NewV[0] = M[0]*V[0] + M[4]*V[1] + M[8]*V[2];
  81   NewV[1] = M[1]*V[0] + M[5]*V[1] + M[9]*V[2];
  82   NewV[2] = M[2]*V[0] + M[6]*V[1] + M[10]*V[2];
  83   return(NewV);
  84 }
  85
  86 float* Mult16fv3fvPerspDiv(float *NewV, const float* M, const float *V)
  87 {
  88   float W = M[3]*V[0] + M[7]*V[1] + M[11]*V[2] + M[15];
  89   NewV[0] = (M[0]*V[0] + M[4]*V[1] + M[8]*V[2] + M[12]) / W;
  90   NewV[1] = (M[1]*V[0] + M[5]*V[1] + M[9]*V[2] + M[13]) / W;
  91   NewV[2] = (M[2]*V[0] + M[6]*V[1] + M[10]*V[2] + M[14]) / W;
  92   return(NewV);
  93 }
  94
  95 float* Mult16fv4fv(float *NewV, const float* M, const float *V)
  96 {
  97   NewV[0] = M[0]*V[0] + M[4]*V[1] + M[8]*V[2] + M[12]*V[3];
  98   NewV[1] = M[1]*V[0] + M[5]*V[1] + M[9]*V[2] + M[13]*V[3];
  99   NewV[2] = M[2]*V[0] + M[6]*V[1] + M[10]*V[2] + M[14]*V[3];
 100   NewV[3] = M[3]*V[0] + M[7]*V[1] + M[11]*V[2] + M[15]*V[3];
 101   return(NewV);
 102 }
 103
 104 float* Identity16fv(float* M)
 105 {
 106   M[0]=M[5]=M[10]=M[15]=1;
 107   M[1]=M[2]=M[3]=M[4]=M[6]=M[7]=M[8]=M[9]=M[11]=M[12]=M[13]=M[14]=0;
 108   return(M);
 109 }
 110
 111 float* Transpose16fv(float* M)
 112 {
 113   #define SWAP(a,b,t) (t)=(a);(a)=(b);(b)=(t);
 114   float t;
 115   SWAP(M[1],M[4],t);
 116   SWAP(M[2],M[8],t);
 117   SWAP(M[6],M[9],t);
 118   SWAP(M[3],M[12],t);
 119   SWAP(M[7],M[13],t);
 120   SWAP(M[11],M[14],t);
 121   return(M);
 122 }
 123
 124 float* Rotate16fv(float *M, float DegAng, const float Axis[3])
 125 {
 126   float RadAng = DegAng * 0.0174532f;
 127   float ca=(float)cos(RadAng),
 128         sa=(float)sin(RadAng);
 129   if (Axis[0]==1 && Axis[1]==0 && Axis[2]==0)  // ABOUT X-AXIS
 130   {
 131    M[0]=1; M[4]=0;  M[8]=0;   M[12]=0;
 132    M[1]=0; M[5]=ca; M[9]=-sa; M[13]=0;
 133    M[2]=0; M[6]=sa; M[10]=ca; M[14]=0;
 134    M[3]=0; M[7]=0;  M[11]=0;  M[15]=1;
 135   }
 136   if (Axis[0]==0 && Axis[1]==1 && Axis[2]==0)  // ABOUT Y-AXIS
 137   {
 138    M[0]=ca;  M[4]=0; M[8]=sa;  M[12]=0;
 139    M[1]=0;   M[5]=1; M[9]=0;   M[13]=0;
 140    M[2]=-sa; M[6]=0; M[10]=ca; M[14]=0;
 141    M[3]=0;   M[7]=0; M[11]=0;  M[15]=1;
 142   }
 143   if (Axis[0]==0 && Axis[1]==0 && Axis[2]==1)  // ABOUT Z-AXIS
 144   {
 145    M[0]=ca; M[4]=-sa; M[8]=0;  M[12]=0;
 146    M[1]=sa; M[5]=ca;  M[9]=0;  M[13]=0;
 147    M[2]=0;  M[6]=0;   M[10]=1; M[14]=0;
 148    M[3]=0;  M[7]=0;   M[11]=0; M[15]=1;
 149   }
 150   else                                      // ARBITRARY AXIS
 151   {
 152    float l = Axis[0]*Axis[0]+Axis[1]*Axis[1]+Axis[2]*Axis[2];
 153    float x, y, z;
 154    x=Axis[0],y=Axis[1],z=Axis[2];
 155    if (l > 1.0001f || l < 0.9999f && l!=0)
 156    {
 157      // needs normalization
 158      l=1.0f/(float)sqrt(l);
 159      x*=l; y*=l; z*=l;
 160    }
 161    float x2=x*x, y2=y*y, z2=z*z;
 162    M[0]=x2+ca*(1-x2); M[4]=(x*y)+ca*(-x*y)+sa*(-z); M[8]=(x*z)+ca*(-x*z)+sa*y;
 163    M[1]=(x*y)+ca*(-x*y)+sa*z; M[5]=y2+ca*(1-y2); M[9]=(y*z)+ca*(-y*z)+sa*(-x);
 164    M[2]=(x*z)+ca*(-x*z)+sa*(-y); M[6]=(y*z)+ca*(-y*z)+sa*x; M[10]=z2+ca*(1-z2);
 165    M[12]=M[13]=M[14]=M[3]=M[7]=M[11]=0;
 166    M[15]=1;
 167   }
 168   return(M);
 169 }
 170
 171 float* invRotate16fv(float *M, float DegAng, const float Axis[3])
 172 {
 173   Rotate16fv(M,DegAng,Axis);
 174   Transpose16fv(M);
 175   return(M);
 176 }
 177
 178 float* Scale16fv(float* M, float sx, float sy, float sz)
 179 {
 180   M[0]=sx; M[4]=0;  M[8]=0;   M[12]=0;
 181   M[1]=0;  M[5]=sy; M[9]=0;   M[13]=0;
 182   M[2]=0;  M[6]=0;  M[10]=sz; M[14]=0;
 183   M[3]=0;  M[7]=0;  M[11]=0;  M[15]=1;
 184   return(M);
 185 }
 186
 187 float* invScale16fv(float* M, float sx, float sy, float sz)
 188 {
 189   M[0]=1/sx; M[4]=0;    M[8]=0;     M[12]=0;
 190   M[1]=0;    M[5]=1/sy; M[9]=0;     M[13]=0;
 191   M[2]=0;    M[6]=0;    M[10]=1/sz; M[14]=0;
 192   M[3]=0;    M[7]=0;    M[11]=0;    M[15]=1;
 193   return(M);
 194 }
 195
 196 float* Translate16fv(float* M, float tx, float ty, float tz)
 197 {
 198   M[0]=1; M[4]=0;  M[8]=0;  M[12]=tx;
 199   M[1]=0; M[5]=1;  M[9]=0;  M[13]=ty;
 200   M[2]=0; M[6]=0;  M[10]=1; M[14]=tz;
 201   M[3]=0; M[7]=0;  M[11]=0; M[15]=1;
 202   return(M);
 203 }
 204
 205 float* invTranslate16fv(float* M, float tx, float ty, float tz)
 206 {
 207   M[0]=1; M[4]=0;  M[8]=0;  M[12]=-tx;
 208   M[1]=0; M[5]=1;  M[9]=0;  M[13]=-ty;
 209   M[2]=0; M[6]=0;  M[10]=1; M[14]=-tz;
 210   M[3]=0; M[7]=0;  M[11]=0; M[15]=1;
 211   return(M);
 212 }
 213
 214 float* LookAt(float* M,
 215               const float Eye[3],
 216               const float LookAtPt[3],
 217               const float ViewUp[3])
 218 {
 219   float X[3], Y[3], Z[3];
 220   Subtract3fv(Z,Eye,LookAtPt);  Normalize3fv(Z);
 221   CrossProd3fv(X,ViewUp,Z);     Normalize3fv(X);
 222   CrossProd3fv(Y,Z,X);          Normalize3fv(Y);
 223   M[0]=X[0];  M[4]=X[1];  M[8]=X[2];  M[12]=-DotProd3fv(X,Eye);  // TRANS->ROT
 224   M[1]=Y[0];  M[5]=Y[1];  M[9]=Y[2];  M[13]=-DotProd3fv(Y,Eye);
 225   M[2]=Z[0];  M[6]=Z[1];  M[10]=Z[2]; M[14]=-DotProd3fv(Z,Eye);
 226   M[3]=0;     M[7]=0;     M[11]=0;    M[15]=1;
 227   return(M);
 228 }
 229
 230 float* invLookAt(float* M,
 231                  const float Eye[3],
 232                  const float LookAtPt[3],
 233                  const float ViewUp[3])
 234 {
 235   float X[3], Y[3], Z[3];
 236   Subtract3fv(Z,Eye,LookAtPt);  Normalize3fv(Z);
 237   CrossProd3fv(X,ViewUp,Z);     Normalize3fv(X);
 238   CrossProd3fv(Y,Z,X);          Normalize3fv(Y);
 239   M[0]=X[0];  M[4]=Y[0];  M[8]=Z[0];  M[12]=Eye[0];  // ROT->TRANS
 240   M[1]=X[1];  M[5]=Y[1];  M[9]=Z[1];  M[13]=Eye[1];
 241   M[2]=X[2];  M[6]=Y[2];  M[10]=Z[2]; M[14]=Eye[2];
 242   M[3]=0;     M[7]=0;     M[11]=0;    M[15]=1;
 243   return(M);
 244 }
 245
 246 float* Frustum16fv(float* M, float l, float r, float b, float t,
 247                    float n, float f)
 248 {
 249   M[0]=(2*n)/(r-l); M[4]=0;           M[8]=(r+l)/(r-l);   M[12]=0;
 250   M[1]=0;           M[5]=(2*n)/(t-b); M[9]=(t+b)/(t-b);   M[13]=0;
 251   M[2]=0;           M[6]=0;           M[10]=-(f+n)/(f-n); M[14]=(-2*f*n)/(f-n);
 252   M[3]=0;           M[7]=0;           M[11]=-1;           M[15]=0;
 253   return(M);
 254 }
 255
 256 float* invFrustum16fv(float* M, float l, float r, float b, float t,
 257                       float n, float f)
 258 {
 259   M[0]=(r-l)/(2*n); M[4]=0;           M[8]=0;               M[12]=(r+l)/(2*n);
 260   M[1]=0;           M[5]=(t-b)/(2*n); M[9]=0;               M[13]=(t+b)/(2*n);
 261   M[2]=0;           M[6]=0;           M[10]=0;              M[14]=-1;
 262   M[3]=0;           M[7]=0;           M[11]=-(f-n)/(2*f*n); M[15]=(f+n)/(2*f*n);
 263   return(M);
 264 }
 265
 266 float* Perspective(float* M, float Yfov, float Aspect,
 267                    float Ndist, float Fdist)
 268 {
 269   Yfov *= 0.0174532f;  // CONVERT TO RADIANS
 270   float wT=(float)tan(Yfov*0.5f)*Ndist, wB=-wT;
 271   float wR=wT*Aspect, wL=-wR;
 272   Frustum16fv(M,wL,wR,wB,wT,Ndist,Fdist);
 273   return(M);
 274 }
 275
 276 float* invPerspective(float* M, float Yfov, float Aspect,
 277                       float Ndist, float Fdist)
 278 {
 279   Yfov *= 0.0174532f;  // CONVERT TO RADIANS
 280   float wT=(float)tan(Yfov*0.5f)*Ndist, wB=-wT;
 281   float wR=wT*Aspect, wL=-wR;
 282   invFrustum16fv(M,wL,wR,wB,wT,Ndist,Fdist);
 283   return(M);
 284 }
 285
 286 float* Viewing16fv(
 287   float* M,
 288   const float X[3], const float Y[3], const float Z[3], const float O[3])
 289 {
 290   M[0]=X[0];  M[4]=X[1];   M[8]=X[2]; M[12]=-DotProd3fv(X,O);
 291   M[1]=Y[0];  M[5]=Y[1];   M[9]=Y[2]; M[13]=-DotProd3fv(Y,O);
 292   M[2]=Z[0];  M[6]=Z[1];  M[10]=Z[2]; M[14]=-DotProd3fv(Z,O);
 293   M[3]=0;    M[7]=0;    M[11]=0;   M[15]=1;
 294   return(M);
 295 }
 296
 297 // THE INVERSE OF Viewing16fv,
 298 // THIS TAKES A VIEW MATRIX AND RETURNS VIEWING AXES.
 299 // MATRIX ASSUMED TO BE ORTHONORMAL
 300 void Viewing2CoordFrame16fv(
 301   const float *M, float X[3], float Y[3], float Z[3], float O[3])
 302 {
 303   X[0]=M[0];  X[1]=M[4];   X[2]=M[8]; O[0]=-DotProd3fv(M,M+12);
 304   Y[0]=M[1];  Y[1]=M[5];   Y[2]=M[9]; O[1]=-DotProd3fv(M+4,M+12);
 305   Z[0]=M[2];  Z[1]=M[6];  Z[2]=M[10]; O[2]=-DotProd3fv(M+8,M+12);
 306 };
 307
 308 float* invViewing16fv(
 309   float* M,
 310   const float X[3], const float Y[3], const float Z[3], const float O[3])
 311 {
 312   M[0]=X[0];  M[4]=Y[0];   M[8]=Z[0]; M[12]=O[0];
 313   M[1]=X[1];  M[5]=Y[1];   M[9]=Z[1]; M[13]=O[1];
 314   M[2]=X[2];  M[6]=Y[2];  M[10]=Z[2]; M[14]=O[2];
 315   M[3]=0;    M[7]=0;    M[11]=0;   M[15]=1;
 316   return(M);
 317 }
 318
 319 float* Viewport16fv(float* M, int WW, int WH)
 320 {
 321   float WW2=(float)WW*0.5f, WH2=(float)WH*0.5f;
 322   M[0]=WW2;  M[4]=0;     M[8]=0;     M[12]=WW2;
 323   M[1]=0;    M[5]=WH2;   M[9]=0;     M[13]=WH2;
 324   M[2]=0;    M[6]=0;     M[10]=0.5f; M[14]=0.5f;
 325   M[3]=0;    M[7]=0;     M[11]=0;    M[15]=1;
 326   return(M);
 327 }
 328
 329 float* invViewport16fv(float* M, int WW, int WH)
 330 {
 331   float WW2=2.0f/(float)WW, WH2=2.0f/(float)WH;
 332   M[0]=WW2;  M[4]=0;     M[8]=0;    M[12]=-1.0;
 333   M[1]=0;    M[5]=WH2;   M[9]=0;    M[13]=-1.0;
 334   M[2]=0;    M[6]=0;     M[10]=2.0; M[14]=-1.0;
 335   M[3]=0;    M[7]=0;     M[11]=0;   M[15]=1;
 336   return(M);
 337 }
 338
 339 //--------------------------------------------------------------------------
 340 // Given the coefficient [A B C D] of a plane in the implicit form
 341 // Ax+By+Cz+D=0 (see plane.hpp), this routine generates a reflection matrix
 342 // that will "reflect" all points/vectors about the given plane.
 343 // NOTE: the plane is assumed to be normalized: normal vector (A,B,C) is
 344 // normalized (unit-length), and D is the negative distance from the origin
 345 // to the plane along the normal.
 346 //--------------------------------------------------------------------------
 347 float* PlanarReflection16fv(float M[16], const float P[4])
 348 {
 349   float AA=P[0]*P[0], AB=P[0]*P[1], AC=P[0]*P[2], AD=P[0]*P[3],
 350         BB=P[1]*P[1], BC=P[1]*P[2], BD=P[1]*P[3],
 351         CC=P[2]*P[2], CD=P[2]*P[3];
 352   M[0]=1-2*AA;  M[4]=-2*AB;   M[8]=-2*AC;   M[12]=-2*AD;
 353   M[1]=-2*AB;   M[5]=1-2*BB;  M[9]=-2*BC;   M[13]=-2*BD;
 354   M[2]=-2*AC;   M[6]=-2*BC;   M[10]=1-2*CC; M[14]=-2*CD;
 355   M[3]=0;       M[7]=0;       M[11]=0;      M[15]=1;
 356   return(M);
 357 }
 358
 359 //--------------------------------------------------------------------------
 360 // Returns a matrix that will xform a point from the given object-space
 361 // coordinate frame to world space
 362 //   A composite matrix is formed as follows:
 363 //   C = Translation * Rotation * Scaling (using column vectors/pre-multiplication)
 364 //   WorldPt = C * ModelPt
 365 // The corresponding inverse Xform is also provided (GetWorld2ObjXform)
 366 //--------------------------------------------------------------------------
 367 float* Obj2WorldXform16fv(
 368   float *M,
 369   const float X[3], const float Y[3], const float Z[3],
 370   const float O[3], float Scale)
 371 {
 372   float sX[3], sY[3], sZ[3]; // CREATE SCALED VERSION OF ROT AXES
 373   ScalarMult3fv(sX,X,Scale);
 374   ScalarMult3fv(sY,Y,Scale);
 375   ScalarMult3fv(sZ,Z,Scale);
 376   M[0]=sX[0]; M[4]=sY[0]; M[8]=sZ[0];  M[12]=O[0];
 377   M[1]=sX[1]; M[5]=sY[1]; M[9]=sZ[1];  M[13]=O[1];
 378   M[2]=sX[2]; M[6]=sY[2]; M[10]=sZ[2]; M[14]=O[2];
 379   M[3]=0;     M[7]=0;     M[11]=0;     M[15]=1;
 380   return(M);
 381 }
 382
 383 float* World2ObjXform16fv(
 384   float *M,
 385   const float X[3], const float Y[3], const float Z[3],
 386   const float O[3], float Scale)
 387 {
 388   if (Scale<=0) { printf("Too small scale!\n"); Scale=1; }
 389   float invScale = 1/Scale;
 390   float sX[3], sY[3], sZ[3];
 391   ScalarMult3fv(sX,X,invScale);
 392   ScalarMult3fv(sY,Y,invScale);
 393   ScalarMult3fv(sZ,Z,invScale);
 394   M[0]=sX[0]; M[4]=sX[1]; M[8]=sX[2];  M[12]=-DotProd3fv(O,sX);
 395   M[1]=sY[0]; M[5]=sY[1]; M[9]=sY[2];  M[13]=-DotProd3fv(O,sY);
 396   M[2]=sZ[0]; M[6]=sZ[1]; M[10]=sZ[2]; M[14]=-DotProd3fv(O,sZ);
 397   M[3]=0;     M[7]=0;     M[11]=0;     M[15]=1;
 398   return(M);
 399 }
 400
 401 // ONLY TRANSLATES, ROTATES, AND SCALES ARE ALLOWED
 402 float XformCoordFrame16fv(
 403   const float *M, float X[3], float Y[3], float Z[3], float O[3])
 404 {
 405   Set3fv(X, X[0]*M[0] + X[1]*M[4] + X[2]*M[8],
 406              X[0]*M[1] + X[1]*M[5] + X[2]*M[9],
 407              X[0]*M[2] + X[1]*M[6] + X[2]*M[10] );
 408   Set3fv(Y, Y[0]*M[0] + Y[1]*M[4] + Y[2]*M[8],
 409              Y[0]*M[1] + Y[1]*M[5] + Y[2]*M[9],
 410              Y[0]*M[2] + Y[1]*M[6] + Y[2]*M[10] );
 411   Set3fv(Z, Z[0]*M[0] + Z[1]*M[4] + Z[2]*M[8],
 412              Z[0]*M[1] + Z[1]*M[5] + Z[2]*M[9],
 413              Z[0]*M[2] + Z[1]*M[6] + Z[2]*M[10] );
 414   Set3fv(O, O[0]*M[0] + O[1]*M[4] + O[2]*M[8] + M[12],
 415              O[0]*M[1] + O[1]*M[5] + O[2]*M[9] + M[13],
 416              O[0]*M[2] + O[1]*M[6] + O[2]*M[10] + M[14] );
 417
 418   // MUST RENORMALIZE AXES TO FIND THE UNIFORM SCALE
 419   float Scale = Length3fv(X);
 420   ScalarDiv3fv(X,Scale);
 421   ScalarDiv3fv(Y,Scale);
 422   ScalarDiv3fv(Z,Scale);
 423
 424   // RETURN UNIFORM SCALING OF AXES (how much coordinate frame was scaled)
 425   return(Scale);
 426 };
 427
 428
 429 //----------------------------------------------------------------------------
 430 // Given a complete definition for a particular view (viewing, projection,
 431 // and viewport), returns the COMPOSITE xform matrix that takes a point in the
 432 // world space to a screen space (pixel) point. The inverse is also provided.
 433 //----------------------------------------------------------------------------
 434 float* Screen2WorldXform16fv(
 435   float* M,
 436   const float X[3], const float Y[3], const float Z[3],  // VIEWING AXES
 437   const float O[3], // VIEWING ORIGIN
 438   float l, float r, float b, float t, float n, float f, // PROJECTION
 439   int WW, int WH) // VIEWPORT
 440 {
 441   // C = InverseModelview * InverseProjection * InverseViewport
 442   float N[16];
 443   invViewing16fv(M,X,Y,Z,O);
 444   invFrustum16fv(N,l,r,b,t,n,f);
 445   Mult16fv(M,M,N);  // M=M*N;
 446   invViewport16fv(N,WW,WH);
 447   Mult16fv(M,M,N);  // M=M*N;
 448   return(M);
 449 }
 450
 451 float* World2ScreenXform16fv(
 452   float* M,
 453   const float X[3], const float Y[3], const float Z[3], // VIEWING AXES
 454   const float O[3], // VIEWING ORIGIN
 455   float l, float r, float b, float t, float n, float f, // PROJECTION
 456   int WW, int WH) // VIEWPORT
 457 {
 458   // C = Viewport * Projection * Modelview
 459   float N[16];
 460   Viewport16fv(M,WW,WH);
 461   Frustum16fv(N,l,r,b,t,n,f);
 462   Mult16fv(M,M,N);  // M=M*N;
 463   Viewing16fv(N,X,Y,Z,O);
 464   Mult16fv(M,M,N);  // M=M*N;
 465   return(M);
 466 }
 467
 468 void Print16fv(const float* M)
 469 {
 470   printf("\n%f %f %f %f\n",M[0],M[4],M[8],M[12]);
 471   printf("%f %f %f %f\n",M[1],M[5],M[9],M[13]);
 472   printf("%f %f %f %f\n",M[2],M[6],M[10],M[14]);
 473   printf("%f %f %f %f\n\n",M[3],M[7],M[11],M[15]);
 474 }
 475