Dogcows Code - chaz/rasterize/blob - mat.h

   1
   2 /*
   3  * CS5600 University of Utah
   4  * Charles McGarvey
   5  * mcgarvey@eng.utah.edu
   6  */
   7
   8 #ifndef _MAT_H_
   9 #define _MAT_H_
  10
  11 #include "common.h"
  12 #include "vec.h"
  13
  14
  15 /*
  16  * A simple matrix class with column-major storage and notation.
  17  */
  18 struct mat
  19 {
  20     vec_t v[4];
  21 };
  22 typedef struct mat mat_t;
  23
  24 /*
  25  * Initialize a matrix with individual components, row by row.
  26  */
  27 INLINE_MAYBE
  28 void mat_init(mat_t* m, scal_t m11, scal_t m12, scal_t m13, scal_t m14,
  29                         scal_t m21, scal_t m22, scal_t m23, scal_t m24,
  30                         scal_t m31, scal_t m32, scal_t m33, scal_t m34,
  31                         scal_t m41, scal_t m42, scal_t m43, scal_t m44)
  32 {
  33     m->v[0] = vec_new2(m11, m21, m31, m41);
  34     m->v[1] = vec_new2(m12, m22, m32, m42);
  35     m->v[2] = vec_new2(m13, m23, m33, m43);
  36     m->v[3] = vec_new2(m14, m24, m34, m44);
  37 }
  38
  39
  40 /*
  41  * Create a new matrix with individual components, row by row.
  42  */
  43 INLINE_MAYBE
  44 mat_t mat_new(scal_t m11, scal_t m12, scal_t m13, scal_t m14,
  45               scal_t m21, scal_t m22, scal_t m23, scal_t m24,
  46               scal_t m31, scal_t m32, scal_t m33, scal_t m34,
  47               scal_t m41, scal_t m42, scal_t m43, scal_t m44)
  48 {
  49     mat_t m;
  50     mat_init(&m, m11, m12, m13, m14,
  51                  m21, m22, m23, m24,
  52                  m31, m32, m33, m34,
  53                  m41, m42, m43, m44);
  54     return m;
  55 }
  56
  57 /*
  58  * Create a new matrix with four column vectors.
  59  */
  60 INLINE_MAYBE
  61 mat_t mat_new2(vec_t a, vec_t b, vec_t c, vec_t d)
  62 {
  63     mat_t m;
  64     m.v[0] = a;
  65     m.v[1] = b;
  66     m.v[2] = c;
  67     m.v[3] = d;
  68     return m;
  69 }
  70
  71 #define MAT_IDENTITY mat_new(S(1.0), S(0.0), S(0.0), S(0.0), \
  72                              S(0.0), S(1.0), S(0.0), S(0.0), \
  73                              S(0.0), S(0.0), S(1.0), S(0.0), \
  74                              S(0.0), S(0.0), S(0.0), S(1.0))
  75
  76
  77 /*
  78  * Get a column vector (can also access the vector array directly).
  79  */
  80 INLINE_MAYBE
  81 vec_t mat_col(mat_t m, int i)
  82 {
  83     return m.v[i];
  84 }
  85
  86 /*
  87  * Get a row vector.
  88  */
  89 INLINE_MAYBE
  90 vec_t mat_row(mat_t m, int i)
  91 {
  92     switch (i) {
  93     case 0:
  94         return vec_new2(m.v[0].x, m.v[1].x, m.v[2].x, m.v[3].x);
  95     case 1:
  96         return vec_new2(m.v[0].y, m.v[1].y, m.v[2].y, m.v[3].y);
  97     case 2:
  98         return vec_new2(m.v[0].z, m.v[1].z, m.v[2].z, m.v[3].z);
  99     case 3:
 100         return vec_new2(m.v[0].w, m.v[1].w, m.v[2].w, m.v[3].w);
 101     }
 102 }
 103
 104
 105 /*
 106  * Print the matrix to stdout.
 107  */
 108 INLINE_MAYBE
 109 void mat_print(mat_t m)
 110 {
 111     printf("|");
 112     vec_print(mat_row(m, 0));
 113     printf("|\n|");
 114     vec_print(mat_row(m, 1));
 115     printf("|\n|");
 116     vec_print(mat_row(m, 2));
 117     printf("|\n|");
 118     vec_print(mat_row(m, 3));
 119     printf("|");
 120 }
 121
 122
 123 /*
 124  * Multiply two matrices together.
 125  */
 126 INLINE_MAYBE
 127 mat_t mat_mult(mat_t a, mat_t b)
 128 {
 129 #define _DOT(I,J) vec_dot2(mat_row(a,I), mat_col(b,J))
 130     return mat_new(_DOT(0,0), _DOT(0,1), _DOT(0,2), _DOT(0,3),
 131                    _DOT(1,0), _DOT(1,1), _DOT(1,2), _DOT(1,3),
 132                    _DOT(2,0), _DOT(2,1), _DOT(2,2), _DOT(2,3),
 133                    _DOT(3,0), _DOT(3,1), _DOT(3,2), _DOT(3,3));
 134 #undef _DOT
 135 }
 136
 137 /*
 138  * Transform a vector using a matrix.
 139  */
 140 INLINE_MAYBE
 141 vec_t mat_apply(mat_t m, vec_t v)
 142 {
 143     return vec_new2(vec_dot2(v,mat_row(m,0)),
 144                     vec_dot2(v,mat_row(m,1)),
 145                     vec_dot2(v,mat_row(m,2)),
 146                     vec_dot2(v,mat_row(m,3)));
 147 }
 148
 149
 150 /*
 151  * Create a new translate matrix.
 152  */
 153 INLINE_MAYBE
 154 mat_t MAT_TRANSLATE(scal_t x, scal_t y, scal_t z)
 155 {
 156     return mat_new(S(1.0), S(0.0), S(0.0),      x,
 157                    S(0.0), S(1.0), S(0.0),      y,
 158                    S(0.0), S(0.0), S(1.0),      z,
 159                    S(0.0), S(0.0), S(0.0), S(1.0));
 160 }
 161
 162 /*
 163  * Create a new translate matrix from a vector.
 164  */
 165 INLINE_MAYBE
 166 mat_t MAT_TRANSLATE2(vec_t v)
 167 {
 168     return MAT_TRANSLATE(v.x, v.y, v.z);
 169 }
 170
 171 /*
 172  * Create a new scale matrix.
 173  */
 174 INLINE_MAYBE
 175 mat_t MAT_SCALE(scal_t x, scal_t y, scal_t z)
 176 {
 177     return mat_new(x, S(0.0), S(0.0), S(0.0),
 178               S(0.0),      y, S(0.0), S(0.0),
 179               S(0.0), S(0.0),      z, S(0.0),
 180               S(0.0), S(0.0), S(0.0), S(1.0));
 181 }
 182
 183 /*
 184  * Create a new scale matrix from a vector.
 185  */
 186 INLINE_MAYBE
 187 mat_t MAT_SCALE2(vec_t v)
 188 {
 189     return MAT_SCALE(v.x, v.y, v.z);
 190 }
 191
 192 /*
 193  * Create a rotation matrix (around the X axis).
 194  */
 195 INLINE_MAYBE
 196 mat_t MAT_ROTATE_X(scal_t theta)
 197 {
 198     scal_t sin_a = scal_sin(theta);
 199     scal_t cos_a = scal_cos(theta);
 200     return mat_new(S(1.0), S(0.0), S(0.0), S(0.0),
 201                    S(0.0),  cos_a, -sin_a, S(0.0),
 202                    S(0.0),  sin_a,  cos_a, S(0.0),
 203                    S(0.0), S(0.0), S(0.0), S(1.0));
 204 }
 205
 206 /*
 207  * Create a rotation matrix (around the Y axis).
 208  */
 209 INLINE_MAYBE
 210 mat_t MAT_ROTATE_Y(scal_t theta)
 211 {
 212     scal_t sin_a = scal_sin(theta);
 213     scal_t cos_a = scal_cos(theta);
 214     return mat_new(cos_a, S(0.0),  sin_a, S(0.0),
 215                   S(0.0), S(1.0), S(0.0), S(0.0),
 216                   -sin_a, S(0.0),  cos_a, S(0.0),
 217                   S(0.0), S(0.0), S(0.0), S(1.0));
 218 }
 219
 220 /*
 221  * Create a rotation matrix (around the Z axis).
 222  */
 223 INLINE_MAYBE
 224 mat_t MAT_ROTATE_Z(scal_t theta)
 225 {
 226     scal_t sin_a = scal_sin(theta);
 227     scal_t cos_a = scal_cos(theta);
 228     return mat_new(cos_a, -sin_a, S(0.0), S(0.0),
 229                    sin_a,  cos_a, S(0.0), S(0.0),
 230                   S(0.0), S(0.0), S(1.0), S(0.0),
 231                   S(0.0), S(0.0), S(0.0), S(1.0));
 232 }
 233
 234 /*
 235  * Create a rotation matrix (around an arbitrary axis).
 236  */
 237 INLINE_MAYBE
 238 mat_t MAT_ROTATE(scal_t theta, scal_t x, scal_t y, scal_t z)
 239 {
 240     /*
 241      * This code is an implementation of an algorithm described by Glenn
 242      * Murray at http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/
 243      */
 244     vec_t v = vec_normalize(vec_new(x, y, z));
 245     x = v.x;
 246     y = v.y;
 247     z = v.z;
 248     scal_t sin_a = scal_sin(theta);
 249     scal_t cos_a = scal_cos(theta);
 250     scal_t x2 = x * x;
 251     scal_t y2 = y * y;
 252     scal_t z2 = z * z;
 253     return mat_new(x2+(S(1.0)-x2)*cos_a, x*y*(S(1.0)-cos_a)-z*sin_a, x*z*(S(1.0)-cos_a)+y*sin_a, S(0.0),
 254                    x*y*(S(1.0)-cos_a)+z*sin_a, y2+(S(1.0)-y2)*cos_a, y*z*(S(1.0)-cos_a)-y*sin_a, S(0.0),
 255                    x*z*(S(1.0)-cos_a)-y*sin_a, y*z*(S(1.0)-cos_a)+x*sin_a, z2+(S(1.0)-z2)*cos_a, S(0.0),
 256                    S(0.0), S(0.0), S(0.0), S(1.0));
 257 }
 258
 259 /*
 260  * Create a view matrix based on eye, spot, and an up vector.
 261  */
 262 INLINE_MAYBE
 263 mat_t MAT_LOOKAT(vec_t eye, vec_t spot, vec_t up)
 264 {
 265     vec_t f = vec_normalize(vec_sub(spot, eye));
 266     vec_t s = vec_normalize(vec_cross(f, up));
 267     vec_t u = vec_cross(s, f);
 268     return mat_mult(mat_new(s.x,    s.y,    s.z, S(0.0),
 269                             u.x,    u.y,    u.z, S(0.0),
 270                            -f.x,   -f.y,   -f.z, S(0.0),
 271                          S(0.0), S(0.0), S(0.0), S(1.0)), MAT_TRANSLATE2(vec_neg(eye)));
 272 }
 273
 274 /*
 275  * Create a 3D orthogonal projection matrix.
 276  */
 277 INLINE_MAYBE
 278 mat_t MAT_ORTHO(scal_t left, scal_t right,
 279                 scal_t bottom, scal_t top,
 280                 scal_t near, scal_t far)
 281 {
 282     scal_t rml = right - left;
 283     scal_t rpl = right + left;
 284     scal_t tmb = top - bottom;
 285     scal_t tpb = top + bottom;
 286     scal_t fmn = far - near;
 287     scal_t fpn = far + near;
 288     return mat_new(S(2.0)/rml,     S(0.0),      S(0.0), -rpl / rml,
 289                        S(0.0), S(2.0)/tmb,      S(0.0), -tpb / tmb,
 290                        S(0.0),     S(0.0), S(-2.0)/fmn, -fpn / fmn,
 291                        S(0.0),     S(0.0),      S(0.0),     S(1.0));
 292 }
 293
 294 /*
 295  * Create a frustum-based projection matrix.
 296  */
 297 INLINE_MAYBE
 298 mat_t MAT_FRUSTUM(scal_t left, scal_t right,
 299                   scal_t bottom, scal_t top,
 300                   scal_t near, scal_t far)
 301 {
 302     scal_t rml = right - left;
 303     scal_t rpl = right + left;
 304     scal_t tmb = top - bottom;
 305     scal_t tpb = top + bottom;
 306     scal_t fmn = far - near;
 307     scal_t fpn = far + near;
 308     scal_t n2  = near * S(2.0);
 309     return mat_new(n2/rml, S(0.0),  rpl/rml,      S(0.0),
 310                    S(0.0), n2/tmb,  tpb/tmb,      S(0.0),
 311                    S(0.0), S(0.0), -fpn/fmn, -n2*far/fmn,
 312                    S(0.0), S(0.0),  S(-1.0),      S(0.0));
 313 }
 314
 315 /*
 316  * Create a perspective projection matrix.
 317  */
 318 INLINE_MAYBE
 319 mat_t MAT_PERSPECTIVE(scal_t fovy, scal_t aspect, scal_t near, scal_t far)
 320 {
 321     scal_t top    = near * scal_tan(fovy * S(0.5));
 322     scal_t bottom = -top;
 323     scal_t left   = bottom * aspect;
 324     scal_t right  = top * aspect;
 325     return MAT_FRUSTUM(left, right, bottom, top, near, far);
 326 }
 327
 328 /*
 329  * Create a viewport matrix.
 330  */
 331 INLINE_MAYBE
 332 mat_t MAT_VIEWPORT(int x, int y, int w, int h)
 333 {
 334     scal_t xs = (scal_t)x;
 335     scal_t ys = (scal_t)y;
 336     scal_t ws = (scal_t)w * S(0.5);
 337     scal_t hs = (scal_t)h * S(0.5);
 338     return mat_new(ws, S(0.0), S(0.0),  ws+xs,
 339                S(0.0),     hs, S(0.0),  hs+ys,
 340                S(0.0), S(0.0), S(1.0), S(0.0),
 341                S(0.0), S(0.0), S(0.0), S(1.0));
 342 }
 343
 344
 345 #endif // _MAT_H_
 346