progs/util/matrix.c

   1 /*
   2  * matrix.c
   3  *
   4  * Some useful matrix functions.
   5  *
   6  * Brian Paul
   7  * 10 Feb 2004
   8  */
   9
  10
  11
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <math.h>
  15
  16
  17 /**
  18  * Pretty-print the given matrix.
  19  */
  20 void
  21 PrintMatrix(const float p[16])
  22 {
  23    printf("[ %6.3f %6.3f %6.3f %6.3f ]\n", p[0], p[4], p[8], p[12]);
  24    printf("[ %6.3f %6.3f %6.3f %6.3f ]\n", p[1], p[5], p[9], p[13]);
  25    printf("[ %6.3f %6.3f %6.3f %6.3f ]\n", p[2], p[6], p[10], p[14]);
  26    printf("[ %6.3f %6.3f %6.3f %6.3f ]\n", p[3], p[7], p[11], p[15]);
  27 }
  28
  29
  30 /**
  31  * Build a glFrustum matrix.
  32  */
  33 void
  34 Frustum(float left, float right, float bottom, float top, float nearZ, float farZ, float *m)
  35 {
  36    float x = (2.0F*nearZ) / (right-left);
  37    float y = (2.0F*nearZ) / (top-bottom);
  38    float a = (right+left) / (right-left);
  39    float b = (top+bottom) / (top-bottom);
  40    float c = -(farZ+nearZ) / ( farZ-nearZ);
  41    float d = -(2.0F*farZ*nearZ) / (farZ-nearZ);
  42
  43 #define M(row,col)  m[col*4+row]
  44    M(0,0) = x;     M(0,1) = 0.0F;  M(0,2) = a;      M(0,3) = 0.0F;
  45    M(1,0) = 0.0F;  M(1,1) = y;     M(1,2) = b;      M(1,3) = 0.0F;
  46    M(2,0) = 0.0F;  M(2,1) = 0.0F;  M(2,2) = c;      M(2,3) = d;
  47    M(3,0) = 0.0F;  M(3,1) = 0.0F;  M(3,2) = -1.0F;  M(3,3) = 0.0F;
  48 #undef M
  49 }
  50
  51
  52 /**
  53  * Build a glOrtho marix.
  54  */
  55 void
  56 Ortho(float left, float right, float bottom, float top, float nearZ, float farZ, float *m)
  57 {
  58 #define M(row,col)  m[col*4+row]
  59    M(0,0) = 2.0F / (right-left);
  60    M(0,1) = 0.0F;
  61    M(0,2) = 0.0F;
  62    M(0,3) = -(right+left) / (right-left);
  63
  64    M(1,0) = 0.0F;
  65    M(1,1) = 2.0F / (top-bottom);
  66    M(1,2) = 0.0F;
  67    M(1,3) = -(top+bottom) / (top-bottom);
  68
  69    M(2,0) = 0.0F;
  70    M(2,1) = 0.0F;
  71    M(2,2) = -2.0F / (farZ-nearZ);
  72    M(2,3) = -(farZ+nearZ) / (farZ-nearZ);
  73
  74    M(3,0) = 0.0F;
  75    M(3,1) = 0.0F;
  76    M(3,2) = 0.0F;
  77    M(3,3) = 1.0F;
  78 #undef M
  79 }
  80
  81
  82 /**
  83  * Decompose a projection matrix to determine original glFrustum or
  84  * glOrtho parameters.
  85  */
  86 void
  87 DecomposeProjection( const float *m,
  88                      int *isPerspective,
  89                      float *leftOut, float *rightOut,
  90                      float *botOut, float *topOut,
  91                      float *nearOut, float *farOut)
  92 {
  93    if (m[15] == 0.0) {
  94       /* perspective */
  95       float p[16];
  96       const float x = m[0];  /* 2N / (R-L) */
  97       const float y = m[5];  /* 2N / (T-B) */
  98       const float a = m[8];  /* (R+L) / (R-L) */
  99       const float b = m[9];  /* (T+B) / (T-B) */
 100       const float c = m[10]; /* -(F+N) / (F-N) */
 101       const float d = m[14]; /* -2FN / (F-N) */
 102
 103       /* These equations found with simple algebra, knowing the arithmetic
 104        * use to set up a typical perspective projection matrix in OpenGL.
 105        */
 106       const float nearZ = -d / (1.0 - c);
 107       const float farZ = (c - 1.0) * nearZ / (c + 1.0);
 108       const float left = nearZ * (a - 1.0) / x;
 109       const float right = 2.0 * nearZ / x + left;
 110       const float bottom = nearZ * (b - 1.0) / y;
 111       const float top = 2.0 * nearZ / y + bottom;
 112
 113       *isPerspective = 1;
 114       *leftOut = left;
 115       *rightOut = right;
 116       *botOut = bottom;
 117       *topOut = top;
 118       *nearOut = nearZ;
 119       *farOut = farZ;
 120    }
 121    else {
 122       /* orthographic */
 123       const float x = m[0];  /*  2 / (R-L) */
 124       const float y = m[5];  /*  2 / (T-B) */
 125       const float z = m[10]; /* -2 / (F-N) */
 126       const float a = m[12]; /* -(R+L) / (R-L) */
 127       const float b = m[13]; /* -(T+B) / (T-B) */
 128       const float c = m[14]; /* -(F+N) / (F-N) */
 129       /* again, simple algebra */
 130       const float right  = -(a - 1.0) / x;
 131       const float left   = right - 2.0 / x;
 132       const float top    = -(b - 1.0) / y;
 133       const float bottom = top - 2.0 / y;
 134       const float farZ   = (c - 1.0) / z;
 135       const float nearZ  = farZ + 2.0 / z;
 136
 137       *isPerspective = 0;
 138       *leftOut = left;
 139       *rightOut = right;
 140       *botOut = bottom;
 141       *topOut = top;
 142       *nearOut = nearZ;
 143       *farOut = farZ;
 144    }
 145 }
 146
 147
 148 #if 0
 149 /* test harness */
 150 int
 151 main(int argc, char *argv[])
 152 {
 153    float m[16], p[16];
 154    float l, r, b, t, n, f;
 155    int persp;
 156    int i;
 157
 158 #if 0
 159    l = -.9;
 160    r = 1.2;
 161    b = -0.5;
 162    t = 1.4;
 163    n = 30;
 164    f = 84;
 165    printf("  Frustum(%f, %f, %f, %f, %f, %f\n",l+1, r+1.2, b+.5, t+.3, n, f);
 166    Frustum(l+1, r+1.2, b+.5, t+.3, n, f, p);
 167    DecomposeProjection(p, &persp, &l, &r, &b, &t, &n, &f);
 168    printf("glFrustum(%f, %f, %f, %f, %f, %f)\n",
 169           l, r, b, t, n, f);
 170    PrintMatrix(p);
 171 #else
 172    printf("Ortho(-1, 1, -1, 1, 10, 84)\n");
 173    Ortho(-1, 1, -1, 1, 10, 84, m);
 174    PrintMatrix(m);
 175    DecomposeProjection(m, &persp, &l, &r, &b, &t, &n, &f);
 176    printf("Ortho(%f, %f, %f, %f, %f, %f) %d\n", l, r, b, t, n, f, persp);
 177 #endif
 178
 179    return 0;
 180 }
 181 #endif