added matrix code

diff --git a/src/csgray.c b/src/csgray.c
index 2c60be3..0713c06 100644 (file)
--- a/src/csgray.c
+++ b/src/csgray.c
@@ -1,4 +1,9 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
 #include "csgray.h"
+#include "matrix.h"
 
 enum {
        OB_NULL,
@@ -49,9 +54,8 @@ struct camera {
        float fov;
 };
 
-static camera cam;
-static object *root;
-static float identity[] = {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1};
+static struct camera cam;
+static csg_object *root;
 
 int csg_init(void)
 {
@@ -97,91 +101,119 @@ int csg_save(const char *fname)
        return 0;       /* TODO */
 }
 
+#define OBPTR(o) ((struct object*)(o))
+
 void csg_add_object(csg_object *parent, csg_object *child)
 {
-       if(parent->clist) {
-               parent->ctail->next = child;
-               parent->ctail = child;
+       if(!parent) {
+               parent = root;
+       }
+
+       if(parent->ob.clist) {
+               parent->ob.ctail->next = OBPTR(child);
+               parent->ob.ctail = OBPTR(child);
        } else {
-               parent->clist = parent->ctail = child;
+               parent->ob.clist = parent->ob.ctail = OBPTR(child);
        }
-       child->parent = parent;
+       child->ob.parent = OBPTR(parent);
 }
 
 void csg_remove_object(csg_object *parent, csg_object *child)
 {
-       csg_object *c = parent->clist;
-       while(c->next) {
-               if(c->next == child) {
-                       c->next = child->next;
-                       child->next = 0;
-                       child->parent = 0;
+       csg_object *c;
+
+       if(!parent) {
+               parent = root;
+       }
+
+       c = (csg_object*)parent->ob.clist;
+       while(c->ob.next) {
+               if(c->ob.next == OBPTR(child)) {
+                       c->ob.next = child->ob.next;
+                       child->ob.next = 0;
+                       child->ob.parent = 0;
                        return;
                }
-               c = c->next;
+               c = (csg_object*)c->ob.next;
        }
 }
 
 void csg_free_object(csg_object *o)
 {
-       csg_object *c = o->clist;
+       csg_object *c = (csg_object*)o->ob.clist;
        while(c) {
                csg_object *tmp = c;
-               c = c->next;
+               c = (csg_object*)c->ob.next;
                csg_free_object(tmp);
        }
        free(o);
 }
 
-static void init_object(union csg_object *o)
+static union csg_object *alloc_object(void)
 {
-       o->ob.type = OBJ_NULL;
-       memcpy(o->ob.xform, identity, sizeof identity);
+       csg_object *o;
+
+       if(!(o = calloc(sizeof *o, 1))) {
+               return 0;
+       }
+       o->ob.type = OB_NULL;
+       mat4_identity(o->ob.xform);
 
        csg_emission(o, 0, 0, 0);
        csg_color(o, 1, 1, 1);
        csg_roughness(o, 1);
        csg_opacity(o, 1);
+
+       return o;
 }
 
 csg_object *csg_null(float x, float y, float z)
 {
+       return alloc_object();
+}
+
+csg_object *csg_sphere(float x, float y, float z, float r)
+{
        csg_object *o;
 
-       if(!(o = calloc(sizeof *o, 1))) {
+       if(!(o = alloc_object())) {
                return 0;
        }
-       init_object(o);
-       return o;
-}
 
-csg_object *csg_sphere(float x, float y, float z, float r)
-{
+       o->sph.rad = r;
+       mat4_translation(o->ob.xform, x, y, z);
+       return o;
 }
 
 csg_object *csg_cylinder(float x0, float y0, float z0, float x1, float y1, float z1, float r)
 {
+       return 0;
 }
 
 csg_object *csg_plane(float x, float y, float z, float nx, float ny, float nz)
 {
+       return 0;
 }
 
 csg_object *csg_box(float x, float y, float z, float xsz, float ysz, float zsz)
 {
+       return 0;
 }
 
 
 csg_object *csg_union(csg_object *a, csg_object *b)
 {
+       return 0;
 }
 
 csg_object *csg_intersection(csg_object *a, csg_object *b)
 {
+       return 0;
 }
 
 csg_object *csg_subtraction(csg_object *a, csg_object *b)
 {
+       return 0;
 }
 
 

diff --git a/src/main.c b/src/main.c
index 9ea196e..9837073 100644 (file)
--- a/src/main.c
+++ b/src/main.c
@@ -32,7 +32,7 @@ int main(int argc, char **argv)
        ob = csg_sphere(0, 1, 0, 0.8);
        oc = csg_intersection(oa, ob);
 
-       csg_add_object(oc);
+       csg_add_object(0, oc);
 
        csg_render_image(pixels, width, height);
        save_image(out_fname, pixels, width, height);

diff --git a/src/matrix.c b/src/matrix.c
new file mode 100644 (file)
index 0000000..f1347ee
--- /dev/null
+++ b/src/matrix.c
@@ -0,0 +1,519 @@
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#define M(r, c) ((r << 2) + c)
+
+static const float idmat[] = {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1};
+
+void mat4_identity(float *m)
+{
+       memcpy(m, idmat, sizeof idmat);
+}
+
+void mat4_copy(float *dest, float *src)
+{
+       memcpy(dest, src, 16 * sizeof(float));
+}
+
+void mat4_mul(float *dest, float *a, float *b)
+{
+       int i, j;
+       float tmp[16];
+
+       if(dest == a) {
+               memcpy(tmp, a, sizeof tmp);
+               a = tmp;
+       } else if(dest == b) {
+               memcpy(tmp, b, sizeof tmp);
+               b = tmp;
+       }
+
+       for(i=0; i<4; i++) {
+               for(j=0; j<4; j++) {
+                       dest[M(i, j)] = a[M(i, 0)] * b[M(0, j)] + a[M(i, 1)] * b[M(1, j)] +
+                               a[M(i, 2)] * b[M(2, j)] + a[M(i, 3)] * b[M(3, j)];
+               }
+       }
+}
+
+
+void mat4_xform3(float *vdest, float *m, float *v)
+{
+       float x = m[0] + v[0] + m[4] * v[1] + m[8] * v[2] + m[12];
+       float y = m[1] + v[0] + m[5] * v[1] + m[9] * v[2] + m[13];
+       vdest[2] = m[2] + v[0] + m[6] * v[1] + m[10] * v[2] + m[14];
+       vdest[0] = x;
+       vdest[1] = y;
+}
+
+void mat4_xform4(float *vdest, float *m, float *v)
+{
+       float x = m[0] + v[0] + m[4] * v[1] + m[8] * v[2] + m[12] * v[3];
+       float y = m[1] + v[0] + m[5] * v[1] + m[9] * v[2] + m[13] * v[3];
+       float z = m[2] + v[0] + m[6] * v[1] + m[10] * v[2] + m[14] * v[3];
+       vdest[3] = m[3] + v[0] + m[7] * v[1] + m[11] * v[2] + m[15] * v[3];
+       vdest[0] = x;
+       vdest[1] = y;
+       vdest[2] = z;
+}
+
+
+float *mat4_row(float *m, int row)
+{
+       return m + (row << 2);
+}
+
+float mat4_elem(float *m, int row, int col)
+{
+       return m[M(row, col)];
+}
+
+
+void mat4_upper3x3(float *m)
+{
+       m[3] = m[7] = m[11] = m[12] = m[13] = m[14] = 0.0f;
+       m[15] = 1.0f;
+}
+
+#define swap(a, b)     \
+       do { \
+               float tmp = a; \
+               a = b; \
+               b = tmp; \
+       } while(0)
+
+void mat4_transpose(float *m)
+{
+       swap(m[1], m[4]);
+       swap(m[2], m[8]);
+       swap(m[3], m[12]);
+       swap(m[6], m[9]);
+       swap(m[7], m[13]);
+       swap(m[11], m[14]);
+}
+
+float mat4_determinant(float *m)
+{
+       float d11 =     (m[M(1, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 2)])) -
+                               (m[M(2, 1)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 2)])) +
+                               (m[M(3, 1)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 2)]));
+
+       float d12 =     (m[M(0, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 2)])) -
+                               (m[M(2, 1)] * (m[M(0, 2)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 2)])) +
+                               (m[M(3, 1)] * (m[M(0, 2)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 2)]));
+
+       float d13 =     (m[M(0, 1)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 2)])) -
+                               (m[M(1, 1)] * (m[M(0, 2)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 2)])) +
+                               (m[M(3, 1)] * (m[M(0, 2)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 2)]));
+
+       float d14 =     (m[M(0, 1)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 2)])) -
+                               (m[M(1, 1)] * (m[M(0, 2)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 2)])) +
+                               (m[M(2, 1)] * (m[M(0, 2)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 2)]));
+
+       return m[M(0, 0)] * d11 - m[M(1, 0)] * d12 + m[M(2, 0)] * d13 - m[M(3, 0)] * d14;
+}
+
+void mat4_adjoint(float *res, float *m)
+{
+       int i, j;
+       float cof[16];
+
+       cof[M(0, 0)] =  (m[M(1, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 2)])) -
+                                       (m[M(2, 1)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 2)])) +
+                                       (m[M(3, 1)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 2)]));
+       cof[M(0, 1)] =  (m[M(0, 1)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 2)])) -
+                                       (m[M(2, 1)] * (m[M(0, 2)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 2)])) +
+                                       (m[M(3, 1)] * (m[M(0, 2)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 2)]));
+       cof[M(0, 2)] =  (m[M(0, 1)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 2)])) -
+                                       (m[M(1, 1)] * (m[M(0, 2)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 2)])) +
+                                       (m[M(3, 1)] * (m[M(0, 2)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 2)]));
+       cof[M(0, 3)] =  (m[M(0, 1)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 2)])) -
+                                       (m[M(1, 1)] * (m[M(0, 2)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 2)])) +
+                                       (m[M(2, 1)] * (m[M(0, 2)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 2)]));
+
+       cof[M(1, 0)] =  (m[M(1, 0)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 2)])) -
+                                       (m[M(2, 0)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 2)])) +
+                                       (m[M(3, 0)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 2)]));
+       cof[M(1, 1)] =  (m[M(0, 0)] * (m[M(2, 2)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 2)])) -
+                                       (m[M(2, 0)] * (m[M(0, 2)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 2)])) +
+                                       (m[M(3, 0)] * (m[M(0, 2)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 2)]));
+       cof[M(1, 2)] =  (m[M(0, 0)] * (m[M(1, 2)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 2)])) -
+                                       (m[M(1, 0)] * (m[M(0, 2)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 2)])) +
+                                       (m[M(3, 0)] * (m[M(0, 2)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 2)]));
+       cof[M(1, 3)] =  (m[M(0, 0)] * (m[M(1, 2)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 2)])) -
+                                       (m[M(1, 0)] * (m[M(0, 2)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 2)])) +
+                                       (m[M(2, 0)] * (m[M(0, 2)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 2)]));
+
+       cof[M(2, 0)] =  (m[M(1, 0)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 1)])) -
+                                       (m[M(2, 0)] * (m[M(1, 1)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 1)])) +
+                                       (m[M(3, 0)] * (m[M(1, 1)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 1)]));
+       cof[M(2, 1)] =  (m[M(0, 0)] * (m[M(2, 1)] * m[M(3, 3)] - m[M(2, 3)] * m[M(3, 1)])) -
+                                       (m[M(2, 0)] * (m[M(0, 1)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 1)])) +
+                                       (m[M(3, 0)] * (m[M(0, 1)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 1)]));
+       cof[M(2, 2)] =  (m[M(0, 0)] * (m[M(1, 1)] * m[M(3, 3)] - m[M(1, 3)] * m[M(3, 1)])) -
+                                       (m[M(1, 0)] * (m[M(0, 1)] * m[M(3, 3)] - m[M(0, 3)] * m[M(3, 1)])) +
+                                       (m[M(3, 0)] * (m[M(0, 1)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 1)]));
+       cof[M(2, 3)] =  (m[M(0, 0)] * (m[M(1, 1)] * m[M(2, 3)] - m[M(1, 3)] * m[M(2, 1)])) -
+                                       (m[M(1, 0)] * (m[M(0, 1)] * m[M(2, 3)] - m[M(0, 3)] * m[M(2, 1)])) +
+                                       (m[M(2, 0)] * (m[M(0, 1)] * m[M(1, 3)] - m[M(0, 3)] * m[M(1, 1)]));
+
+       cof[M(3, 0)] =  (m[M(1, 0)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(2, 2)] * m[M(3, 1)])) -
+                                       (m[M(2, 0)] * (m[M(1, 1)] * m[M(3, 2)] - m[M(1, 2)] * m[M(3, 1)])) +
+                                       (m[M(3, 0)] * (m[M(1, 1)] * m[M(2, 2)] - m[M(1, 2)] * m[M(2, 1)]));
+       cof[M(3, 1)] =  (m[M(0, 0)] * (m[M(2, 1)] * m[M(3, 2)] - m[M(2, 2)] * m[M(3, 1)])) -
+                                       (m[M(2, 0)] * (m[M(0, 1)] * m[M(3, 2)] - m[M(0, 2)] * m[M(3, 1)])) +
+                                       (m[M(3, 0)] * (m[M(0, 1)] * m[M(2, 2)] - m[M(0, 2)] * m[M(2, 1)]));
+       cof[M(3, 2)] =  (m[M(0, 0)] * (m[M(1, 1)] * m[M(3, 2)] - m[M(1, 2)] * m[M(3, 1)])) -
+                                       (m[M(1, 0)] * (m[M(0, 1)] * m[M(3, 2)] - m[M(0, 2)] * m[M(3, 1)])) +
+                                       (m[M(3, 0)] * (m[M(0, 1)] * m[M(1, 2)] - m[M(0, 2)] * m[M(1, 1)]));
+       cof[M(3, 3)] =  (m[M(0, 0)] * (m[M(1, 1)] * m[M(2, 2)] - m[M(1, 2)] * m[M(2, 1)])) -
+                                       (m[M(1, 0)] * (m[M(0, 1)] * m[M(2, 2)] - m[M(0, 2)] * m[M(2, 1)])) +
+                                       (m[M(2, 0)] * (m[M(0, 1)] * m[M(1, 2)] - m[M(0, 2)] * m[M(1, 1)]));
+
+       for(i=0; i<4; i++) {
+               for(j=0; j<4; j++) {
+                       float val = j % 2 ? -cof[M(j, i)] : cof[M(j, i)];
+                       res[M(j, i)] = (i % 2) ? -val : val;
+               }
+       }
+}
+
+int mat4_inverse(float *m)
+{
+       int i, j;
+       float adj[16];
+       float det;
+
+       mat4_adjoint(adj, m);
+       det = mat4_determinant(m);
+       if(det == 0.0f) {
+               return -1;
+       }
+
+       for(i=0; i<4; i++) {
+               for(j=0; j<4; j++) {
+                       m[M(j, i)] = adj[M(j, i)] / det;
+               }
+       }
+       return 0;
+}
+
+
+void mat4_translation(float *m, float x, float y, float z)
+{
+       mat4_identity(m);
+       m[12] = x;
+       m[13] = y;
+       m[14] = z;
+}
+
+void mat4_rotation_x(float *m, float angle)
+{
+       float sa, ca;
+
+       mat4_identity(m);
+       sa = sin(angle);
+       ca = cos(angle);
+       m[5] = ca;
+       m[6] = sa;
+       m[9] = -sa;
+       m[10] = ca;
+}
+
+void mat4_rotation_y(float *m, float angle)
+{
+       float sa, ca;
+
+       mat4_identity(m);
+       sa = sin(angle);
+       ca = cos(angle);
+       m[0] = ca;
+       m[2] = -sa;
+       m[8] = sa;
+       m[10] = ca;
+}
+
+void mat4_rotation_z(float *m, float angle)
+{
+       float sa, ca;
+
+       mat4_identity(m);
+       sa = sin(angle);
+       ca = cos(angle);
+       m[0] = ca;
+       m[1] = sa;
+       m[4] = -sa;
+       m[5] = ca;
+}
+
+void mat4_rotation(float *m, float angle, float x, float y, float z)
+{
+       float sa = sin(angle);
+       float ca = cos(angle);
+       float invca = 1.0f - ca;
+       float xsq = x * x;
+       float ysq = y * y;
+       float zsq = z * z;
+
+       mat4_identity(m);
+       m[0] = xsq + (1.0f - xsq) * ca;
+       m[4] = x * y * invca - z * sa;
+       m[8] = x * z * invca + y * sa;
+
+       m[1] = x * y * invca + z * sa;
+       m[5] = ysq + (1.0f - ysq) * ca;
+       m[9] = y * z * invca - x * sa;
+
+       m[2] = x * z * invca - y * sa;
+       m[6] = y * z * invca + x * sa;
+       m[10] = zsq + (1.0f - zsq) * ca;
+}
+
+void mat4_scaling(float *m, float sx, float sy, float sz)
+{
+       mat4_identity(m);
+       m[0] = sx;
+       m[5] = sy;
+       m[10] = sz;
+}
+
+
+void mat4_translate(float *m, float x, float y, float z)
+{
+       float tmp[16];
+       mat4_translation(tmp, x, y, z);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_rotate_x(float *m, float angle)
+{
+       float tmp[16];
+       mat4_rotation_x(tmp, angle);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_rotate_y(float *m, float angle)
+{
+       float tmp[16];
+       mat4_rotation_y(tmp, angle);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_rotate_z(float *m, float angle)
+{
+       float tmp[16];
+       mat4_rotation_z(tmp, angle);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_rotate(float *m, float angle, float x, float y, float z)
+{
+       float tmp[16];
+       mat4_rotation(tmp, angle, x, y, z);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_scale(float *m, float sx, float sy, float sz)
+{
+       float tmp[16];
+       mat4_scaling(tmp, sx, sy, sz);
+       mat4_mul(m, m, tmp);
+}
+
+
+void mat4_pre_translate(float *m, float x, float y, float z)
+{
+       float tmp[16];
+       mat4_translation(tmp, x, y, z);
+       mat4_mul(m, tmp, m);
+}
+
+void mat4_pre_rotate_x(float *m, float angle)
+{
+       float tmp[16];
+       mat4_rotation_x(tmp, angle);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_pre_rotate_y(float *m, float angle)
+{
+       float tmp[16];
+       mat4_rotation_y(tmp, angle);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_pre_rotate_z(float *m, float angle)
+{
+       float tmp[16];
+       mat4_rotation_z(tmp, angle);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_pre_rotate(float *m, float angle, float x, float y, float z)
+{
+       float tmp[16];
+       mat4_rotation(tmp, angle, x, y, z);
+       mat4_mul(m, m, tmp);
+}
+
+void mat4_pre_scale(float *m, float sx, float sy, float sz)
+{
+       float tmp[16];
+       mat4_scaling(tmp, sx, sy, sz);
+       mat4_mul(m, m, tmp);
+}
+
+static void normalize(float *v)
+{
+       float len = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+       if(len != 0.0f) {
+               float s = 1.0f / len;
+               v[0] *= s;
+               v[1] *= s;
+               v[2] *= s;
+       }
+}
+
+static void cross(float *res, float *a, float *b)
+{
+       res[0] = a[1] * b[2] - a[2] * b[1];
+       res[1] = a[2] * b[0] - a[0] * b[2];
+       res[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+void mat4_lookat(float *m, float x, float y, float z, float tx, float ty, float tz, float ux, float uy, float uz)
+{
+       int i;
+       float dir[3], right[3], up[3];
+
+       dir[0] = tx - x;
+       dir[1] = ty - y;
+       dir[2] = tz - z;
+       normalize(dir);
+
+       up[0] = ux;
+       up[1] = uy;
+       up[2] = uz;
+       cross(right, dir, up);
+       normalize(right);
+
+       cross(up, right, dir);
+       normalize(up);
+
+       mat4_identity(m);
+
+       for(i=0; i<3; i++) {
+               m[i] = right[i];
+               m[i + 4] = up[i];
+               m[i + 8] = -dir[i];
+       }
+
+       mat4_translate(m, x, y, z);
+}
+
+void mat4_inv_lookat(float *m, float x, float y, float z, float tx, float ty, float tz, float ux, float uy, float uz)
+{
+       int i;
+       float dir[3], right[3], up[3];
+
+       dir[0] = tx - x;
+       dir[1] = ty - y;
+       dir[2] = tz - z;
+       normalize(dir);
+
+       up[0] = ux;
+       up[1] = uy;
+       up[2] = uz;
+       cross(right, dir, up);
+       normalize(right);
+
+       cross(up, right, dir);
+       normalize(up);
+
+       mat4_identity(m);
+
+       for(i=0; i<3; i++) {
+               m[i << 2] = right[i];
+               m[(i << 2) + 1] = up[i];
+               m[(i << 2) + 2] = -dir[i];
+       }
+
+       mat4_pre_translate(m, -x, -y, -z);
+}
+
+void mat4_ortho(float *m, float left, float right, float bottom, float top, float znear, float zfar)
+{
+       float dx = right - left;
+       float dy = top - bottom;
+       float dz = zfar - znear;
+
+       mat4_identity(m);
+       m[0] = 2.0f / dx;
+       m[5] = 2.0f / dy;
+       m[10] = -2.0f / dz;
+       m[12] = -(right + left) / dx;
+       m[13] = -(top + bottom) / dy;
+       m[14] = -(zfar + znear) / dz;
+}
+
+void mat4_frustum(float *m, float left, float right, float bottom, float top, float znear, float zfar)
+{
+       float dx = right - left;
+       float dy = top - bottom;
+       float dz = zfar - znear;
+
+       memset(m, 0, 16 * sizeof *m);
+       m[0] = 2.0f * znear / dx;
+       m[5] = 2.0f * znear / dy;
+       m[8] = (right + left) / dx;
+       m[9] = (top + bottom) / dy;
+       m[10] = -(zfar + znear) / dz;
+       m[14] = -2.0f * zfar * znear / dz;
+       m[11] = -1.0f;
+}
+
+void mat4_perspective(float *m, float fov, float aspect, float znear, float zfar)
+{
+       float s = 1.0f / tan(fov / 2.0f);
+       float range = znear - zfar;
+
+       memset(m, 0, 16 * sizeof *m);
+       m[0] = s / aspect;
+       m[5] = s;
+       m[10] = (znear + zfar) / range;
+       m[14] = 2.0f * znear * zfar / range;
+       m[11] = -1.0f;
+}
+
+
+void mat4_mirror(float *m, float a, float b, float c, float d)
+{
+       m[0] = 1.0f - 2.0f * a * a;
+       m[5] = 1.0f - 2.0f * b * b;
+       m[10] = 1.0f - 2.0f * c * c;
+       m[15] = 1.0f;
+
+       m[1] = m[4] = -2.0f * a * b;
+       m[2] = m[8] = -2.0f * a * c;
+       m[6] = m[9] = -2.0f * b * c;
+
+       m[12] = -2.0f * a * d;
+       m[13] = -2.0f * b * d;
+       m[14] = -2.0f * c * d;
+
+       m[3] = m[7] = m[11] = 0.0f;
+}
+
+
+void mat4_print(float *m, FILE *fp)
+{
+       int i;
+
+       for(i=0; i<4; i++) {
+               fprintf(fp, "[ %4.4g %4.4g %4.4g %4.4g ]\n", m[0], m[1], m[2], m[3]);
+               m += 4;
+       }
+       fputc('\n', fp);
+}
+

diff --git a/src/matrix.h b/src/matrix.h
new file mode 100644 (file)
index 0000000..193fc39
--- /dev/null
+++ b/src/matrix.h
@@ -0,0 +1,52 @@
+#ifndef MATRIX_H_
+#define MATRIX_H_
+
+#include <stdio.h>
+
+void mat4_identity(float *m);
+void mat4_copy(float *dest, float *src);
+void mat4_mul(float *dest, float *a, float *b);
+
+void mat4_xform3(float *vdest, float *m, float *v);
+void mat4_xform4(float *vdest, float *m, float *v);
+
+float *mat4_row(float *m, int row);
+float mat4_elem(float *m, int row, int col);
+
+void mat4_upper3x3(float *m);
+
+void mat4_transpose(float *m);
+int mat4_inverse(float *m);
+
+void mat4_translation(float *m, float x, float y, float z);
+void mat4_rotation_x(float *m, float angle);
+void mat4_rotation_y(float *m, float angle);
+void mat4_rotation_z(float *m, float angle);
+void mat4_rotation(float *m, float angle, float x, float y, float z);
+void mat4_scaling(float *m, float sx, float sy, float sz);
+
+void mat4_translate(float *m, float x, float y, float z);
+void mat4_rotate_x(float *m, float angle);
+void mat4_rotate_y(float *m, float angle);
+void mat4_rotate_z(float *m, float angle);
+void mat4_rotate(float *m, float angle, float x, float y, float z);
+void mat4_scale(float *m, float sx, float sy, float sz);
+
+void mat4_pre_translate(float *m, float x, float y, float z);
+void mat4_pre_rotate_x(float *m, float angle);
+void mat4_pre_rotate_y(float *m, float angle);
+void mat4_pre_rotate_z(float *m, float angle);
+void mat4_pre_rotate(float *m, float angle, float x, float y, float z);
+void mat4_pre_scale(float *m, float sx, float sy, float sz);
+
+void mat4_lookat(float *m, float x, float y, float z, float tx, float ty, float tz, float ux, float uy, float uz);
+void mat4_inv_lookat(float *m, float x, float y, float z, float tx, float ty, float tz, float ux, float uy, float uz);
+void mat4_ortho(float *m, float left, float right, float bottom, float top, float znear, float zfar);
+void mat4_frustum(float *m, float *left, float right, float bottom, float top, float znear, float zfar);
+void mat4_perspective(float *m, float fov, float aspect, float znear, float zfar);
+
+void mat4_mirror(float *m, float a, float b, float c, float d);
+
+void mat4_print(float *m, FILE *fp);
+
+#endif /* MATRIX_H_ */