# HG changeset patch # User John Tsiombikas # Date 1340852750 -10800 # Node ID 96de911d05d4ffc58b3d03ba4880348793825b05 # Parent aebcd71cc3cfddc51cf0fb83b87874932ed2da13 started a rough prototype diff -r aebcd71cc3cf -r 96de911d05d4 prototype/Makefile --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/Makefile Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,30 @@ +csrc = $(wildcard src/*.c) $(wildcard vmath/*.c) +ccsrc = $(wildcard src/*.cc) $(wildcard vmath/*.cc) +obj = $(csrc:.c=.o) $(ccsrc:.cc=.o) +dep = $(obj:.o=.d) +bin = proto + +CFLAGS = -pedantic -Wall -g -Ivmath +CXXFLAGS = $(CFLAGS) +LDFLAGS = $(libgl) -lm + +ifeq ($(shell uname -s), Darwin) + libgl = -framework OpenGL -framework GLUT -lglew +else + libgl = -lGL -lGLU -lglut -lGLEW +endif + +$(bin): $(obj) + $(CXX) -o $@ $(obj) $(LDFLAGS) + +-include $(dep) + +%.d: %.c + @$(CPP) $(CFLAGS) $< -MM -MT $(@:.d=.o) >$@ + +%.d: %.cc + @$(CPP) $(CXXFLAGS) $< -MM -MT $(@:.d=.o) >$@ + +.PHONY: clean +clean: + rm -f $(obj) $(bin) $(dep) diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/camera.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/camera.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,177 @@ +#include "opengl.h" +#include "camera.h" + +Camera::Camera() +{ + inval_cache(); +} + +Camera::~Camera() +{ +} + +void Camera::calc_inv_matrix(Matrix4x4 *mat) const +{ + *mat = matrix().inverse(); +} + +void Camera::set_glmat(const Matrix4x4 &mat) const +{ +#ifdef SINGLE_PRECISION_MATH + if(glLoadTransposeMatrixfARB) { + glLoadTransposeMatrixfARB((float*)&mat); + } else { + Matrix4x4 tmat = mat.transposed(); + glLoadMatrixf((float*)&tmat); + } +#else + if(glLoadTransposeMatrixdARB) { + glLoadTransposeMatrixdARB((double*)&mat); + } else { + Matrix4x4 tmat = mat.transposed(); + glLoadMatrixd((double*)&tmat); + } +#endif +} + +const Matrix4x4 &Camera::matrix() const +{ + if(!mcache.valid) { + calc_matrix(&mcache.mat); + mcache.valid = true; + } + return mcache.mat; +} + +const Matrix4x4 &Camera::inv_matrix() const +{ + if(!mcache_inv.valid) { + calc_inv_matrix(&mcache_inv.mat); + mcache_inv.valid = true; + } + return mcache_inv.mat; +} + +void Camera::use() const +{ + set_glmat(matrix()); +} + +void Camera::use_inverse() const +{ + set_glmat(inv_matrix()); +} + +void Camera::input_move(float x, float y, float z) +{ +} + +void Camera::input_rotate(float x, float y, float z) +{ +} + +void Camera::input_zoom(float factor) +{ +} + + +// ---- orbit camera ---- + +OrbitCamera::OrbitCamera() +{ + theta = 0.0; + phi = 0.0; + rad = 10.0; +} + +OrbitCamera::~OrbitCamera() +{ +} + +void OrbitCamera::calc_matrix(Matrix4x4 *mat) const +{ + mat->reset_identity(); + mat->translate(Vector3(0, 0, -rad)); + mat->rotate(Vector3(phi, 0, 0)); + mat->rotate(Vector3(0, theta, 0)); +} + +void OrbitCamera::calc_inv_matrix(Matrix4x4 *mat) const +{ + mat->reset_identity(); + mat->rotate(Vector3(0, theta, 0)); + mat->rotate(Vector3(phi, 0, 0)); + mat->translate(Vector3(0, 0, -rad)); +} + +void OrbitCamera::input_rotate(float x, float y, float z) +{ + theta += x; + phi += y; + + if(phi < -M_PI / 2.0) + phi = -M_PI / 2.0; + if(phi > M_PI / 2.0) + phi = M_PI / 2.0; + + inval_cache(); +} + +void OrbitCamera::input_zoom(float factor) +{ + rad += factor; + if(rad < 0.0) + rad = 0.0; + + inval_cache(); +} + + +FlyCamera::FlyCamera() +{ + pos.z = 10.0f; +} + +void FlyCamera::calc_matrix(Matrix4x4 *mat) const +{ + Matrix3x3 rmat = rot.get_rotation_matrix().transposed(); + Matrix4x4 tmat; + tmat.set_translation(pos); + *mat = tmat * Matrix4x4(rmat); +} + +/*void FlyCamera::calc_inv_matrix(Matrix4x4 *mat) const +{ +}*/ + +const Vector3 &FlyCamera::get_position() const +{ + return pos; +} + +const Quaternion &FlyCamera::get_rotation() const +{ + return rot; +} + +void FlyCamera::input_move(float x, float y, float z) +{ + static const Vector3 vfwd(0, 0, 1), vright(1, 0, 0); + + Vector3 k = vfwd.transformed(rot); + Vector3 i = vright.transformed(rot); + Vector3 j = cross_product(k, i); + + pos += i * x + j * y + k * z; + inval_cache(); +} + +void FlyCamera::input_rotate(float x, float y, float z) +{ + Vector3 axis(x, y, z); + float axis_len = axis.length(); + rot.rotate(axis / axis_len, axis_len); + rot.normalize(); + + inval_cache(); +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/camera.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/camera.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,69 @@ +#ifndef CAMERA_H_ +#define CAMERA_H_ + +#include "vmath/vmath.h" + +class Camera { +protected: + mutable struct { + bool valid; + Matrix4x4 mat; + } mcache, mcache_inv; + + virtual void calc_matrix(Matrix4x4 *mat) const = 0; + virtual void calc_inv_matrix(Matrix4x4 *mat) const; + + void inval_cache() { mcache.valid = mcache_inv.valid = false; } + void set_glmat(const Matrix4x4 &m) const; + +public: + Camera(); + virtual ~Camera(); + + const Matrix4x4 &matrix() const; + const Matrix4x4 &inv_matrix() const; + + void use() const; + void use_inverse() const; + + // these do nothing, override to provide input handling + virtual void input_move(float x, float y, float z); + virtual void input_rotate(float x, float y, float z); + virtual void input_zoom(float factor); +}; + +class OrbitCamera : public Camera { +private: + float theta, phi, rad; + + void calc_matrix(Matrix4x4 *mat) const; + void calc_inv_matrix(Matrix4x4 *mat) const; + +public: + OrbitCamera(); + virtual ~OrbitCamera(); + + void input_rotate(float x, float y, float z); + void input_zoom(float factor); +}; + +class FlyCamera : public Camera { +private: + Vector3 pos; + Quaternion rot; + + void calc_matrix(Matrix4x4 *mat) const; + //void calc_inv_matrix(Matrix4x4 *mat) const; + +public: + FlyCamera(); + + const Vector3 &get_position() const; + const Quaternion &get_rotation() const; + + void input_move(float x, float y, float z); + void input_rotate(float x, float y, float z); +}; + + +#endif // CAMERA_H_ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/level.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/level.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,114 @@ +#include +#include +#include "opengl.h" +#include "level.h" +#include "tile.h" + +Level::Level() +{ + cell_size = 1.0; + + cells = 0; + xsz = ysz = 0; +} + +Level::~Level() +{ + delete [] cells; +} + +bool Level::load(const char *fname) +{ + xsz = ysz = 64; + + cells = new GridCell*[xsz * ysz]; + memset(cells, 0, xsz * ysz * sizeof *cells); + + GridCell *g = new GridCell; + g->add_tile(new Tile); + cells[ysz / 2 * xsz + xsz / 2] = g; + + return true; +} + +bool Level::save(const char *fname) const +{ + return false; +} + +const GridCell *Level::get_cell(int x, int y) const +{ + if(x < 0 || x >= xsz || y < 0 || y >= ysz) { + return 0; + } + return cells[y * xsz + x]; +} + +Vector3 Level::get_cell_pos(int x, int y) const +{ + float posx = (x - xsz / 2) * cell_size; + float posy = (y - ysz / 2) * cell_size; + return Vector3(posx, 0, posy); +} + +void Level::draw() const +{ + glMatrixMode(GL_MODELVIEW); + + draw_grid(); + + for(int i=0; idraw(); + glPopMatrix(); + } + } + } +} + +void Level::draw_grid() const +{ + float xlen = xsz * cell_size; + float ylen = ysz * cell_size; + + glPushAttrib(GL_ENABLE_BIT); + glDisable(GL_LIGHTING); + + glBegin(GL_LINES); + glColor3f(0.4, 0.4, 0.4); + + float y = -ylen / 2.0 - cell_size / 2.0; + for(int i=0; idraw(); + } +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/level.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/level.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,44 @@ +#ifndef LEVEL_H_ +#define LEVEL_H_ + +#include +#include "vmath/vmath.h" + +class GridCell; +class Tile; + +class Level { +private: + // cells are stored as a linear array of pointers to GridCells + // null pointers mean unpopulated cells. + GridCell **cells; + int xsz, ysz; + float cell_size; + + void draw_grid() const; + +public: + Level(); + ~Level(); + + bool load(const char *fname); + bool save(const char *fname) const; + + const GridCell *get_cell(int x, int y) const; + Vector3 get_cell_pos(int x, int y) const; + + void draw() const; +}; + +class GridCell { +private: + // each grid-cell might contain multiple tiles. + std::vector tiles; + +public: + void add_tile(const Tile *tile); + + void draw() const; +}; + +#endif // LEVEL_H_ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/main.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/main.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,105 @@ +#include +#include +#include +#include "opengl.h" +#include "level.h" +#include "camera.h" + +void disp(); +void reshape(int x, int y); +void keyb(unsigned char key, int x, int y); +void mouse(int bn, int state, int x, int y); +void motion(int x, int y); + +static Level *level; +static OrbitCamera cam; + +int main(int argc, char **argv) +{ + glutInit(&argc, argv); + glutInitWindowSize(800, 600); + glutInitDisplayMode(GLUT_RGB | GLUT_DEPTH | GLUT_DOUBLE | GLUT_MULTISAMPLE); + glutCreateWindow("prototype"); + + glutDisplayFunc(disp); + glutReshapeFunc(reshape); + glutKeyboardFunc(keyb); + glutMouseFunc(mouse); + glutMotionFunc(motion); + + glewInit(); + + glEnable(GL_LIGHTING); + glEnable(GL_LIGHT0); + float ldir[] = {-1, 1, 2, 0}; + glLightfv(GL_LIGHT0, GL_POSITION, ldir); + glEnable(GL_NORMALIZE); + + glEnable(GL_DEPTH_TEST); + glEnable(GL_CULL_FACE); + glEnable(GL_MULTISAMPLE); + + level = new Level; + if(!level->load("foobar")) { + return 1; + } + + glutMainLoop(); +} + +void disp() +{ + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + cam.use(); + + level->draw(); + + glutSwapBuffers(); + assert(glGetError() == GL_NO_ERROR); +} + +void reshape(int x, int y) +{ + glViewport(0, 0, x, y); + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluPerspective(45.0, (float)x / (float)y, 1.0, 1000.0); +} + +void keyb(unsigned char key, int x, int y) +{ + switch(key) { + case 27: + exit(0); + } +} + +static int prev_x, prev_y; +static bool bnstate[32]; + +void mouse(int bn, int state, int x, int y) +{ + prev_x = x; + prev_y = y; + bnstate[bn - GLUT_LEFT_BUTTON] = state == GLUT_DOWN; +} + +void motion(int x, int y) +{ + int dx = x - prev_x; + int dy = y - prev_y; + prev_x = x; + prev_y = y; + + if(bnstate[0]) { + cam.input_rotate(dx * 0.01, dy * 0.01, 0); + glutPostRedisplay(); + } + if(bnstate[2]) { + cam.input_zoom(dy * 0.1); + glutPostRedisplay(); + } +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/opengl.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/opengl.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,12 @@ +#ifndef OPENGL_H_ +#define OPENGL_H_ + +#include + +#ifndef __APPLE__ +#include +#else +#include +#endif + +#endif /* OPENGL_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/tile.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/tile.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,23 @@ +#include +#include "opengl.h" +#include "tile.h" + +bool Tile::load(const char *fname) +{ + return true; +} + +void Tile::draw() const +{ + float color[] = {1, 0, 0, 1}; + float white[] = {1, 1, 1, 1}; + + glPushAttrib(GL_LIGHTING_BIT); + + glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, color); + glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, white); + glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 60.0); + glutSolidSphere(0.5, 16, 8); + + glPopAttrib(); +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/src/tile.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/src/tile.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,11 @@ +#ifndef TILE_H_ +#define TILE_H_ + +class Tile { +public: + bool load(const char *fname); + + void draw() const; +}; + +#endif // TILE_H_ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/basis.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/basis.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,63 @@ +#include "basis.h" +#include "vmath.h" + +Basis::Basis() +{ + i = Vector3(1, 0, 0); + j = Vector3(0, 1, 0); + k = Vector3(0, 0, 1); +} + +Basis::Basis(const Vector3 &i, const Vector3 &j, const Vector3 &k) +{ + this->i = i; + this->j = j; + this->k = k; +} + +Basis::Basis(const Vector3 &dir, bool left_handed) +{ + k = dir; + j = Vector3(0, 1, 0); + i = cross_product(j, k); + j = cross_product(k, i); +} + +/** Rotate with euler angles */ +void Basis::rotate(scalar_t x, scalar_t y, scalar_t z) { + Matrix4x4 RotMat; + RotMat.set_rotation(Vector3(x, y, z)); + i.transform(RotMat); + j.transform(RotMat); + k.transform(RotMat); +} + +/** Rotate by axis-angle specification */ +void Basis::rotate(const Vector3 &axis, scalar_t angle) { + Quaternion q; + q.set_rotation(axis, angle); + i.transform(q); + j.transform(q); + k.transform(q); +} + +/** Rotate with a 4x4 matrix */ +void Basis::rotate(const Matrix4x4 &mat) { + i.transform(mat); + j.transform(mat); + k.transform(mat); +} + +/** Rotate with a quaternion */ +void Basis::rotate(const Quaternion &quat) { + i.transform(quat); + j.transform(quat); + k.transform(quat); +} + +/** Extract a rotation matrix from the orthogonal basis */ +Matrix3x3 Basis::create_rotation_matrix() const { + return Matrix3x3( i.x, j.x, k.x, + i.y, j.y, k.y, + i.z, j.z, k.z); +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/basis.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/basis.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,57 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_BASIS_H_ +#define VMATH_BASIS_H_ + +#include "vector.h" +#include "matrix.h" + +enum { + LEFT_HANDED, + RIGHT_HANDED +}; + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +void basis_matrix(mat4_t res, vec3_t i, vec3_t j, vec3_t k); +void basis_matrix_dir(mat4_t res, vec3_t dir); + +#ifdef __cplusplus +} /* extern "C" */ + +class Basis { +public: + Vector3 i, j, k; + + Basis(); + Basis(const Vector3 &i, const Vector3 &j, const Vector3 &k); + Basis(const Vector3 &dir, bool left_handed = true); + + void rotate(scalar_t x, scalar_t y, scalar_t z); + void rotate(const Vector3 &axis, scalar_t angle); + void rotate(const Matrix4x4 &mat); + void rotate(const Quaternion &quat); + + Matrix3x3 create_rotation_matrix() const; +}; +#endif /* __cplusplus */ + +#endif /* VMATH_BASIS_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/basis_c.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/basis_c.c Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,37 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include "basis.h" +#include "matrix.h" + +void basis_matrix(mat4_t res, vec3_t i, vec3_t j, vec3_t k) +{ + m4_identity(res); + m4_set_column(res, v4_cons(i.x, i.y, i.z, 1.0), 0); + m4_set_column(res, v4_cons(j.x, j.y, j.z, 1.0), 1); + m4_set_column(res, v4_cons(k.x, k.y, k.z, 1.0), 2); +} + +void basis_matrix_dir(mat4_t res, vec3_t dir) +{ + vec3_t k = v3_normalize(dir); + vec3_t j = {0, 1, 0}; + vec3_t i = v3_cross(j, k); + j = v3_cross(k, i); + basis_matrix(res, i, j, k); +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/geom.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/geom.c Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,150 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + + +#include +#include "geom.h" +#include "vector.h" + +plane_t plane_cons(scalar_t nx, scalar_t ny, scalar_t nz, scalar_t d) +{ + plane_t p; + p.norm.x = nx; + p.norm.y = ny; + p.norm.z = nz; + p.d = d; + return p; +} + +plane_t plane_poly(vec3_t v0, vec3_t v1, vec3_t v2) +{ + vec3_t a, b, norm; + + a = v3_sub(v1, v0); + b = v3_sub(v2, v0); + norm = v3_cross(a, b); + norm = v3_normalize(norm); + + return plane_ptnorm(v0, norm); +} + +plane_t plane_ptnorm(vec3_t pt, vec3_t normal) +{ + plane_t plane; + + plane.norm = normal; + plane.d = v3_dot(pt, normal); + + return plane; +} + +plane_t plane_invert(plane_t p) +{ + p.norm = v3_neg(p.norm); + p.d = -p.d; + return p; +} + +scalar_t plane_signed_dist(plane_t plane, vec3_t pt) +{ + vec3_t pp = plane_point(plane); + vec3_t pptopt = v3_sub(pt, pp); + return v3_dot(pptopt, plane.norm); +} + +scalar_t plane_dist(plane_t plane, vec3_t pt) +{ + return fabs(plane_signed_dist(plane, pt)); +} + +vec3_t plane_point(plane_t plane) +{ + return v3_scale(plane.norm, plane.d); +} + +int plane_ray_intersect(ray_t ray, plane_t plane, scalar_t *pos) +{ + vec3_t pt, orig_to_pt; + scalar_t ndotdir; + + pt = plane_point(plane); + ndotdir = v3_dot(plane.norm, ray.dir); + + if(fabs(ndotdir) < 1e-7) { + return 0; + } + + if(pos) { + orig_to_pt = v3_sub(pt, ray.origin); + *pos = v3_dot(plane.norm, orig_to_pt) / ndotdir; + } + return 1; +} + +sphere_t sphere_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t rad) +{ + sphere_t sph; + sph.pos.x = x; + sph.pos.y = y; + sph.pos.z = z; + sph.rad = rad; + return sph; +} + +int sphere_ray_intersect(ray_t ray, sphere_t sph, scalar_t *pos) +{ + scalar_t a, b, c, d, sqrt_d, t1, t2, t; + + a = v3_dot(ray.dir, ray.dir); + b = 2.0 * ray.dir.x * (ray.origin.x - sph.pos.x) + + 2.0 * ray.dir.y * (ray.origin.y - sph.pos.y) + + 2.0 * ray.dir.z * (ray.origin.z - sph.pos.z); + c = v3_dot(sph.pos, sph.pos) + v3_dot(ray.origin, ray.origin) + + 2.0 * v3_dot(v3_neg(sph.pos), ray.origin) - sph.rad * sph.rad; + + d = b * b - 4.0 * a * c; + if(d < 0.0) { + return 0; + } + + sqrt_d = sqrt(d); + t1 = (-b + sqrt_d) / (2.0 * a); + t2 = (-b - sqrt_d) / (2.0 * a); + + if(t1 < 1e-7 || t1 > 1.0) { + t1 = t2; + } + if(t2 < 1e-7 || t2 > 1.0) { + t2 = t1; + } + t = t1 < t2 ? t1 : t2; + + if(t < 1e-7 || t > 1.0) { + return 0; + } + + if(pos) { + *pos = t; + } + return 1; +} + +int sphere_sphere_intersect(sphere_t sph1, sphere_t sph2, scalar_t *pos, scalar_t *rad) +{ + return -1; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/geom.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/geom.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,72 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2012 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ +#ifndef LIBVMATH_GEOM_H_ +#define LIBVMATH_GEOM_H_ + +#include "vector.h" +#include "ray.h" + +typedef struct { + vec3_t norm; + scalar_t d; +} plane_t; + +typedef struct { + vec3_t pos; + scalar_t rad; +} sphere_t; + +typedef struct { + vec3_t min, max; +} aabox_t; + +#ifdef __cplusplus +extern "C" { +#endif + +/* planes are good... you need planes, yes you do */ +plane_t plane_cons(scalar_t nx, scalar_t ny, scalar_t nz, scalar_t d); +plane_t plane_poly(vec3_t v0, vec3_t v1, vec3_t v2); +plane_t plane_ptnorm(vec3_t pt, vec3_t normal); + +plane_t plane_invert(plane_t p); + +scalar_t plane_signed_dist(plane_t plane, vec3_t pt); +scalar_t plane_dist(plane_t plane, vec3_t pt); +vec3_t plane_point(plane_t plane); + +int plane_ray_intersect(ray_t ray, plane_t plane, scalar_t *pos); + +/* spheres always come in handy */ +sphere_t sphere_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t rad); + +int sphere_ray_intersect(ray_t ray, sphere_t sph, scalar_t *pos); +int sphere_sphere_intersect(sphere_t sph1, sphere_t sph2, scalar_t *pos, scalar_t *rad); + +#ifdef __cplusplus +} + +/* TODO +class Plane : public plane_t { +public: +}; +*/ + +#endif + +#endif /* LIBVMATH_GEOM_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/matrix.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/matrix.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,743 @@ +#include +#include +#include "matrix.h" +#include "vector.h" +#include "quat.h" + +using namespace std; + +// ----------- Matrix3x3 -------------- + +Matrix3x3 Matrix3x3::identity = Matrix3x3(1, 0, 0, 0, 1, 0, 0, 0, 1); + +Matrix3x3::Matrix3x3() +{ + *this = identity; +} + +Matrix3x3::Matrix3x3( scalar_t m11, scalar_t m12, scalar_t m13, + scalar_t m21, scalar_t m22, scalar_t m23, + scalar_t m31, scalar_t m32, scalar_t m33) +{ + m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; + m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; + m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; +} + +Matrix3x3::Matrix3x3(const Vector3 &ivec, const Vector3 &jvec, const Vector3 &kvec) +{ + set_row_vector(ivec, 0); + set_row_vector(jvec, 1); + set_row_vector(kvec, 2); +} + +Matrix3x3::Matrix3x3(const mat3_t cmat) +{ + memcpy(m, cmat, sizeof(mat3_t)); +} + +Matrix3x3::Matrix3x3(const Matrix4x4 &mat4x4) +{ + for(int i=0; i<3; i++) { + for(int j=0; j<3; j++) { + m[i][j] = mat4x4[i][j]; + } + } +} + +Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2) +{ + Matrix3x3 res; + const scalar_t *op1 = m1.m[0], *op2 = m2.m[0]; + scalar_t *dest = res.m[0]; + + for(int i=0; i<9; i++) { + *dest++ = *op1++ + *op2++; + } + return res; +} + +Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2) +{ + Matrix3x3 res; + const scalar_t *op1 = m1.m[0], *op2 = m2.m[0]; + scalar_t *dest = res.m[0]; + + for(int i=0; i<9; i++) { + *dest++ = *op1++ - *op2++; + } + return res; +} + +Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2) +{ + Matrix3x3 res; + for(int i=0; i<3; i++) { + for(int j=0; j<3; j++) { + res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j]; + } + } + return res; +} + +void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2) +{ + scalar_t *op1 = m1.m[0]; + const scalar_t *op2 = m2.m[0]; + + for(int i=0; i<9; i++) { + *op1++ += *op2++; + } +} + +void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2) +{ + scalar_t *op1 = m1.m[0]; + const scalar_t *op2 = m2.m[0]; + + for(int i=0; i<9; i++) { + *op1++ -= *op2++; + } +} + +void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2) +{ + Matrix3x3 res; + for(int i=0; i<3; i++) { + for(int j=0; j<3; j++) { + res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j]; + } + } + memcpy(m1.m, res.m, 9 * sizeof(scalar_t)); +} + +Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar) +{ + Matrix3x3 res; + const scalar_t *mptr = mat.m[0]; + scalar_t *dptr = res.m[0]; + + for(int i=0; i<9; i++) { + *dptr++ = *mptr++ * scalar; + } + return res; +} + +Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat) +{ + Matrix3x3 res; + const scalar_t *mptr = mat.m[0]; + scalar_t *dptr = res.m[0]; + + for(int i=0; i<9; i++) { + *dptr++ = *mptr++ * scalar; + } + return res; +} + +void operator *=(Matrix3x3 &mat, scalar_t scalar) +{ + scalar_t *mptr = mat.m[0]; + + for(int i=0; i<9; i++) { + *mptr++ *= scalar; + } +} + +void Matrix3x3::translate(const Vector2 &trans) +{ + Matrix3x3 tmat(1, 0, trans.x, 0, 1, trans.y, 0, 0, 1); + *this *= tmat; +} + +void Matrix3x3::set_translation(const Vector2 &trans) +{ + *this = Matrix3x3(1, 0, trans.x, 0, 1, trans.y, 0, 0, 1); +} + +void Matrix3x3::rotate(scalar_t angle) +{ + scalar_t cos_a = cos(angle); + scalar_t sin_a = sin(angle); + Matrix3x3 rmat( cos_a, -sin_a, 0, + sin_a, cos_a, 0, + 0, 0, 1); + *this *= rmat; +} + +void Matrix3x3::set_rotation(scalar_t angle) +{ + scalar_t cos_a = cos(angle); + scalar_t sin_a = sin(angle); + *this = Matrix3x3(cos_a, -sin_a, 0, sin_a, cos_a, 0, 0, 0, 1); +} + +void Matrix3x3::rotate(const Vector3 &euler_angles) +{ + Matrix3x3 xrot, yrot, zrot; + + xrot = Matrix3x3( 1, 0, 0, + 0, cos(euler_angles.x), -sin(euler_angles.x), + 0, sin(euler_angles.x), cos(euler_angles.x)); + + yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y), + 0, 1, 0, + -sin(euler_angles.y), 0, cos(euler_angles.y)); + + zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0, + sin(euler_angles.z), cos(euler_angles.z), 0, + 0, 0, 1); + + *this *= xrot * yrot * zrot; +} + +void Matrix3x3::set_rotation(const Vector3 &euler_angles) +{ + Matrix3x3 xrot, yrot, zrot; + + xrot = Matrix3x3( 1, 0, 0, + 0, cos(euler_angles.x), -sin(euler_angles.x), + 0, sin(euler_angles.x), cos(euler_angles.x)); + + yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y), + 0, 1, 0, + -sin(euler_angles.y), 0, cos(euler_angles.y)); + + zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0, + sin(euler_angles.z), cos(euler_angles.z), 0, + 0, 0, 1); + + *this = xrot * yrot * zrot; +} + +void Matrix3x3::rotate(const Vector3 &axis, scalar_t angle) +{ + scalar_t sina = (scalar_t)sin(angle); + scalar_t cosa = (scalar_t)cos(angle); + scalar_t invcosa = 1-cosa; + scalar_t nxsq = axis.x * axis.x; + scalar_t nysq = axis.y * axis.y; + scalar_t nzsq = axis.z * axis.z; + + Matrix3x3 xform; + xform.m[0][0] = nxsq + (1-nxsq) * cosa; + xform.m[0][1] = axis.x * axis.y * invcosa - axis.z * sina; + xform.m[0][2] = axis.x * axis.z * invcosa + axis.y * sina; + + xform.m[1][0] = axis.x * axis.y * invcosa + axis.z * sina; + xform.m[1][1] = nysq + (1-nysq) * cosa; + xform.m[1][2] = axis.y * axis.z * invcosa - axis.x * sina; + + xform.m[2][0] = axis.x * axis.z * invcosa - axis.y * sina; + xform.m[2][1] = axis.y * axis.z * invcosa + axis.x * sina; + xform.m[2][2] = nzsq + (1-nzsq) * cosa; + + *this *= xform; +} + +void Matrix3x3::set_rotation(const Vector3 &axis, scalar_t angle) +{ + scalar_t sina = (scalar_t)sin(angle); + scalar_t cosa = (scalar_t)cos(angle); + scalar_t invcosa = 1-cosa; + scalar_t nxsq = axis.x * axis.x; + scalar_t nysq = axis.y * axis.y; + scalar_t nzsq = axis.z * axis.z; + + reset_identity(); + m[0][0] = nxsq + (1-nxsq) * cosa; + m[0][1] = axis.x * axis.y * invcosa - axis.z * sina; + m[0][2] = axis.x * axis.z * invcosa + axis.y * sina; + m[1][0] = axis.x * axis.y * invcosa + axis.z * sina; + m[1][1] = nysq + (1-nysq) * cosa; + m[1][2] = axis.y * axis.z * invcosa - axis.x * sina; + m[2][0] = axis.x * axis.z * invcosa - axis.y * sina; + m[2][1] = axis.y * axis.z * invcosa + axis.x * sina; + m[2][2] = nzsq + (1-nzsq) * cosa; +} + +void Matrix3x3::scale(const Vector3 &scale_vec) +{ + Matrix3x3 smat( scale_vec.x, 0, 0, + 0, scale_vec.y, 0, + 0, 0, scale_vec.z); + *this *= smat; +} + +void Matrix3x3::set_scaling(const Vector3 &scale_vec) +{ + *this = Matrix3x3( scale_vec.x, 0, 0, + 0, scale_vec.y, 0, + 0, 0, scale_vec.z); +} + +void Matrix3x3::set_column_vector(const Vector3 &vec, unsigned int col_index) +{ + m[0][col_index] = vec.x; + m[1][col_index] = vec.y; + m[2][col_index] = vec.z; +} + +void Matrix3x3::set_row_vector(const Vector3 &vec, unsigned int row_index) +{ + m[row_index][0] = vec.x; + m[row_index][1] = vec.y; + m[row_index][2] = vec.z; +} + +Vector3 Matrix3x3::get_column_vector(unsigned int col_index) const +{ + return Vector3(m[0][col_index], m[1][col_index], m[2][col_index]); +} + +Vector3 Matrix3x3::get_row_vector(unsigned int row_index) const +{ + return Vector3(m[row_index][0], m[row_index][1], m[row_index][2]); +} + +void Matrix3x3::transpose() +{ + Matrix3x3 tmp = *this; + for(int i=0; i<3; i++) { + for(int j=0; j<3; j++) { + m[i][j] = tmp[j][i]; + } + } +} + +Matrix3x3 Matrix3x3::transposed() const +{ + Matrix3x3 res; + for(int i=0; i<3; i++) { + for(int j=0; j<3; j++) { + res[i][j] = m[j][i]; + } + } + return res; +} + +scalar_t Matrix3x3::determinant() const +{ + return m[0][0] * (m[1][1]*m[2][2] - m[1][2]*m[2][1]) - + m[0][1] * (m[1][0]*m[2][2] - m[1][2]*m[2][0]) + + m[0][2] * (m[1][0]*m[2][1] - m[1][1]*m[2][0]); +} + +Matrix3x3 Matrix3x3::inverse() const +{ + // TODO: implement 3x3 inverse + return *this; +} + +ostream &operator <<(ostream &out, const Matrix3x3 &mat) +{ + for(int i=0; i<3; i++) { + char str[100]; + sprintf(str, "[ %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2]); + out << str; + } + return out; +} + + + +/* ----------------- Matrix4x4 implementation --------------- */ + +Matrix4x4 Matrix4x4::identity = Matrix4x4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); + +Matrix4x4::Matrix4x4() +{ + *this = identity; +} + +Matrix4x4::Matrix4x4( scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14, + scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24, + scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34, + scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44) +{ + m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14; + m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24; + m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34; + m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44; +} + +Matrix4x4::Matrix4x4(const mat4_t cmat) +{ + memcpy(m, cmat, sizeof(mat4_t)); +} + +Matrix4x4::Matrix4x4(const Matrix3x3 &mat3x3) +{ + reset_identity(); + for(int i=0; i<3; i++) { + for(int j=0; j<3; j++) { + m[i][j] = mat3x3[i][j]; + } + } +} + +Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2) +{ + Matrix4x4 res; + const scalar_t *op1 = m1.m[0], *op2 = m2.m[0]; + scalar_t *dest = res.m[0]; + + for(int i=0; i<16; i++) { + *dest++ = *op1++ + *op2++; + } + return res; +} + +Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2) +{ + Matrix4x4 res; + const scalar_t *op1 = m1.m[0], *op2 = m2.m[0]; + scalar_t *dest = res.m[0]; + + for(int i=0; i<16; i++) { + *dest++ = *op1++ - *op2++; + } + return res; +} + +/* +Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2) { + Matrix4x4 res; + + for(int i=0; i<4; i++) { + for(int j=0; j<4; j++) { + res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j]; + } + } + + return res; +} +*/ + +void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2) +{ + scalar_t *op1 = m1.m[0]; + const scalar_t *op2 = m2.m[0]; + + for(int i=0; i<16; i++) { + *op1++ += *op2++; + } +} + +void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2) +{ + scalar_t *op1 = m1.m[0]; + const scalar_t *op2 = m2.m[0]; + + for(int i=0; i<16; i++) { + *op1++ -= *op2++; + } +} + +void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2) +{ + Matrix4x4 res; + for(int i=0; i<4; i++) { + for(int j=0; j<4; j++) { + res.m[i][j] = m1.m[i][0] * m2.m[0][j] + m1.m[i][1] * m2.m[1][j] + m1.m[i][2] * m2.m[2][j] + m1.m[i][3] * m2.m[3][j]; + } + } + memcpy(m1.m, res.m, 16 * sizeof(scalar_t)); +} + +Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar) +{ + Matrix4x4 res; + const scalar_t *mptr = mat.m[0]; + scalar_t *dptr = res.m[0]; + + for(int i=0; i<16; i++) { + *dptr++ = *mptr++ * scalar; + } + return res; +} + +Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat) +{ + Matrix4x4 res; + const scalar_t *mptr = mat.m[0]; + scalar_t *dptr = res.m[0]; + + for(int i=0; i<16; i++) { + *dptr++ = *mptr++ * scalar; + } + return res; +} + +void operator *=(Matrix4x4 &mat, scalar_t scalar) +{ + scalar_t *mptr = mat.m[0]; + + for(int i=0; i<16; i++) { + *mptr++ *= scalar; + } +} + +void Matrix4x4::translate(const Vector3 &trans) +{ + Matrix4x4 tmat(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1); + *this *= tmat; +} + +void Matrix4x4::set_translation(const Vector3 &trans) +{ + *this = Matrix4x4(1, 0, 0, trans.x, 0, 1, 0, trans.y, 0, 0, 1, trans.z, 0, 0, 0, 1); +} + +void Matrix4x4::rotate(const Vector3 &euler_angles) +{ + Matrix3x3 xrot, yrot, zrot; + + xrot = Matrix3x3( 1, 0, 0, + 0, cos(euler_angles.x), -sin(euler_angles.x), + 0, sin(euler_angles.x), cos(euler_angles.x)); + + yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y), + 0, 1, 0, + -sin(euler_angles.y), 0, cos(euler_angles.y)); + + zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0, + sin(euler_angles.z), cos(euler_angles.z), 0, + 0, 0, 1); + + *this *= Matrix4x4(xrot * yrot * zrot); +} + +void Matrix4x4::set_rotation(const Vector3 &euler_angles) +{ + Matrix3x3 xrot, yrot, zrot; + + xrot = Matrix3x3( 1, 0, 0, + 0, cos(euler_angles.x), -sin(euler_angles.x), + 0, sin(euler_angles.x), cos(euler_angles.x)); + + yrot = Matrix3x3( cos(euler_angles.y), 0, sin(euler_angles.y), + 0, 1, 0, + -sin(euler_angles.y), 0, cos(euler_angles.y)); + + zrot = Matrix3x3( cos(euler_angles.z), -sin(euler_angles.z), 0, + sin(euler_angles.z), cos(euler_angles.z), 0, + 0, 0, 1); + + *this = Matrix4x4(xrot * yrot * zrot); +} + +void Matrix4x4::rotate(const Vector3 &axis, scalar_t angle) +{ + scalar_t sina = (scalar_t)sin(angle); + scalar_t cosa = (scalar_t)cos(angle); + scalar_t invcosa = 1-cosa; + scalar_t nxsq = axis.x * axis.x; + scalar_t nysq = axis.y * axis.y; + scalar_t nzsq = axis.z * axis.z; + + Matrix3x3 xform; + xform[0][0] = nxsq + (1-nxsq) * cosa; + xform[0][1] = axis.x * axis.y * invcosa - axis.z * sina; + xform[0][2] = axis.x * axis.z * invcosa + axis.y * sina; + xform[1][0] = axis.x * axis.y * invcosa + axis.z * sina; + xform[1][1] = nysq + (1-nysq) * cosa; + xform[1][2] = axis.y * axis.z * invcosa - axis.x * sina; + xform[2][0] = axis.x * axis.z * invcosa - axis.y * sina; + xform[2][1] = axis.y * axis.z * invcosa + axis.x * sina; + xform[2][2] = nzsq + (1-nzsq) * cosa; + + *this *= Matrix4x4(xform); +} + +void Matrix4x4::set_rotation(const Vector3 &axis, scalar_t angle) +{ + scalar_t sina = (scalar_t)sin(angle); + scalar_t cosa = (scalar_t)cos(angle); + scalar_t invcosa = 1-cosa; + scalar_t nxsq = axis.x * axis.x; + scalar_t nysq = axis.y * axis.y; + scalar_t nzsq = axis.z * axis.z; + + reset_identity(); + m[0][0] = nxsq + (1-nxsq) * cosa; + m[0][1] = axis.x * axis.y * invcosa - axis.z * sina; + m[0][2] = axis.x * axis.z * invcosa + axis.y * sina; + m[1][0] = axis.x * axis.y * invcosa + axis.z * sina; + m[1][1] = nysq + (1-nysq) * cosa; + m[1][2] = axis.y * axis.z * invcosa - axis.x * sina; + m[2][0] = axis.x * axis.z * invcosa - axis.y * sina; + m[2][1] = axis.y * axis.z * invcosa + axis.x * sina; + m[2][2] = nzsq + (1-nzsq) * cosa; +} + +void Matrix4x4::scale(const Vector4 &scale_vec) +{ + Matrix4x4 smat( scale_vec.x, 0, 0, 0, + 0, scale_vec.y, 0, 0, + 0, 0, scale_vec.z, 0, + 0, 0, 0, scale_vec.w); + *this *= smat; +} + +void Matrix4x4::set_scaling(const Vector4 &scale_vec) +{ + *this = Matrix4x4( scale_vec.x, 0, 0, 0, + 0, scale_vec.y, 0, 0, + 0, 0, scale_vec.z, 0, + 0, 0, 0, scale_vec.w); +} + +void Matrix4x4::set_column_vector(const Vector4 &vec, unsigned int col_index) +{ + m[0][col_index] = vec.x; + m[1][col_index] = vec.y; + m[2][col_index] = vec.z; + m[3][col_index] = vec.w; +} + +void Matrix4x4::set_row_vector(const Vector4 &vec, unsigned int row_index) +{ + m[row_index][0] = vec.x; + m[row_index][1] = vec.y; + m[row_index][2] = vec.z; + m[row_index][3] = vec.w; +} + +Vector4 Matrix4x4::get_column_vector(unsigned int col_index) const +{ + return Vector4(m[0][col_index], m[1][col_index], m[2][col_index], m[3][col_index]); +} + +Vector4 Matrix4x4::get_row_vector(unsigned int row_index) const +{ + return Vector4(m[row_index][0], m[row_index][1], m[row_index][2], m[row_index][3]); +} + +void Matrix4x4::transpose() +{ + Matrix4x4 tmp = *this; + for(int i=0; i<4; i++) { + for(int j=0; j<4; j++) { + m[i][j] = tmp[j][i]; + } + } +} + +Matrix4x4 Matrix4x4::transposed() const +{ + Matrix4x4 res; + for(int i=0; i<4; i++) { + for(int j=0; j<4; j++) { + res[i][j] = m[j][i]; + } + } + return res; +} + +scalar_t Matrix4x4::determinant() const +{ + scalar_t det11 = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + + (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); + + scalar_t det12 = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); + + scalar_t det13 = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - + (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + scalar_t det14 = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - + (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + + (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + return m[0][0] * det11 - m[0][1] * det12 + m[0][2] * det13 - m[0][3] * det14; +} + + +Matrix4x4 Matrix4x4::adjoint() const +{ + Matrix4x4 coef; + + coef.m[0][0] = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + + (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); + coef.m[0][1] = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); + coef.m[0][2] = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - + (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + coef.m[0][3] = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - + (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + + (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + coef.m[1][0] = (m[0][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[0][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + + (m[0][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); + coef.m[1][1] = (m[0][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[0][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[0][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); + coef.m[1][2] = (m[0][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - + (m[0][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[0][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + coef.m[1][3] = (m[0][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - + (m[0][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + + (m[0][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + coef.m[2][0] = (m[0][1] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) - + (m[0][2] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) + + (m[0][3] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])); + coef.m[2][1] = (m[0][0] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) - + (m[0][2] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])); + coef.m[2][2] = (m[0][0] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) - + (m[0][1] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[3][1] - m[3][0] * m[1][1])); + coef.m[2][3] = (m[0][0] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])) - + (m[0][1] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])) + + (m[0][2] * (m[1][0] * m[3][1] - m[3][0] * m[1][1])); + + coef.m[3][0] = (m[0][1] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) - + (m[0][2] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) + + (m[0][3] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])); + coef.m[3][1] = (m[0][0] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) - + (m[0][2] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])); + coef.m[3][2] = (m[0][0] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) - + (m[0][1] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[2][1] - m[2][0] * m[1][1])); + coef.m[3][3] = (m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])) - + (m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])) + + (m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1])); + + coef.transpose(); + + for(int i=0; i<4; i++) { + for(int j=0; j<4; j++) { + coef.m[i][j] = j%2 ? -coef.m[i][j] : coef.m[i][j]; + if(i%2) coef.m[i][j] = -coef.m[i][j]; + } + } + + return coef; +} + +Matrix4x4 Matrix4x4::inverse() const +{ + Matrix4x4 adj = adjoint(); + + return adj * (1.0f / determinant()); +} + +ostream &operator <<(ostream &out, const Matrix4x4 &mat) +{ + for(int i=0; i<4; i++) { + char str[100]; + sprintf(str, "[ %12.5f %12.5f %12.5f %12.5f ]\n", (float)mat.m[i][0], (float)mat.m[i][1], (float)mat.m[i][2], (float)mat.m[i][3]); + out << str; + } + return out; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/matrix.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/matrix.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,206 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_MATRIX_H_ +#define VMATH_MATRIX_H_ + +#include +#include "vmath_types.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +/* C matrix 3x3 functions */ +static inline void m3_identity(mat3_t m); +static inline void m3_cons(mat3_t m, + scalar_t m11, scalar_t m12, scalar_t m13, + scalar_t m21, scalar_t m22, scalar_t m23, + scalar_t m31, scalar_t m32, scalar_t m33); +static inline void m3_copy(mat3_t dest, mat3_t src); +void m3_to_m4(mat4_t dest, mat3_t src); + +void m3_print(FILE *fp, mat3_t m); + +/* C matrix 4x4 functions */ +static inline void m4_identity(mat4_t m); +static inline void m4_cons(mat4_t m, + scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14, + scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24, + scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34, + scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44); +static inline void m4_copy(mat4_t dest, mat4_t src); +void m4_to_m3(mat3_t dest, mat4_t src); + +static inline void m4_mult(mat4_t res, mat4_t m1, mat4_t m2); + +void m4_translate(mat4_t m, scalar_t x, scalar_t y, scalar_t z); +void m4_rotate(mat4_t m, scalar_t x, scalar_t y, scalar_t z); +void m4_rotate_x(mat4_t m, scalar_t angle); +void m4_rotate_y(mat4_t m, scalar_t angle); +void m4_rotate_z(mat4_t m, scalar_t angle); +void m4_rotate_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z); +void m4_rotate_quat(mat4_t m, quat_t q); +void m4_scale(mat4_t m, scalar_t x, scalar_t y, scalar_t z); +static inline void m4_set_column(mat4_t m, vec4_t v, int idx); +static inline void m4_set_row(mat4_t m, vec4_t v, int idx); + +void m4_transpose(mat4_t res, mat4_t m); +scalar_t m4_determinant(mat4_t m); +void m4_adjoint(mat4_t res, mat4_t m); +void m4_inverse(mat4_t res, mat4_t m); + +void m4_print(FILE *fp, mat4_t m); + +#ifdef __cplusplus +} + +/* when included from C++ source files, also define the matrix classes */ +#include + +/** 3x3 matrix */ +class Matrix3x3 { +private: + scalar_t m[3][3]; + +public: + + static Matrix3x3 identity; + + Matrix3x3(); + Matrix3x3( scalar_t m11, scalar_t m12, scalar_t m13, + scalar_t m21, scalar_t m22, scalar_t m23, + scalar_t m31, scalar_t m32, scalar_t m33); + Matrix3x3(const Vector3 &ivec, const Vector3 &jvec, const Vector3 &kvec); + Matrix3x3(const mat3_t cmat); + + Matrix3x3(const Matrix4x4 &mat4x4); + + /* binary operations matrix (op) matrix */ + friend Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2); + friend Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2); + friend Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2); + + friend void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2); + friend void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2); + friend void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2); + + /* binary operations matrix (op) scalar and scalar (op) matrix */ + friend Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar); + friend Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat); + + friend void operator *=(Matrix3x3 &mat, scalar_t scalar); + + inline scalar_t *operator [](int index); + inline const scalar_t *operator [](int index) const; + + inline void reset_identity(); + + void translate(const Vector2 &trans); + void set_translation(const Vector2 &trans); + + void rotate(scalar_t angle); /* 2d rotation */ + void rotate(const Vector3 &euler_angles); /* 3d rotation with euler angles */ + void rotate(const Vector3 &axis, scalar_t angle); /* 3d axis/angle rotation */ + void set_rotation(scalar_t angle); + void set_rotation(const Vector3 &euler_angles); + void set_rotation(const Vector3 &axis, scalar_t angle); + + void scale(const Vector3 &scale_vec); + void set_scaling(const Vector3 &scale_vec); + + void set_column_vector(const Vector3 &vec, unsigned int col_index); + void set_row_vector(const Vector3 &vec, unsigned int row_index); + Vector3 get_column_vector(unsigned int col_index) const; + Vector3 get_row_vector(unsigned int row_index) const; + + void transpose(); + Matrix3x3 transposed() const; + scalar_t determinant() const; + Matrix3x3 inverse() const; + + friend std::ostream &operator <<(std::ostream &out, const Matrix3x3 &mat); +}; + +/** 4x4 matrix */ +class Matrix4x4 { +private: + scalar_t m[4][4]; + +public: + + static Matrix4x4 identity; + + Matrix4x4(); + Matrix4x4( scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14, + scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24, + scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34, + scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44); + Matrix4x4(const mat4_t cmat); + + Matrix4x4(const Matrix3x3 &mat3x3); + + /* binary operations matrix (op) matrix */ + friend Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2); + friend Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2); + friend Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2); + + friend void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2); + friend void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2); + friend inline void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2); + + /* binary operations matrix (op) scalar and scalar (op) matrix */ + friend Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar); + friend Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat); + + friend void operator *=(Matrix4x4 &mat, scalar_t scalar); + + inline scalar_t *operator [](int index); + inline const scalar_t *operator [](int index) const; + + inline void reset_identity(); + + void translate(const Vector3 &trans); + void set_translation(const Vector3 &trans); + + void rotate(const Vector3 &euler_angles); /* 3d rotation with euler angles */ + void rotate(const Vector3 &axis, scalar_t angle); /* 3d axis/angle rotation */ + void set_rotation(const Vector3 &euler_angles); + void set_rotation(const Vector3 &axis, scalar_t angle); + + void scale(const Vector4 &scale_vec); + void set_scaling(const Vector4 &scale_vec); + + void set_column_vector(const Vector4 &vec, unsigned int col_index); + void set_row_vector(const Vector4 &vec, unsigned int row_index); + Vector4 get_column_vector(unsigned int col_index) const; + Vector4 get_row_vector(unsigned int row_index) const; + + void transpose(); + Matrix4x4 transposed() const; + scalar_t determinant() const; + Matrix4x4 adjoint() const; + Matrix4x4 inverse() const; + + friend std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat); +}; +#endif /* __cplusplus */ + +#include "matrix.inl" + +#endif /* VMATH_MATRIX_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/matrix.inl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/matrix.inl Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,165 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +/* C matrix 3x3 functions */ +static inline void m3_identity(mat3_t m) +{ + static const mat3_t id = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}; + memcpy(m, id, sizeof id); +} + +static inline void m3_cons(mat3_t m, + scalar_t m11, scalar_t m12, scalar_t m13, + scalar_t m21, scalar_t m22, scalar_t m23, + scalar_t m31, scalar_t m32, scalar_t m33) +{ + m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; + m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; + m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; +} + +static inline void m3_copy(mat3_t dest, mat3_t src) +{ + memcpy(dest, src, sizeof(mat3_t)); +} + + +/* C matrix 4x4 functions */ +static inline void m4_identity(mat4_t m) +{ + static const mat4_t id = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}; + memcpy(m, id, sizeof id); +} + +static inline void m4_cons(mat4_t m, + scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14, + scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24, + scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34, + scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44) +{ + m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[0][3] = m14; + m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[1][3] = m24; + m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; m[2][3] = m34; + m[3][0] = m41; m[3][1] = m42; m[3][2] = m43; m[3][3] = m44; +} + +static inline void m4_copy(mat4_t dest, mat4_t src) +{ + memcpy(dest, src, sizeof(mat4_t)); +} + +static inline void m4_mult(mat4_t res, mat4_t m1, mat4_t m2) +{ + res[0][0] = m1[0][0] * m2[0][0] + m1[0][1] * m2[1][0] + m1[0][2] * m2[2][0] + m1[0][3] * m2[3][0]; + res[0][1] = m1[0][0] * m2[0][1] + m1[0][1] * m2[1][1] + m1[0][2] * m2[2][1] + m1[0][3] * m2[3][1]; + res[0][2] = m1[0][0] * m2[0][2] + m1[0][1] * m2[1][2] + m1[0][2] * m2[2][2] + m1[0][3] * m2[3][2]; + res[0][3] = m1[0][0] * m2[0][3] + m1[0][1] * m2[1][3] + m1[0][2] * m2[2][3] + m1[0][3] * m2[3][3]; + + res[1][0] = m1[1][0] * m2[0][0] + m1[1][1] * m2[1][0] + m1[1][2] * m2[2][0] + m1[1][3] * m2[3][0]; + res[1][1] = m1[1][0] * m2[0][1] + m1[1][1] * m2[1][1] + m1[1][2] * m2[2][1] + m1[1][3] * m2[3][1]; + res[1][2] = m1[1][0] * m2[0][2] + m1[1][1] * m2[1][2] + m1[1][2] * m2[2][2] + m1[1][3] * m2[3][2]; + res[1][3] = m1[1][0] * m2[0][3] + m1[1][1] * m2[1][3] + m1[1][2] * m2[2][3] + m1[1][3] * m2[3][3]; + + res[2][0] = m1[2][0] * m2[0][0] + m1[2][1] * m2[1][0] + m1[2][2] * m2[2][0] + m1[2][3] * m2[3][0]; + res[2][1] = m1[2][0] * m2[0][1] + m1[2][1] * m2[1][1] + m1[2][2] * m2[2][1] + m1[2][3] * m2[3][1]; + res[2][2] = m1[2][0] * m2[0][2] + m1[2][1] * m2[1][2] + m1[2][2] * m2[2][2] + m1[2][3] * m2[3][2]; + res[2][3] = m1[2][0] * m2[0][3] + m1[2][1] * m2[1][3] + m1[2][2] * m2[2][3] + m1[2][3] * m2[3][3]; + + res[3][0] = m1[3][0] * m2[0][0] + m1[3][1] * m2[1][0] + m1[3][2] * m2[2][0] + m1[3][3] * m2[3][0]; + res[3][1] = m1[3][0] * m2[0][1] + m1[3][1] * m2[1][1] + m1[3][2] * m2[2][1] + m1[3][3] * m2[3][1]; + res[3][2] = m1[3][0] * m2[0][2] + m1[3][1] * m2[1][2] + m1[3][2] * m2[2][2] + m1[3][3] * m2[3][2]; + res[3][3] = m1[3][0] * m2[0][3] + m1[3][1] * m2[1][3] + m1[3][2] * m2[2][3] + m1[3][3] * m2[3][3]; +} + +static inline void m4_set_column(mat4_t m, vec4_t v, int idx) +{ + m[0][idx] = v.x; + m[1][idx] = v.y; + m[2][idx] = v.z; + m[3][idx] = v.w; +} + +static inline void m4_set_row(mat4_t m, vec4_t v, int idx) +{ + m[idx][0] = v.x; + m[idx][1] = v.y; + m[idx][2] = v.z; + m[idx][3] = v.w; +} + +#ifdef __cplusplus +} /* extern "C" */ + + +/* unrolled to hell and inline */ +inline Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2) { + Matrix4x4 res; + + res.m[0][0] = m1.m[0][0] * m2.m[0][0] + m1.m[0][1] * m2.m[1][0] + m1.m[0][2] * m2.m[2][0] + m1.m[0][3] * m2.m[3][0]; + res.m[0][1] = m1.m[0][0] * m2.m[0][1] + m1.m[0][1] * m2.m[1][1] + m1.m[0][2] * m2.m[2][1] + m1.m[0][3] * m2.m[3][1]; + res.m[0][2] = m1.m[0][0] * m2.m[0][2] + m1.m[0][1] * m2.m[1][2] + m1.m[0][2] * m2.m[2][2] + m1.m[0][3] * m2.m[3][2]; + res.m[0][3] = m1.m[0][0] * m2.m[0][3] + m1.m[0][1] * m2.m[1][3] + m1.m[0][2] * m2.m[2][3] + m1.m[0][3] * m2.m[3][3]; + + res.m[1][0] = m1.m[1][0] * m2.m[0][0] + m1.m[1][1] * m2.m[1][0] + m1.m[1][2] * m2.m[2][0] + m1.m[1][3] * m2.m[3][0]; + res.m[1][1] = m1.m[1][0] * m2.m[0][1] + m1.m[1][1] * m2.m[1][1] + m1.m[1][2] * m2.m[2][1] + m1.m[1][3] * m2.m[3][1]; + res.m[1][2] = m1.m[1][0] * m2.m[0][2] + m1.m[1][1] * m2.m[1][2] + m1.m[1][2] * m2.m[2][2] + m1.m[1][3] * m2.m[3][2]; + res.m[1][3] = m1.m[1][0] * m2.m[0][3] + m1.m[1][1] * m2.m[1][3] + m1.m[1][2] * m2.m[2][3] + m1.m[1][3] * m2.m[3][3]; + + res.m[2][0] = m1.m[2][0] * m2.m[0][0] + m1.m[2][1] * m2.m[1][0] + m1.m[2][2] * m2.m[2][0] + m1.m[2][3] * m2.m[3][0]; + res.m[2][1] = m1.m[2][0] * m2.m[0][1] + m1.m[2][1] * m2.m[1][1] + m1.m[2][2] * m2.m[2][1] + m1.m[2][3] * m2.m[3][1]; + res.m[2][2] = m1.m[2][0] * m2.m[0][2] + m1.m[2][1] * m2.m[1][2] + m1.m[2][2] * m2.m[2][2] + m1.m[2][3] * m2.m[3][2]; + res.m[2][3] = m1.m[2][0] * m2.m[0][3] + m1.m[2][1] * m2.m[1][3] + m1.m[2][2] * m2.m[2][3] + m1.m[2][3] * m2.m[3][3]; + + res.m[3][0] = m1.m[3][0] * m2.m[0][0] + m1.m[3][1] * m2.m[1][0] + m1.m[3][2] * m2.m[2][0] + m1.m[3][3] * m2.m[3][0]; + res.m[3][1] = m1.m[3][0] * m2.m[0][1] + m1.m[3][1] * m2.m[1][1] + m1.m[3][2] * m2.m[2][1] + m1.m[3][3] * m2.m[3][1]; + res.m[3][2] = m1.m[3][0] * m2.m[0][2] + m1.m[3][1] * m2.m[1][2] + m1.m[3][2] * m2.m[2][2] + m1.m[3][3] * m2.m[3][2]; + res.m[3][3] = m1.m[3][0] * m2.m[0][3] + m1.m[3][1] * m2.m[1][3] + m1.m[3][2] * m2.m[2][3] + m1.m[3][3] * m2.m[3][3]; + + return res; +} + +inline scalar_t *Matrix3x3::operator [](int index) { + return m[index]; +} + +inline const scalar_t *Matrix3x3::operator [](int index) const { + return m[index]; +} + +inline void Matrix3x3::reset_identity() { + memcpy(this->m, identity.m, 9 * sizeof(scalar_t)); +} + +inline scalar_t *Matrix4x4::operator [](int index) { + return m[index]; +} + +inline const scalar_t *Matrix4x4::operator [](int index) const { + return m[index]; +} + +inline void Matrix4x4::reset_identity() { + memcpy(this->m, identity.m, 16 * sizeof(scalar_t)); +} +#endif /* __cplusplus */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/matrix_c.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/matrix_c.c Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,265 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + + +#include +#include "matrix.h" +#include "vector.h" +#include "quat.h" + +void m3_to_m4(mat4_t dest, mat3_t src) +{ + int i, j; + + memset(dest, 0, sizeof(mat4_t)); + for(i=0; i<3; i++) { + for(j=0; j<3; j++) { + dest[i][j] = src[i][j]; + } + } + dest[3][3] = 1.0; +} + +void m3_print(FILE *fp, mat3_t m) +{ + int i; + for(i=0; i<3; i++) { + fprintf(fp, "[ %12.5f %12.5f %12.5f ]\n", (float)m[i][0], (float)m[i][1], (float)m[i][2]); + } +} + +/* C matrix 4x4 functions */ +void m4_to_m3(mat3_t dest, mat4_t src) +{ + int i, j; + for(i=0; i<3; i++) { + for(j=0; j<3; j++) { + dest[i][j] = src[i][j]; + } + } +} + +void m4_translate(mat4_t m, scalar_t x, scalar_t y, scalar_t z) +{ + mat4_t tm; + m4_identity(tm); + tm[0][3] = x; + tm[1][3] = y; + tm[2][3] = z; + m4_mult(m, m, tm); +} + +void m4_rotate(mat4_t m, scalar_t x, scalar_t y, scalar_t z) +{ + m4_rotate_x(m, x); + m4_rotate_y(m, y); + m4_rotate_z(m, z); +} + +void m4_rotate_x(mat4_t m, scalar_t angle) +{ + mat4_t rm; + m4_identity(rm); + rm[1][1] = cos(angle); rm[1][2] = -sin(angle); + rm[2][1] = sin(angle); rm[2][2] = cos(angle); + m4_mult(m, m, rm); +} + +void m4_rotate_y(mat4_t m, scalar_t angle) +{ + mat4_t rm; + m4_identity(rm); + rm[0][0] = cos(angle); rm[0][2] = sin(angle); + rm[2][0] = -sin(angle); rm[2][2] = cos(angle); + m4_mult(m, m, rm); +} + +void m4_rotate_z(mat4_t m, scalar_t angle) +{ + mat4_t rm; + m4_identity(rm); + rm[0][0] = cos(angle); rm[0][1] = -sin(angle); + rm[1][0] = sin(angle); rm[1][1] = cos(angle); + m4_mult(m, m, rm); +} + +void m4_rotate_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z) +{ + mat4_t xform; + scalar_t sina = sin(angle); + scalar_t cosa = cos(angle); + scalar_t one_minus_cosa = 1.0 - cosa; + scalar_t nxsq = x * x; + scalar_t nysq = y * y; + scalar_t nzsq = z * z; + + m4_identity(xform); + xform[0][0] = nxsq + (1.0 - nxsq) * cosa; + xform[0][1] = x * y * one_minus_cosa - z * sina; + xform[0][2] = x * z * one_minus_cosa + y * sina; + xform[1][0] = x * y * one_minus_cosa + z * sina; + xform[1][1] = nysq + (1.0 - nysq) * cosa; + xform[1][2] = y * z * one_minus_cosa - x * sina; + xform[2][0] = x * z * one_minus_cosa - y * sina; + xform[2][1] = y * z * one_minus_cosa + x * sina; + xform[2][2] = nzsq + (1.0 - nzsq) * cosa; + + m4_mult(m, m, xform); +} + +void m4_rotate_quat(mat4_t m, quat_t q) +{ + mat4_t rm; + quat_to_mat4(rm, q); + m4_mult(m, m, rm); +} + +void m4_scale(mat4_t m, scalar_t x, scalar_t y, scalar_t z) +{ + mat4_t sm; + m4_identity(sm); + sm[0][0] = x; + sm[1][1] = y; + sm[2][2] = z; + m4_mult(m, m, sm); +} + +void m4_transpose(mat4_t res, mat4_t m) +{ + int i, j; + mat4_t tmp; + m4_copy(tmp, m); + + for(i=0; i<4; i++) { + for(j=0; j<4; j++) { + res[i][j] = tmp[j][i]; + } + } +} + +scalar_t m4_determinant(mat4_t m) +{ + scalar_t det11 = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + + (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); + + scalar_t det12 = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); + + scalar_t det13 = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - + (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + scalar_t det14 = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - + (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + + (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + return m[0][0] * det11 - m[0][1] * det12 + m[0][2] * det13 - m[0][3] * det14; +} + +void m4_adjoint(mat4_t res, mat4_t m) +{ + int i, j; + mat4_t coef; + + coef[0][0] = (m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + + (m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); + coef[0][1] = (m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); + coef[0][2] = (m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - + (m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + coef[0][3] = (m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - + (m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + + (m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + coef[1][0] = (m[0][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[0][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) + + (m[0][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])); + coef[1][1] = (m[0][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) - + (m[0][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[0][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])); + coef[1][2] = (m[0][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) - + (m[0][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) + + (m[0][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + coef[1][3] = (m[0][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) - + (m[0][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) + + (m[0][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1])); + + coef[2][0] = (m[0][1] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) - + (m[0][2] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) + + (m[0][3] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])); + coef[2][1] = (m[0][0] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) - + (m[0][2] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])); + coef[2][2] = (m[0][0] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) - + (m[0][1] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[3][1] - m[3][0] * m[1][1])); + coef[2][3] = (m[0][0] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])) - + (m[0][1] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])) + + (m[0][2] * (m[1][0] * m[3][1] - m[3][0] * m[1][1])); + + coef[3][0] = (m[0][1] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) - + (m[0][2] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) + + (m[0][3] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])); + coef[3][1] = (m[0][0] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) - + (m[0][2] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])); + coef[3][2] = (m[0][0] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) - + (m[0][1] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) + + (m[0][3] * (m[1][0] * m[2][1] - m[2][0] * m[1][1])); + coef[3][3] = (m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])) - + (m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])) + + (m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1])); + + m4_transpose(res, coef); + + for(i=0; i<4; i++) { + for(j=0; j<4; j++) { + res[i][j] = j % 2 ? -res[i][j] : res[i][j]; + if(i % 2) res[i][j] = -res[i][j]; + } + } +} + +void m4_inverse(mat4_t res, mat4_t m) +{ + int i, j; + mat4_t adj; + scalar_t det; + + m4_adjoint(adj, m); + det = m4_determinant(m); + + for(i=0; i<4; i++) { + for(j=0; j<4; j++) { + res[i][j] = adj[i][j] / det; + } + } +} + +void m4_print(FILE *fp, mat4_t m) +{ + int i; + for(i=0; i<4; i++) { + fprintf(fp, "[ %12.5f %12.5f %12.5f %12.5f ]\n", (float)m[i][0], (float)m[i][1], (float)m[i][2], (float)m[i][3]); + } +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/quat.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/quat.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,159 @@ +#include "quat.h" + +Quaternion::Quaternion() { + s = 1.0; + v.x = v.y = v.z = 0.0; +} + +Quaternion::Quaternion(scalar_t s, const Vector3 &v) { + this->s = s; + this->v = v; +} + +Quaternion::Quaternion(scalar_t s, scalar_t x, scalar_t y, scalar_t z) { + v.x = x; + v.y = y; + v.z = z; + this->s = s; +} + +Quaternion::Quaternion(const Vector3 &axis, scalar_t angle) { + set_rotation(axis, angle); +} + +Quaternion::Quaternion(const quat_t &quat) +{ + v.x = quat.x; + v.y = quat.y; + v.z = quat.z; + s = quat.w; +} + +Quaternion Quaternion::operator +(const Quaternion &quat) const { + return Quaternion(s + quat.s, v + quat.v); +} + +Quaternion Quaternion::operator -(const Quaternion &quat) const { + return Quaternion(s - quat.s, v - quat.v); +} + +Quaternion Quaternion::operator -() const { + return Quaternion(-s, -v); +} + +/** Quaternion Multiplication: + * Q1*Q2 = [s1*s2 - v1.v2, s1*v2 + s2*v1 + v1(x)v2] + */ +Quaternion Quaternion::operator *(const Quaternion &quat) const { + Quaternion newq; + newq.s = s * quat.s - dot_product(v, quat.v); + newq.v = quat.v * s + v * quat.s + cross_product(v, quat.v); + return newq; +} + +void Quaternion::operator +=(const Quaternion &quat) { + *this = Quaternion(s + quat.s, v + quat.v); +} + +void Quaternion::operator -=(const Quaternion &quat) { + *this = Quaternion(s - quat.s, v - quat.v); +} + +void Quaternion::operator *=(const Quaternion &quat) { + *this = *this * quat; +} + +void Quaternion::reset_identity() { + s = 1.0; + v.x = v.y = v.z = 0.0; +} + +Quaternion Quaternion::conjugate() const { + return Quaternion(s, -v); +} + +scalar_t Quaternion::length() const { + return (scalar_t)sqrt(v.x*v.x + v.y*v.y + v.z*v.z + s*s); +} + +/** Q * ~Q = ||Q||^2 */ +scalar_t Quaternion::length_sq() const { + return v.x*v.x + v.y*v.y + v.z*v.z + s*s; +} + +void Quaternion::normalize() { + scalar_t len = (scalar_t)sqrt(v.x*v.x + v.y*v.y + v.z*v.z + s*s); + v.x /= len; + v.y /= len; + v.z /= len; + s /= len; +} + +Quaternion Quaternion::normalized() const { + Quaternion nq = *this; + scalar_t len = (scalar_t)sqrt(v.x*v.x + v.y*v.y + v.z*v.z + s*s); + nq.v.x /= len; + nq.v.y /= len; + nq.v.z /= len; + nq.s /= len; + return nq; +} + +/** Quaternion Inversion: Q^-1 = ~Q / ||Q||^2 */ +Quaternion Quaternion::inverse() const { + Quaternion inv = conjugate(); + scalar_t lensq = length_sq(); + inv.v /= lensq; + inv.s /= lensq; + + return inv; +} + + +void Quaternion::set_rotation(const Vector3 &axis, scalar_t angle) { + scalar_t half_angle = angle / 2.0; + s = cos(half_angle); + v = axis * sin(half_angle); +} + +void Quaternion::rotate(const Vector3 &axis, scalar_t angle) { + Quaternion q; + scalar_t half_angle = angle / 2.0; + q.s = cos(half_angle); + q.v = axis * sin(half_angle); + + *this *= q; +} + +void Quaternion::rotate(const Quaternion &q) { + *this = q * *this * q.conjugate(); +} + +Matrix3x3 Quaternion::get_rotation_matrix() const { + return Matrix3x3( 1.0 - 2.0 * v.y*v.y - 2.0 * v.z*v.z, 2.0 * v.x * v.y + 2.0 * s * v.z, 2.0 * v.z * v.x - 2.0 * s * v.y, + 2.0 * v.x * v.y - 2.0 * s * v.z, 1.0 - 2.0 * v.x*v.x - 2.0 * v.z*v.z, 2.0 * v.y * v.z + 2.0 * s * v.x, + 2.0 * v.z * v.x + 2.0 * s * v.y, 2.0 * v.y * v.z - 2.0 * s * v.x, 1.0 - 2.0 * v.x*v.x - 2.0 * v.y*v.y); +} + + +/** Spherical linear interpolation (slerp) */ +Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t) { + scalar_t angle = acos(q1.s * q2.s + q1.v.x * q2.v.x + q1.v.y * q2.v.y + q1.v.z * q2.v.z); + scalar_t a = sin((1.0f - t) * angle); + scalar_t b = sin(t * angle); + scalar_t c = sin(angle); + + scalar_t x = (q1.v.x * a + q2.v.x * b) / c; + scalar_t y = (q1.v.y * a + q2.v.y * b) / c; + scalar_t z = (q1.v.z * a + q2.v.z * b) / c; + scalar_t s = (q1.s * a + q2.s * b) / c; + + return Quaternion(s, Vector3(x, y, z)).normalized(); +} + + + +std::ostream &operator <<(std::ostream &out, const Quaternion &q) { + out << "(" << q.s << ", " << q.v << ")"; + return out; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/quat.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/quat.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,119 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_QUATERNION_H_ +#define VMATH_QUATERNION_H_ + +#include +#include "vmath_types.h" +#include "vector.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +#define quat_cons(s, x, y, z) v4_cons(x, y, z, s) +#define quat_vec(q) v3_cons((q).x, (q).y, (q).z) +#define quat_s(q) ((q).w) +#define quat_identity() quat_cons(1.0, 0.0, 0.0, 0.0) +void quat_print(FILE *fp, quat_t q); + +#define quat_add v4_add +#define quat_sub v4_sub +#define quat_neg v4_neg + +static inline quat_t quat_mul(quat_t q1, quat_t q2); + +static inline quat_t quat_conjugate(quat_t q); + +#define quat_length v4_length +#define quat_length_sq v4_length_sq + +#define quat_normalize v4_normalize +static inline quat_t quat_inverse(quat_t q); + +quat_t quat_rotate(quat_t q, scalar_t angle, scalar_t x, scalar_t y, scalar_t z); +quat_t quat_rotate_quat(quat_t q, quat_t rotq); + +static inline void quat_to_mat3(mat3_t res, quat_t q); +static inline void quat_to_mat4(mat4_t res, quat_t q); + +#define quat_lerp quat_slerp +quat_t quat_slerp(quat_t q1, quat_t q2, scalar_t t); + + +#ifdef __cplusplus +} /* extern "C" */ + +#include + +/* Quaternion */ +class Quaternion { +public: + scalar_t s; + Vector3 v; + + Quaternion(); + Quaternion(scalar_t s, const Vector3 &v); + Quaternion(scalar_t s, scalar_t x, scalar_t y, scalar_t z); + Quaternion(const Vector3 &axis, scalar_t angle); + Quaternion(const quat_t &quat); + + Quaternion operator +(const Quaternion &quat) const; + Quaternion operator -(const Quaternion &quat) const; + Quaternion operator -() const; + Quaternion operator *(const Quaternion &quat) const; + + void operator +=(const Quaternion &quat); + void operator -=(const Quaternion &quat); + void operator *=(const Quaternion &quat); + + void reset_identity(); + + Quaternion conjugate() const; + + scalar_t length() const; + scalar_t length_sq() const; + + void normalize(); + Quaternion normalized() const; + + Quaternion inverse() const; + + void set_rotation(const Vector3 &axis, scalar_t angle); + void rotate(const Vector3 &axis, scalar_t angle); + /* note: this is a totally different operation from the above + * this treats the quaternion as signifying direction and rotates + * it by a rotation quaternion by rot * q * rot' + */ + void rotate(const Quaternion &q); + + Matrix3x3 get_rotation_matrix() const; + + friend Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t); + + friend std::ostream &operator <<(std::ostream &out, const Quaternion &q); +}; + +Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t); +inline Quaternion lerp(const Quaternion &q1, const Quaternion &q2, scalar_t t); +#endif /* __cplusplus */ + +#include "quat.inl" + +#endif /* VMATH_QUATERNION_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/quat.inl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/quat.inl Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,81 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include "vector.h" +#include "matrix.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +static inline quat_t quat_mul(quat_t q1, quat_t q2) +{ + quat_t res; + vec3_t v1 = quat_vec(q1); + vec3_t v2 = quat_vec(q2); + + res.w = q1.w * q2.w - v3_dot(v1, v2); + /* resvec = v2 * q1 + v1 * q2 + cross(v1, v2) */ + res.x = v2.x * q1.w + v1.x * q2.w + (v1.y * v2.z - v1.z * v2.y); + res.y = v2.y * q1.w + v1.y * q2.w + (v1.z * v2.x - v1.x * v2.z); + res.z = v2.z * q1.w + v1.z * q2.w + (v1.x * v2.y - v1.y * v2.x); + return res; +} + +static inline quat_t quat_conjugate(quat_t q) +{ + q.x = -q.x; + q.y = -q.y; + q.z = -q.z; + return q; +} + +static inline quat_t quat_inverse(quat_t q) +{ + scalar_t lensq = quat_length_sq(q); + q = quat_conjugate(q); + q.x /= lensq; + q.y /= lensq; + q.z /= lensq; + q.w /= lensq; + return q; +} + +static inline void quat_to_mat3(mat3_t res, quat_t q) +{ + m3_cons(res, 1.0 - 2.0 * q.y*q.y - 2.0 * q.z*q.z, 2.0 * q.x * q.y + 2.0 * q.w * q.z, 2.0 * q.z * q.x - 2.0 * q.w * q.y, + 2.0 * q.x * q.y - 2.0 * q.w * q.z, 1.0 - 2.0 * q.x*q.x - 2.0 * q.z*q.z, 2.0 * q.y * q.z + 2.0 * q.w * q.x, + 2.0 * q.z * q.x + 2.0 * q.w * q.y, 2.0 * q.y * q.z - 2.0 * q.w * q.x, 1.0 - 2.0 * q.x*q.x - 2.0 * q.y*q.y); +} + +static inline void quat_to_mat4(mat4_t res, quat_t q) +{ + m4_cons(res, 1.0 - 2.0 * q.y*q.y - 2.0 * q.z*q.z, 2.0 * q.x * q.y + 2.0 * q.w * q.z, 2.0 * q.z * q.x - 2.0 * q.w * q.y, 0, + 2.0 * q.x * q.y - 2.0 * q.w * q.z, 1.0 - 2.0 * q.x*q.x - 2.0 * q.z*q.z, 2.0 * q.y * q.z + 2.0 * q.w * q.x, 0, + 2.0 * q.z * q.x + 2.0 * q.w * q.y, 2.0 * q.y * q.z - 2.0 * q.w * q.x, 1.0 - 2.0 * q.x*q.x - 2.0 * q.y*q.y, 0, + 0, 0, 0, 1); +} + +#ifdef __cplusplus +} /* extern "C" */ + +inline Quaternion lerp(const Quaternion &a, const Quaternion &b, scalar_t t) +{ + return slerp(a, b, t); +} +#endif /* __cplusplus */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/quat_c.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/quat_c.c Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,61 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + + +#include +#include +#include "quat.h" + +void quat_print(FILE *fp, quat_t q) +{ + fprintf(fp, "([ %.4f %.4f %.4f ] %.4f)", q.x, q.y, q.z, q.w); +} + +quat_t quat_rotate(quat_t q, scalar_t angle, scalar_t x, scalar_t y, scalar_t z) +{ + quat_t rq; + scalar_t half_angle = angle * 0.5; + scalar_t sin_half = sin(half_angle); + + rq.w = cos(half_angle); + rq.x = x * sin_half; + rq.y = y * sin_half; + rq.z = z * sin_half; + + return quat_mul(q, rq); +} + +quat_t quat_rotate_quat(quat_t q, quat_t rotq) +{ + return quat_mul(quat_mul(rotq, q), quat_conjugate(rotq)); +} + +quat_t quat_slerp(quat_t q1, quat_t q2, scalar_t t) +{ + quat_t res; + scalar_t angle = acos(q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z); + scalar_t a = sin((1.0f - t) * angle); + scalar_t b = sin(t * angle); + scalar_t c = sin(angle); + + res.x = (q1.x * a + q2.x * b) / c; + res.y = (q1.y * a + q2.y * b) / c; + res.z = (q1.z * a + q2.z * b) / c; + res.w = (q1.w * a + q2.w * b) / c; + return quat_normalize(res); +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/ray.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/ray.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,66 @@ +#include "ray.h" +#include "vector.h" + +scalar_t Ray::env_ior = 1.0; + +Ray::Ray() { + ior = env_ior; + energy = 1.0; + time = 0; + iter = 0; +} + +Ray::Ray(const Vector3 &origin, const Vector3 &dir) { + this->origin = origin; + this->dir = dir; + ior = env_ior; + energy = 1.0; + time = 0; + iter = 0; +} + +void Ray::transform(const Matrix4x4 &xform) { + Matrix4x4 upper; + Vector3 dir; + + upper = xform; + upper[0][3] = upper[1][3] = upper[2][3] = upper[3][0] = upper[3][1] = upper[3][2] = 0.0; + upper[3][3] = 1.0; + + dir = (this->dir - origin).transformed(upper); + origin.transform(xform); + this->dir = dir + origin; +} + +Ray Ray::transformed(const Matrix4x4 &xform) const { + Ray foo = *this; + foo.transform(xform); + return foo; +} + +void Ray::enter(scalar_t new_ior) { + ior = new_ior; + ior_stack.push(ior); +} + +void Ray::leave() { + if(ior_stack.empty()) { + //std::cerr << "empty ior stack?\n"; + return; + } + ior_stack.pop(); + ior = ior_stack.empty() ? env_ior : ior_stack.top(); +} + +scalar_t Ray::calc_ior(bool entering, scalar_t mat_ior) const +{ + scalar_t from_ior = this->ior; + + if(entering) { + return from_ior / mat_ior; + } + + Ray tmp = *this; + tmp.leave(); + return from_ior / tmp.ior; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/ray.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/ray.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,73 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_RAY_H_ +#define VMATH_RAY_H_ + +#include "matrix.h" +#include "vector.h" + +typedef struct { + vec3_t origin, dir; +} ray_t; + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +static inline ray_t ray_cons(vec3_t origin, vec3_t dir); +ray_t ray_transform(ray_t r, mat4_t m); + +#ifdef __cplusplus +} /* __cplusplus */ + +#include + +class Ray { +private: + std::stack ior_stack; + +public: + /* enviornmental index of refraction, normally 1.0 */ + static scalar_t env_ior; + + Vector3 origin, dir; + scalar_t energy; + int iter; + scalar_t ior; + long time; + + Ray(); + Ray(const Vector3 &origin, const Vector3 &dir); + + void transform(const Matrix4x4 &xform); + Ray transformed(const Matrix4x4 &xform) const; + + void enter(scalar_t new_ior); + void leave(); + + scalar_t calc_ior(bool entering, scalar_t mat_ior = 1.0) const; +}; + +inline Ray reflect_ray(const Ray &inray, const Vector3 &norm); +inline Ray refract_ray(const Ray &inray, const Vector3 &norm, scalar_t ior, bool entering, scalar_t ray_mag = -1.0); +#endif /* __cplusplus */ + +#include "ray.inl" + +#endif /* VMATH_RAY_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/ray.inl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/ray.inl Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,70 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +static inline ray_t ray_cons(vec3_t origin, vec3_t dir) +{ + ray_t r; + r.origin = origin; + r.dir = dir; + return r; +} + +#ifdef __cplusplus +} + +inline Ray reflect_ray(const Ray &inray, const Vector3 &norm) +{ + Ray ray = inray; + ray.iter--; + ray.dir = ray.dir.reflection(norm); + return ray; +} + +inline Ray refract_ray(const Ray &inray, const Vector3 &norm, scalar_t mat_ior, bool entering, scalar_t ray_mag) +{ + Ray ray = inray; + ray.iter--; + + scalar_t ior = ray.calc_ior(entering, mat_ior); + + if(entering) { + ray.enter(mat_ior); + } else { + ray.leave(); + } + + if(ray_mag < 0.0) { + ray_mag = ray.dir.length(); + } + ray.dir = (ray.dir / ray_mag).refraction(norm, ior) * ray_mag; + + /* check TIR */ + if(dot_product(ray.dir, norm) > 0.0) { + if(entering) { + ray.leave(); + } else { + ray.enter(mat_ior); + } + } + return ray; +} +#endif /* __cplusplus */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/ray_c.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/ray_c.c Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,36 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include "ray.h" +#include "vector.h" + +ray_t ray_transform(ray_t r, mat4_t xform) +{ + mat4_t upper; + vec3_t dir; + + m4_copy(upper, xform); + upper[0][3] = upper[1][3] = upper[2][3] = upper[3][0] = upper[3][1] = upper[3][2] = 0.0; + upper[3][3] = 1.0; + + dir = v3_sub(r.dir, r.origin); + dir = v3_transform(dir, upper); + r.origin = v3_transform(r.origin, xform); + r.dir = v3_add(dir, r.origin); + return r; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/sphvec.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/sphvec.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,27 @@ +#include "sphvec.h" +#include "vector.h" + +/* theta: 0 <= theta <= 2pi, the angle around Y axis. + * phi: 0 <= phi <= pi, the angle from Y axis. + * r: radius. + */ +SphVector::SphVector(scalar_t theta, scalar_t phi, scalar_t r) { + this->theta = theta; + this->phi = phi; + this->r = r; +} + +/* Constructs a spherical coordinate vector from a cartesian vector */ +SphVector::SphVector(const Vector3 &cvec) { + *this = cvec; +} + +/* Assignment operator that converts cartesian to spherical coords */ +SphVector &SphVector::operator =(const Vector3 &cvec) { + r = cvec.length(); + //theta = atan2(cvec.y, cvec.x); + theta = atan2(cvec.z, cvec.x); + //phi = acos(cvec.z / r); + phi = acos(cvec.y / r); + return *this; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/sphvec.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/sphvec.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,36 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_SPHVEC_H_ +#define VMATH_SPHVEC_H_ + +#include "vmath_types.h" + +#ifdef __cplusplus +/* Vector in spherical coordinates */ +class SphVector { +public: + scalar_t theta, phi, r; + + SphVector(scalar_t theta = 0.0, scalar_t phi = 0.0, scalar_t r = 1.0); + SphVector(const Vector3 &cvec); + SphVector &operator =(const Vector3 &cvec); +}; +#endif /* __cplusplus */ + +#endif /* VMATH_SPHVEC_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vector.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vector.cc Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,326 @@ +#include "vector.h" +#include "vmath.h" + +// ---------- Vector2 ----------- + +Vector2::Vector2(scalar_t x, scalar_t y) +{ + this->x = x; + this->y = y; +} + +Vector2::Vector2(const vec2_t &vec) +{ + x = vec.x; + y = vec.y; +} + +Vector2::Vector2(const Vector3 &vec) +{ + x = vec.x; + y = vec.y; +} + +Vector2::Vector2(const Vector4 &vec) +{ + x = vec.x; + y = vec.y; +} + +void Vector2::normalize() +{ + scalar_t len = length(); + x /= len; + y /= len; +} + +Vector2 Vector2::normalized() const +{ + scalar_t len = length(); + return Vector2(x / len, y / len); +} + +void Vector2::transform(const Matrix3x3 &mat) +{ + scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2]; + y = mat[1][0] * x + mat[1][1] * y + mat[1][2]; + x = nx; +} + +Vector2 Vector2::transformed(const Matrix3x3 &mat) const +{ + Vector2 vec; + vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2]; + vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2]; + return vec; +} + +void Vector2::rotate(scalar_t angle) +{ + *this = Vector2(cos(angle) * x - sin(angle) * y, sin(angle) * x + cos(angle) * y); +} + +Vector2 Vector2::rotated(scalar_t angle) const +{ + return Vector2(cos(angle) * x - sin(angle) * y, sin(angle) * x + cos(angle) * y); +} + +Vector2 Vector2::reflection(const Vector2 &normal) const +{ + return 2.0 * dot_product(*this, normal) * normal - *this; +} + +Vector2 Vector2::refraction(const Vector2 &normal, scalar_t src_ior, scalar_t dst_ior) const +{ + // quick and dirty implementation :) + Vector3 v3refr = Vector3(this->x, this->y, 1.0).refraction(Vector3(this->x, this->y, 1), src_ior, dst_ior); + return Vector2(v3refr.x, v3refr.y); +} + +std::ostream &operator <<(std::ostream &out, const Vector2 &vec) +{ + out << "[" << vec.x << " " << vec.y << "]"; + return out; +} + + + +// --------- Vector3 ---------- + +Vector3::Vector3(scalar_t x, scalar_t y, scalar_t z) +{ + this->x = x; + this->y = y; + this->z = z; +} + +Vector3::Vector3(const vec3_t &vec) +{ + x = vec.x; + y = vec.y; + z = vec.z; +} + +Vector3::Vector3(const Vector2 &vec) +{ + x = vec.x; + y = vec.y; + z = 1; +} + +Vector3::Vector3(const Vector4 &vec) +{ + x = vec.x; + y = vec.y; + z = vec.z; +} + +Vector3::Vector3(const SphVector &sph) +{ + *this = sph; +} + +Vector3 &Vector3::operator =(const SphVector &sph) +{ + x = sph.r * cos(sph.theta) * sin(sph.phi); + z = sph.r * sin(sph.theta) * sin(sph.phi); + y = sph.r * cos(sph.phi); + return *this; +} + +void Vector3::normalize() +{ + scalar_t len = length(); + x /= len; + y /= len; + z /= len; +} + +Vector3 Vector3::normalized() const +{ + scalar_t len = length(); + return Vector3(x / len, y / len, z / len); +} + +Vector3 Vector3::reflection(const Vector3 &normal) const +{ + return 2.0 * dot_product(*this, normal) * normal - *this; +} + +Vector3 Vector3::refraction(const Vector3 &normal, scalar_t src_ior, scalar_t dst_ior) const +{ + return refraction(normal, src_ior / dst_ior); +} + +Vector3 Vector3::refraction(const Vector3 &normal, scalar_t ior) const +{ + scalar_t cos_inc = dot_product(*this, -normal); + + scalar_t radical = 1.0 + SQ(ior) * (SQ(cos_inc) - 1.0); + + if(radical < 0.0) { // total internal reflection + return -reflection(normal); + } + + scalar_t beta = ior * cos_inc - sqrt(radical); + + return *this * ior + normal * beta; +} + +void Vector3::transform(const Matrix3x3 &mat) +{ + scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z; + scalar_t ny = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z; + z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z; + x = nx; + y = ny; +} + +Vector3 Vector3::transformed(const Matrix3x3 &mat) const +{ + Vector3 vec; + vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z; + vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z; + vec.z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z; + return vec; +} + +void Vector3::transform(const Matrix4x4 &mat) +{ + scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3]; + scalar_t ny = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3]; + z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3]; + x = nx; + y = ny; +} + +Vector3 Vector3::transformed(const Matrix4x4 &mat) const +{ + Vector3 vec; + vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3]; + vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3]; + vec.z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3]; + return vec; +} + +void Vector3::transform(const Quaternion &quat) +{ + Quaternion vq(0.0f, *this); + vq = quat * vq * quat.inverse(); + *this = vq.v; +} + +Vector3 Vector3::transformed(const Quaternion &quat) const +{ + Quaternion vq(0.0f, *this); + vq = quat * vq * quat.inverse(); + return vq.v; +} + +void Vector3::rotate(const Vector3 &euler) +{ + Matrix4x4 rot; + rot.set_rotation(euler); + transform(rot); +} + +Vector3 Vector3::rotated(const Vector3 &euler) const +{ + Matrix4x4 rot; + rot.set_rotation(euler); + return transformed(rot); +} + +std::ostream &operator <<(std::ostream &out, const Vector3 &vec) +{ + out << "[" << vec.x << " " << vec.y << " " << vec.z << "]"; + return out; +} + + +// -------------- Vector4 -------------- +Vector4::Vector4(scalar_t x, scalar_t y, scalar_t z, scalar_t w) +{ + this->x = x; + this->y = y; + this->z = z; + this->w = w; +} + +Vector4::Vector4(const vec4_t &vec) +{ + x = vec.x; + y = vec.y; + z = vec.z; + w = vec.w; +} + +Vector4::Vector4(const Vector2 &vec) +{ + x = vec.x; + y = vec.y; + z = 1; + w = 1; +} + +Vector4::Vector4(const Vector3 &vec) +{ + x = vec.x; + y = vec.y; + z = vec.z; + w = 1; +} + +void Vector4::normalize() +{ + scalar_t len = (scalar_t)sqrt(x*x + y*y + z*z + w*w); + x /= len; + y /= len; + z /= len; + w /= len; +} + +Vector4 Vector4::normalized() const +{ + scalar_t len = (scalar_t)sqrt(x*x + y*y + z*z + w*w); + return Vector4(x / len, y / len, z / len, w / len); +} + +void Vector4::transform(const Matrix4x4 &mat) +{ + scalar_t nx = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3] * w; + scalar_t ny = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3] * w; + scalar_t nz = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3] * w; + w = mat[3][0] * x + mat[3][1] * y + mat[3][2] * z + mat[3][3] * w; + x = nx; + y = ny; + z = nz; +} + +Vector4 Vector4::transformed(const Matrix4x4 &mat) const +{ + Vector4 vec; + vec.x = mat[0][0] * x + mat[0][1] * y + mat[0][2] * z + mat[0][3] * w; + vec.y = mat[1][0] * x + mat[1][1] * y + mat[1][2] * z + mat[1][3] * w; + vec.z = mat[2][0] * x + mat[2][1] * y + mat[2][2] * z + mat[2][3] * w; + vec.w = mat[3][0] * x + mat[3][1] * y + mat[3][2] * z + mat[3][3] * w; + return vec; +} + +// TODO: implement 4D vector reflection +Vector4 Vector4::reflection(const Vector4 &normal) const +{ + return *this; +} + +// TODO: implement 4D vector refraction +Vector4 Vector4::refraction(const Vector4 &normal, scalar_t src_ior, scalar_t dst_ior) const +{ + return *this; +} + +std::ostream &operator <<(std::ostream &out, const Vector4 &vec) +{ + out << "[" << vec.x << " " << vec.y << " " << vec.z << " " << vec.w << "]"; + return out; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vector.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vector.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,292 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_VECTOR_H_ +#define VMATH_VECTOR_H_ + +#include +#include "vmath_types.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +/* C 2D vector functions */ +static inline vec2_t v2_cons(scalar_t x, scalar_t y); +static inline void v2_print(FILE *fp, vec2_t v); + +static inline vec2_t v2_add(vec2_t v1, vec2_t v2); +static inline vec2_t v2_sub(vec2_t v1, vec2_t v2); +static inline vec2_t v2_scale(vec2_t v, scalar_t s); +static inline scalar_t v2_dot(vec2_t v1, vec2_t v2); +static inline scalar_t v2_length(vec2_t v); +static inline scalar_t v2_length_sq(vec2_t v); +static inline vec2_t v2_normalize(vec2_t v); + +static inline vec2_t v2_lerp(vec2_t v1, vec2_t v2, scalar_t t); + +/* C 3D vector functions */ +static inline vec3_t v3_cons(scalar_t x, scalar_t y, scalar_t z); +static inline void v3_print(FILE *fp, vec3_t v); + +static inline vec3_t v3_add(vec3_t v1, vec3_t v2); +static inline vec3_t v3_sub(vec3_t v1, vec3_t v2); +static inline vec3_t v3_neg(vec3_t v); +static inline vec3_t v3_mul(vec3_t v1, vec3_t v2); +static inline vec3_t v3_scale(vec3_t v1, scalar_t s); +static inline scalar_t v3_dot(vec3_t v1, vec3_t v2); +static inline vec3_t v3_cross(vec3_t v1, vec3_t v2); +static inline scalar_t v3_length(vec3_t v); +static inline scalar_t v3_length_sq(vec3_t v); +static inline vec3_t v3_normalize(vec3_t v); +static inline vec3_t v3_transform(vec3_t v, mat4_t m); + +static inline vec3_t v3_rotate(vec3_t v, scalar_t x, scalar_t y, scalar_t z); +static inline vec3_t v3_rotate_axis(vec3_t v, scalar_t angle, scalar_t x, scalar_t y, scalar_t z); +static inline vec3_t v3_rotate_quat(vec3_t v, quat_t q); + +static inline vec3_t v3_reflect(vec3_t v, vec3_t n); + +static inline vec3_t v3_lerp(vec3_t v1, vec3_t v2, scalar_t t); + +/* C 4D vector functions */ +static inline vec4_t v4_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t w); +static inline void v4_print(FILE *fp, vec4_t v); + +static inline vec4_t v4_add(vec4_t v1, vec4_t v2); +static inline vec4_t v4_sub(vec4_t v1, vec4_t v2); +static inline vec4_t v4_neg(vec4_t v); +static inline vec4_t v4_mul(vec4_t v1, vec4_t v2); +static inline vec4_t v4_scale(vec4_t v, scalar_t s); +static inline scalar_t v4_dot(vec4_t v1, vec4_t v2); +static inline scalar_t v4_length(vec4_t v); +static inline scalar_t v4_length_sq(vec4_t v); +static inline vec4_t v4_normalize(vec4_t v); +static inline vec4_t v4_transform(vec4_t v, mat4_t m); + +#ifdef __cplusplus +} /* extern "C" */ + +/* when included from C++ source files, also define the vector classes */ +#include + +/** 2D Vector */ +class Vector2 { +public: + scalar_t x, y; + + Vector2(scalar_t x = 0.0, scalar_t y = 0.0); + Vector2(const vec2_t &vec); + Vector2(const Vector3 &vec); + Vector2(const Vector4 &vec); + + inline scalar_t &operator [](int elem); + inline const scalar_t &operator [](int elem) const; + + inline scalar_t length() const; + inline scalar_t length_sq() const; + void normalize(); + Vector2 normalized() const; + + void transform(const Matrix3x3 &mat); + Vector2 transformed(const Matrix3x3 &mat) const; + + void rotate(scalar_t angle); + Vector2 rotated(scalar_t angle) const; + + Vector2 reflection(const Vector2 &normal) const; + Vector2 refraction(const Vector2 &normal, scalar_t src_ior, scalar_t dst_ior) const; +}; + +/* unary operations */ +inline Vector2 operator -(const Vector2 &vec); + +/* binary vector (op) vector operations */ +inline scalar_t dot_product(const Vector2 &v1, const Vector2 &v2); + +inline Vector2 operator +(const Vector2 &v1, const Vector2 &v2); +inline Vector2 operator -(const Vector2 &v1, const Vector2 &v2); +inline Vector2 operator *(const Vector2 &v1, const Vector2 &v2); +inline Vector2 operator /(const Vector2 &v1, const Vector2 &v2); +inline bool operator ==(const Vector2 &v1, const Vector2 &v2); + +inline void operator +=(Vector2 &v1, const Vector2 &v2); +inline void operator -=(Vector2 &v1, const Vector2 &v2); +inline void operator *=(Vector2 &v1, const Vector2 &v2); +inline void operator /=(Vector2 &v1, const Vector2 &v2); + +/* binary vector (op) scalar and scalar (op) vector operations */ +inline Vector2 operator +(const Vector2 &vec, scalar_t scalar); +inline Vector2 operator +(scalar_t scalar, const Vector2 &vec); +inline Vector2 operator -(const Vector2 &vec, scalar_t scalar); +inline Vector2 operator *(const Vector2 &vec, scalar_t scalar); +inline Vector2 operator *(scalar_t scalar, const Vector2 &vec); +inline Vector2 operator /(const Vector2 &vec, scalar_t scalar); + +inline void operator +=(Vector2 &vec, scalar_t scalar); +inline void operator -=(Vector2 &vec, scalar_t scalar); +inline void operator *=(Vector2 &vec, scalar_t scalar); +inline void operator /=(Vector2 &vec, scalar_t scalar); + +std::ostream &operator <<(std::ostream &out, const Vector2 &vec); + +inline Vector2 lerp(const Vector2 &a, const Vector2 &b, scalar_t t); +inline Vector2 catmull_rom_spline(const Vector2 &v0, const Vector2 &v1, + const Vector2 &v2, const Vector2 &v3, scalar_t t); + +/* 3D Vector */ +class Vector3 { +public: + scalar_t x, y, z; + + Vector3(scalar_t x = 0.0, scalar_t y = 0.0, scalar_t z = 0.0); + Vector3(const vec3_t &vec); + Vector3(const Vector2 &vec); + Vector3(const Vector4 &vec); + Vector3(const SphVector &sph); + + Vector3 &operator =(const SphVector &sph); + + inline scalar_t &operator [](int elem); + inline const scalar_t &operator [](int elem) const; + + inline scalar_t length() const; + inline scalar_t length_sq() const; + void normalize(); + Vector3 normalized() const; + + void transform(const Matrix3x3 &mat); + Vector3 transformed(const Matrix3x3 &mat) const; + void transform(const Matrix4x4 &mat); + Vector3 transformed(const Matrix4x4 &mat) const; + void transform(const Quaternion &quat); + Vector3 transformed(const Quaternion &quat) const; + + void rotate(const Vector3 &euler); + Vector3 rotated(const Vector3 &euler) const; + + Vector3 reflection(const Vector3 &normal) const; + Vector3 refraction(const Vector3 &normal, scalar_t src_ior, scalar_t dst_ior) const; + Vector3 refraction(const Vector3 &normal, scalar_t ior) const; +}; + +/* unary operations */ +inline Vector3 operator -(const Vector3 &vec); + +/* binary vector (op) vector operations */ +inline scalar_t dot_product(const Vector3 &v1, const Vector3 &v2); +inline Vector3 cross_product(const Vector3 &v1, const Vector3 &v2); + +inline Vector3 operator +(const Vector3 &v1, const Vector3 &v2); +inline Vector3 operator -(const Vector3 &v1, const Vector3 &v2); +inline Vector3 operator *(const Vector3 &v1, const Vector3 &v2); +inline Vector3 operator /(const Vector3 &v1, const Vector3 &v2); +inline bool operator ==(const Vector3 &v1, const Vector3 &v2); + +inline void operator +=(Vector3 &v1, const Vector3 &v2); +inline void operator -=(Vector3 &v1, const Vector3 &v2); +inline void operator *=(Vector3 &v1, const Vector3 &v2); +inline void operator /=(Vector3 &v1, const Vector3 &v2); + +/* binary vector (op) scalar and scalar (op) vector operations */ +inline Vector3 operator +(const Vector3 &vec, scalar_t scalar); +inline Vector3 operator +(scalar_t scalar, const Vector3 &vec); +inline Vector3 operator -(const Vector3 &vec, scalar_t scalar); +inline Vector3 operator *(const Vector3 &vec, scalar_t scalar); +inline Vector3 operator *(scalar_t scalar, const Vector3 &vec); +inline Vector3 operator /(const Vector3 &vec, scalar_t scalar); + +inline void operator +=(Vector3 &vec, scalar_t scalar); +inline void operator -=(Vector3 &vec, scalar_t scalar); +inline void operator *=(Vector3 &vec, scalar_t scalar); +inline void operator /=(Vector3 &vec, scalar_t scalar); + +std::ostream &operator <<(std::ostream &out, const Vector3 &vec); + +inline Vector3 lerp(const Vector3 &a, const Vector3 &b, scalar_t t); +inline Vector3 catmull_rom_spline(const Vector3 &v0, const Vector3 &v1, + const Vector3 &v2, const Vector3 &v3, scalar_t t); + +/* 4D Vector */ +class Vector4 { +public: + scalar_t x, y, z, w; + + Vector4(scalar_t x = 0.0, scalar_t y = 0.0, scalar_t z = 0.0, scalar_t w = 0.0); + Vector4(const vec4_t &vec); + Vector4(const Vector2 &vec); + Vector4(const Vector3 &vec); + + inline scalar_t &operator [](int elem); + inline const scalar_t &operator [](int elem) const; + + inline scalar_t length() const; + inline scalar_t length_sq() const; + void normalize(); + Vector4 normalized() const; + + void transform(const Matrix4x4 &mat); + Vector4 transformed(const Matrix4x4 &mat) const; + + Vector4 reflection(const Vector4 &normal) const; + Vector4 refraction(const Vector4 &normal, scalar_t src_ior, scalar_t dst_ior) const; +}; + + +/* unary operations */ +inline Vector4 operator -(const Vector4 &vec); + +/* binary vector (op) vector operations */ +inline scalar_t dot_product(const Vector4 &v1, const Vector4 &v2); +inline Vector4 cross_product(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3); + +inline Vector4 operator +(const Vector4 &v1, const Vector4 &v2); +inline Vector4 operator -(const Vector4 &v1, const Vector4 &v2); +inline Vector4 operator *(const Vector4 &v1, const Vector4 &v2); +inline Vector4 operator /(const Vector4 &v1, const Vector4 &v2); +inline bool operator ==(const Vector4 &v1, const Vector4 &v2); + +inline void operator +=(Vector4 &v1, const Vector4 &v2); +inline void operator -=(Vector4 &v1, const Vector4 &v2); +inline void operator *=(Vector4 &v1, const Vector4 &v2); +inline void operator /=(Vector4 &v1, const Vector4 &v2); + +/* binary vector (op) scalar and scalar (op) vector operations */ +inline Vector4 operator +(const Vector4 &vec, scalar_t scalar); +inline Vector4 operator +(scalar_t scalar, const Vector4 &vec); +inline Vector4 operator -(const Vector4 &vec, scalar_t scalar); +inline Vector4 operator *(const Vector4 &vec, scalar_t scalar); +inline Vector4 operator *(scalar_t scalar, const Vector4 &vec); +inline Vector4 operator /(const Vector4 &vec, scalar_t scalar); + +inline void operator +=(Vector4 &vec, scalar_t scalar); +inline void operator -=(Vector4 &vec, scalar_t scalar); +inline void operator *=(Vector4 &vec, scalar_t scalar); +inline void operator /=(Vector4 &vec, scalar_t scalar); + +std::ostream &operator <<(std::ostream &out, const Vector4 &vec); + +inline Vector4 lerp(const Vector4 &v0, const Vector4 &v1, scalar_t t); +inline Vector4 catmull_rom_spline(const Vector4 &v0, const Vector4 &v1, + const Vector4 &v2, const Vector4 &v3, scalar_t t); + +#endif /* __cplusplus */ + +#include "vector.inl" + +#endif /* VMATH_VECTOR_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vector.inl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vector.inl Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,761 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +/* C 2D vector functions */ +static inline vec2_t v2_cons(scalar_t x, scalar_t y) +{ + vec2_t v; + v.x = x; + v.y = y; + return v; +} + +static inline void v2_print(FILE *fp, vec2_t v) +{ + fprintf(fp, "[ %.4f %.4f ]", v.x, v.y); +} + +static inline vec2_t v2_add(vec2_t v1, vec2_t v2) +{ + vec2_t res; + res.x = v1.x + v2.x; + res.y = v1.y + v2.y; + return res; +} + +static inline vec2_t v2_sub(vec2_t v1, vec2_t v2) +{ + vec2_t res; + res.x = v1.x - v2.x; + res.y = v1.y - v2.y; + return res; +} + +static inline vec2_t v2_scale(vec2_t v, scalar_t s) +{ + vec2_t res; + res.x = v.x * s; + res.y = v.y * s; + return res; +} + +static inline scalar_t v2_dot(vec2_t v1, vec2_t v2) +{ + return v1.x * v2.x + v1.y * v2.y; +} + +static inline scalar_t v2_length(vec2_t v) +{ + return sqrt(v.x * v.x + v.y * v.y); +} + +static inline scalar_t v2_length_sq(vec2_t v) +{ + return v.x * v.x + v.y * v.y; +} + +static inline vec2_t v2_normalize(vec2_t v) +{ + scalar_t len = (scalar_t)sqrt(v.x * v.x + v.y * v.y); + v.x /= len; + v.y /= len; + return v; +} + +static inline vec2_t v2_lerp(vec2_t v1, vec2_t v2, scalar_t t) +{ + vec2_t res; + res.x = v1.x + (v2.x - v1.x) * t; + res.y = v1.y + (v2.y - v1.y) * t; + return res; +} + + +/* C 3D vector functions */ +static inline vec3_t v3_cons(scalar_t x, scalar_t y, scalar_t z) +{ + vec3_t v; + v.x = x; + v.y = y; + v.z = z; + return v; +} + +static inline void v3_print(FILE *fp, vec3_t v) +{ + fprintf(fp, "[ %.4f %.4f %.4f ]", v.x, v.y, v.z); +} + +static inline vec3_t v3_add(vec3_t v1, vec3_t v2) +{ + v1.x += v2.x; + v1.y += v2.y; + v1.z += v2.z; + return v1; +} + +static inline vec3_t v3_sub(vec3_t v1, vec3_t v2) +{ + v1.x -= v2.x; + v1.y -= v2.y; + v1.z -= v2.z; + return v1; +} + +static inline vec3_t v3_neg(vec3_t v) +{ + v.x = -v.x; + v.y = -v.y; + v.z = -v.z; + return v; +} + +static inline vec3_t v3_mul(vec3_t v1, vec3_t v2) +{ + v1.x *= v2.x; + v1.y *= v2.y; + v1.z *= v2.z; + return v1; +} + +static inline vec3_t v3_scale(vec3_t v1, scalar_t s) +{ + v1.x *= s; + v1.y *= s; + v1.z *= s; + return v1; +} + +static inline scalar_t v3_dot(vec3_t v1, vec3_t v2) +{ + return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; +} + +static inline vec3_t v3_cross(vec3_t v1, vec3_t v2) +{ + vec3_t v; + v.x = v1.y * v2.z - v1.z * v2.y; + v.y = v1.z * v2.x - v1.x * v2.z; + v.z = v1.x * v2.y - v1.y * v2.x; + return v; +} + +static inline scalar_t v3_length(vec3_t v) +{ + return sqrt(v.x * v.x + v.y * v.y + v.z * v.z); +} + +static inline scalar_t v3_length_sq(vec3_t v) +{ + return v.x * v.x + v.y * v.y + v.z * v.z; +} + +static inline vec3_t v3_normalize(vec3_t v) +{ + scalar_t len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z); + v.x /= len; + v.y /= len; + v.z /= len; + return v; +} + +static inline vec3_t v3_transform(vec3_t v, mat4_t m) +{ + vec3_t res; + res.x = m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z + m[0][3]; + res.y = m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z + m[1][3]; + res.z = m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z + m[2][3]; + return res; +} + +static inline vec3_t v3_rotate(vec3_t v, scalar_t x, scalar_t y, scalar_t z) +{ + void m4_rotate(mat4_t, scalar_t, scalar_t, scalar_t); + + mat4_t m = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}; + m4_rotate(m, x, y, z); + return v3_transform(v, m); +} + +static inline vec3_t v3_rotate_axis(vec3_t v, scalar_t angle, scalar_t x, scalar_t y, scalar_t z) +{ + void m4_rotate_axis(mat4_t, scalar_t, scalar_t, scalar_t, scalar_t); + + mat4_t m = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}}; + m4_rotate_axis(m, angle, x, y, z); + return v3_transform(v, m); +} + +static inline vec3_t v3_rotate_quat(vec3_t v, quat_t q) +{ + quat_t quat_rotate_quat(quat_t, quat_t); + + quat_t vq = v4_cons(v.x, v.y, v.z, 0.0); + quat_t res = quat_rotate_quat(vq, q); + return v3_cons(res.x, res.y, res.z); +} + +static inline vec3_t v3_reflect(vec3_t v, vec3_t n) +{ + scalar_t dot = v3_dot(v, n); + return v3_sub(v3_scale(n, dot * 2.0), v); +} + +static inline vec3_t v3_lerp(vec3_t v1, vec3_t v2, scalar_t t) +{ + v1.x += (v2.x - v1.x) * t; + v1.y += (v2.y - v1.y) * t; + v1.z += (v2.z - v1.z) * t; + return v1; +} + +/* C 4D vector functions */ +static inline vec4_t v4_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t w) +{ + vec4_t v; + v.x = x; + v.y = y; + v.z = z; + v.w = w; + return v; +} + +static inline void v4_print(FILE *fp, vec4_t v) +{ + fprintf(fp, "[ %.4f %.4f %.4f %.4f ]", v.x, v.y, v.z, v.w); +} + +static inline vec4_t v4_add(vec4_t v1, vec4_t v2) +{ + v1.x += v2.x; + v1.y += v2.y; + v1.z += v2.z; + v1.w += v2.w; + return v1; +} + +static inline vec4_t v4_sub(vec4_t v1, vec4_t v2) +{ + v1.x -= v2.x; + v1.y -= v2.y; + v1.z -= v2.z; + v1.w -= v2.w; + return v1; +} + +static inline vec4_t v4_neg(vec4_t v) +{ + v.x = -v.x; + v.y = -v.y; + v.z = -v.z; + v.w = -v.w; + return v; +} + +static inline vec4_t v4_mul(vec4_t v1, vec4_t v2) +{ + v1.x *= v2.x; + v1.y *= v2.y; + v1.z *= v2.z; + v1.w *= v2.w; + return v1; +} + +static inline vec4_t v4_scale(vec4_t v, scalar_t s) +{ + v.x *= s; + v.y *= s; + v.z *= s; + v.w *= s; + return v; +} + +static inline scalar_t v4_dot(vec4_t v1, vec4_t v2) +{ + return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + v1.w * v2.w; +} + +static inline scalar_t v4_length(vec4_t v) +{ + return sqrt(v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w); +} + +static inline scalar_t v4_length_sq(vec4_t v) +{ + return v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w; +} + +static inline vec4_t v4_normalize(vec4_t v) +{ + scalar_t len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z + v.w * v.w); + v.x /= len; + v.y /= len; + v.z /= len; + v.w /= len; + return v; +} + +static inline vec4_t v4_transform(vec4_t v, mat4_t m) +{ + vec4_t res; + res.x = m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z + m[0][3] * v.w; + res.y = m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z + m[1][3] * v.w; + res.z = m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z + m[2][3] * v.w; + res.w = m[3][0] * v.x + m[3][1] * v.y + m[3][2] * v.z + m[3][3] * v.w; + return res; +} + +#ifdef __cplusplus +} /* extern "C" */ + + +/* --------------- C++ part -------------- */ + +inline scalar_t &Vector2::operator [](int elem) { + return elem ? y : x; +} + +inline const scalar_t &Vector2::operator [](int elem) const { + return elem ? y : x; +} + +inline Vector2 operator -(const Vector2 &vec) { + return Vector2(-vec.x, -vec.y); +} + +inline scalar_t dot_product(const Vector2 &v1, const Vector2 &v2) { + return v1.x * v2.x + v1.y * v2.y; +} + +inline Vector2 operator +(const Vector2 &v1, const Vector2 &v2) { + return Vector2(v1.x + v2.x, v1.y + v2.y); +} + +inline Vector2 operator -(const Vector2 &v1, const Vector2 &v2) { + return Vector2(v1.x - v2.x, v1.y - v2.y); +} + +inline Vector2 operator *(const Vector2 &v1, const Vector2 &v2) { + return Vector2(v1.x * v2.x, v1.y * v2.y); +} + +inline Vector2 operator /(const Vector2 &v1, const Vector2 &v2) { + return Vector2(v1.x / v2.x, v1.y / v2.y); +} + +inline bool operator ==(const Vector2 &v1, const Vector2 &v2) { + return (fabs(v1.x - v2.x) < XSMALL_NUMBER) && (fabs(v1.y - v2.x) < XSMALL_NUMBER); +} + +inline void operator +=(Vector2 &v1, const Vector2 &v2) { + v1.x += v2.x; + v1.y += v2.y; +} + +inline void operator -=(Vector2 &v1, const Vector2 &v2) { + v1.x -= v2.x; + v1.y -= v2.y; +} + +inline void operator *=(Vector2 &v1, const Vector2 &v2) { + v1.x *= v2.x; + v1.y *= v2.y; +} + +inline void operator /=(Vector2 &v1, const Vector2 &v2) { + v1.x /= v2.x; + v1.y /= v2.y; +} + +inline Vector2 operator +(const Vector2 &vec, scalar_t scalar) { + return Vector2(vec.x + scalar, vec.y + scalar); +} + +inline Vector2 operator +(scalar_t scalar, const Vector2 &vec) { + return Vector2(vec.x + scalar, vec.y + scalar); +} + +inline Vector2 operator -(scalar_t scalar, const Vector2 &vec) { + return Vector2(vec.x - scalar, vec.y - scalar); +} + +inline Vector2 operator *(const Vector2 &vec, scalar_t scalar) { + return Vector2(vec.x * scalar, vec.y * scalar); +} + +inline Vector2 operator *(scalar_t scalar, const Vector2 &vec) { + return Vector2(vec.x * scalar, vec.y * scalar); +} + +inline Vector2 operator /(const Vector2 &vec, scalar_t scalar) { + return Vector2(vec.x / scalar, vec.y / scalar); +} + +inline void operator +=(Vector2 &vec, scalar_t scalar) { + vec.x += scalar; + vec.y += scalar; +} + +inline void operator -=(Vector2 &vec, scalar_t scalar) { + vec.x -= scalar; + vec.y -= scalar; +} + +inline void operator *=(Vector2 &vec, scalar_t scalar) { + vec.x *= scalar; + vec.y *= scalar; +} + +inline void operator /=(Vector2 &vec, scalar_t scalar) { + vec.x /= scalar; + vec.y /= scalar; +} + +inline scalar_t Vector2::length() const { + return sqrt(x*x + y*y); +} + +inline scalar_t Vector2::length_sq() const { + return x*x + y*y; +} + +inline Vector2 lerp(const Vector2 &a, const Vector2 &b, scalar_t t) +{ + return a + (b - a) * t; +} + +inline Vector2 catmull_rom_spline(const Vector2 &v0, const Vector2 &v1, + const Vector2 &v2, const Vector2 &v3, scalar_t t) +{ + scalar_t spline(scalar_t, scalar_t, scalar_t, scalar_t, scalar_t); + scalar_t x = spline(v0.x, v1.x, v2.x, v3.x, t); + scalar_t y = spline(v0.y, v1.y, v2.y, v3.y, t); + return Vector2(x, y); +} + + +/* ------------- Vector3 -------------- */ + +inline scalar_t &Vector3::operator [](int elem) { + return elem ? (elem == 1 ? y : z) : x; +} + +inline const scalar_t &Vector3::operator [](int elem) const { + return elem ? (elem == 1 ? y : z) : x; +} + +/* unary operations */ +inline Vector3 operator -(const Vector3 &vec) { + return Vector3(-vec.x, -vec.y, -vec.z); +} + +/* binary vector (op) vector operations */ +inline scalar_t dot_product(const Vector3 &v1, const Vector3 &v2) { + return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; +} + +inline Vector3 cross_product(const Vector3 &v1, const Vector3 &v2) { + return Vector3(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x); +} + + +inline Vector3 operator +(const Vector3 &v1, const Vector3 &v2) { + return Vector3(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z); +} + +inline Vector3 operator -(const Vector3 &v1, const Vector3 &v2) { + return Vector3(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z); +} + +inline Vector3 operator *(const Vector3 &v1, const Vector3 &v2) { + return Vector3(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z); +} + +inline Vector3 operator /(const Vector3 &v1, const Vector3 &v2) { + return Vector3(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z); +} + +inline bool operator ==(const Vector3 &v1, const Vector3 &v2) { + return (fabs(v1.x - v2.x) < XSMALL_NUMBER) && (fabs(v1.y - v2.y) < XSMALL_NUMBER) && (fabs(v1.z - v2.z) < XSMALL_NUMBER); +} + +inline void operator +=(Vector3 &v1, const Vector3 &v2) { + v1.x += v2.x; + v1.y += v2.y; + v1.z += v2.z; +} + +inline void operator -=(Vector3 &v1, const Vector3 &v2) { + v1.x -= v2.x; + v1.y -= v2.y; + v1.z -= v2.z; +} + +inline void operator *=(Vector3 &v1, const Vector3 &v2) { + v1.x *= v2.x; + v1.y *= v2.y; + v1.z *= v2.z; +} + +inline void operator /=(Vector3 &v1, const Vector3 &v2) { + v1.x /= v2.x; + v1.y /= v2.y; + v1.z /= v2.z; +} +/* binary vector (op) scalar and scalar (op) vector operations */ +inline Vector3 operator +(const Vector3 &vec, scalar_t scalar) { + return Vector3(vec.x + scalar, vec.y + scalar, vec.z + scalar); +} + +inline Vector3 operator +(scalar_t scalar, const Vector3 &vec) { + return Vector3(vec.x + scalar, vec.y + scalar, vec.z + scalar); +} + +inline Vector3 operator -(const Vector3 &vec, scalar_t scalar) { + return Vector3(vec.x - scalar, vec.y - scalar, vec.z - scalar); +} + +inline Vector3 operator *(const Vector3 &vec, scalar_t scalar) { + return Vector3(vec.x * scalar, vec.y * scalar, vec.z * scalar); +} + +inline Vector3 operator *(scalar_t scalar, const Vector3 &vec) { + return Vector3(vec.x * scalar, vec.y * scalar, vec.z * scalar); +} + +inline Vector3 operator /(const Vector3 &vec, scalar_t scalar) { + return Vector3(vec.x / scalar, vec.y / scalar, vec.z / scalar); +} + +inline void operator +=(Vector3 &vec, scalar_t scalar) { + vec.x += scalar; + vec.y += scalar; + vec.z += scalar; +} + +inline void operator -=(Vector3 &vec, scalar_t scalar) { + vec.x -= scalar; + vec.y -= scalar; + vec.z -= scalar; +} + +inline void operator *=(Vector3 &vec, scalar_t scalar) { + vec.x *= scalar; + vec.y *= scalar; + vec.z *= scalar; +} + +inline void operator /=(Vector3 &vec, scalar_t scalar) { + vec.x /= scalar; + vec.y /= scalar; + vec.z /= scalar; +} + +inline scalar_t Vector3::length() const { + return sqrt(x*x + y*y + z*z); +} +inline scalar_t Vector3::length_sq() const { + return x*x + y*y + z*z; +} + +inline Vector3 lerp(const Vector3 &a, const Vector3 &b, scalar_t t) { + return a + (b - a) * t; +} + +inline Vector3 catmull_rom_spline(const Vector3 &v0, const Vector3 &v1, + const Vector3 &v2, const Vector3 &v3, scalar_t t) +{ + scalar_t spline(scalar_t, scalar_t, scalar_t, scalar_t, scalar_t); + scalar_t x = spline(v0.x, v1.x, v2.x, v3.x, t); + scalar_t y = spline(v0.y, v1.y, v2.y, v3.y, t); + scalar_t z = spline(v0.z, v1.z, v2.z, v3.z, t); + return Vector3(x, y, z); +} + +/* ----------- Vector4 ----------------- */ + +inline scalar_t &Vector4::operator [](int elem) { + return elem ? (elem == 1 ? y : (elem == 2 ? z : w)) : x; +} + +inline const scalar_t &Vector4::operator [](int elem) const { + return elem ? (elem == 1 ? y : (elem == 2 ? z : w)) : x; +} + +inline Vector4 operator -(const Vector4 &vec) { + return Vector4(-vec.x, -vec.y, -vec.z, -vec.w); +} + +inline scalar_t dot_product(const Vector4 &v1, const Vector4 &v2) { + return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z + v1.w * v2.w; +} + +inline Vector4 cross_product(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3) { + scalar_t a, b, c, d, e, f; /* Intermediate Values */ + Vector4 result; + + /* Calculate intermediate values. */ + a = (v2.x * v3.y) - (v2.y * v3.x); + b = (v2.x * v3.z) - (v2.z * v3.x); + c = (v2.x * v3.w) - (v2.w * v3.x); + d = (v2.y * v3.z) - (v2.z * v3.y); + e = (v2.y * v3.w) - (v2.w * v3.y); + f = (v2.z * v3.w) - (v2.w * v3.z); + + /* Calculate the result-vector components. */ + result.x = (v1.y * f) - (v1.z * e) + (v1.w * d); + result.y = - (v1.x * f) + (v1.z * c) - (v1.w * b); + result.z = (v1.x * e) - (v1.y * c) + (v1.w * a); + result.w = - (v1.x * d) + (v1.y * b) - (v1.z * a); + return result; +} + +inline Vector4 operator +(const Vector4 &v1, const Vector4 &v2) { + return Vector4(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z, v1.w + v2.w); +} + +inline Vector4 operator -(const Vector4 &v1, const Vector4 &v2) { + return Vector4(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z, v1.w - v2.w); +} + +inline Vector4 operator *(const Vector4 &v1, const Vector4 &v2) { + return Vector4(v1.x * v2.x, v1.y * v2.y, v1.z * v2.z, v1.w * v2.w); +} + +inline Vector4 operator /(const Vector4 &v1, const Vector4 &v2) { + return Vector4(v1.x / v2.x, v1.y / v2.y, v1.z / v2.z, v1.w / v2.w); +} + +inline bool operator ==(const Vector4 &v1, const Vector4 &v2) { + return (fabs(v1.x - v2.x) < XSMALL_NUMBER) && + (fabs(v1.y - v2.y) < XSMALL_NUMBER) && + (fabs(v1.z - v2.z) < XSMALL_NUMBER) && + (fabs(v1.w - v2.w) < XSMALL_NUMBER); +} + +inline void operator +=(Vector4 &v1, const Vector4 &v2) { + v1.x += v2.x; + v1.y += v2.y; + v1.z += v2.z; + v1.w += v2.w; +} + +inline void operator -=(Vector4 &v1, const Vector4 &v2) { + v1.x -= v2.x; + v1.y -= v2.y; + v1.z -= v2.z; + v1.w -= v2.w; +} + +inline void operator *=(Vector4 &v1, const Vector4 &v2) { + v1.x *= v2.x; + v1.y *= v2.y; + v1.z *= v2.z; + v1.w *= v2.w; +} + +inline void operator /=(Vector4 &v1, const Vector4 &v2) { + v1.x /= v2.x; + v1.y /= v2.y; + v1.z /= v2.z; + v1.w /= v2.w; +} + +/* binary vector (op) scalar and scalar (op) vector operations */ +inline Vector4 operator +(const Vector4 &vec, scalar_t scalar) { + return Vector4(vec.x + scalar, vec.y + scalar, vec.z + scalar, vec.w + scalar); +} + +inline Vector4 operator +(scalar_t scalar, const Vector4 &vec) { + return Vector4(vec.x + scalar, vec.y + scalar, vec.z + scalar, vec.w + scalar); +} + +inline Vector4 operator -(const Vector4 &vec, scalar_t scalar) { + return Vector4(vec.x - scalar, vec.y - scalar, vec.z - scalar, vec.w - scalar); +} + +inline Vector4 operator *(const Vector4 &vec, scalar_t scalar) { + return Vector4(vec.x * scalar, vec.y * scalar, vec.z * scalar, vec.w * scalar); +} + +inline Vector4 operator *(scalar_t scalar, const Vector4 &vec) { + return Vector4(vec.x * scalar, vec.y * scalar, vec.z * scalar, vec.w * scalar); +} + +inline Vector4 operator /(const Vector4 &vec, scalar_t scalar) { + return Vector4(vec.x / scalar, vec.y / scalar, vec.z / scalar, vec.w / scalar); +} + +inline void operator +=(Vector4 &vec, scalar_t scalar) { + vec.x += scalar; + vec.y += scalar; + vec.z += scalar; + vec.w += scalar; +} + +inline void operator -=(Vector4 &vec, scalar_t scalar) { + vec.x -= scalar; + vec.y -= scalar; + vec.z -= scalar; + vec.w -= scalar; +} + +inline void operator *=(Vector4 &vec, scalar_t scalar) { + vec.x *= scalar; + vec.y *= scalar; + vec.z *= scalar; + vec.w *= scalar; +} + +inline void operator /=(Vector4 &vec, scalar_t scalar) { + vec.x /= scalar; + vec.y /= scalar; + vec.z /= scalar; + vec.w /= scalar; +} + +inline scalar_t Vector4::length() const { + return sqrt(x*x + y*y + z*z + w*w); +} +inline scalar_t Vector4::length_sq() const { + return x*x + y*y + z*z + w*w; +} + +inline Vector4 lerp(const Vector4 &v0, const Vector4 &v1, scalar_t t) +{ + return v0 + (v1 - v0) * t; +} + +inline Vector4 catmull_rom_spline(const Vector4 &v0, const Vector4 &v1, + const Vector4 &v2, const Vector4 &v3, scalar_t t) +{ + scalar_t spline(scalar_t, scalar_t, scalar_t, scalar_t, scalar_t); + scalar_t x = spline(v0.x, v1.x, v2.x, v3.x, t); + scalar_t y = spline(v0.y, v1.y, v2.y, v3.y, t); + scalar_t z = spline(v0.z, v1.z, v2.z, v3.z, t); + scalar_t w = spline(v0.w, v1.w, v2.w, v3.w, t); + return Vector4(x, y, z, w); +} + +#endif /* __cplusplus */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vmath.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vmath.c Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,335 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include +#include +#include "vmath.h" + +/** Numerical calculation of integrals using simpson's rule */ +scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples) +{ + int i; + scalar_t h = (high - low) / (scalar_t)samples; + scalar_t sum = 0.0; + + for(i=0; i + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_H_ +#define VMATH_H_ + +#include +#include "vmath_types.h" + +#ifndef M_PI +#define M_PI PI +#endif + +#ifndef M_E +#define M_E 2.718281828459045 +#endif + +#define PI 3.141592653589793 +#define HALF_PI 1.570796326794897 +#define QUARTER_PI 0.785398163397448 +#define TWO_PI 6.283185307179586 + + +#define RAD_TO_DEG(a) ((((scalar_t)a) * 360.0) / TWO_PI) +#define DEG_TO_RAD(a) (((scalar_t)a) * (PI / 180.0)) + +#define SQ(x) ((x) * (x)) + +#define MIN(a, b) ((a) < (b) ? (a) : (b)) +#define MAX(a, b) ((a) > (b) ? (a) : (b)) + +#ifndef __GNUC__ +#define round(x) ((x) >= 0 ? (x) + 0.5 : (x) - 0.5) +#endif + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +static inline scalar_t frand(scalar_t range); +static inline vec3_t sphrand(scalar_t rad); + +scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples); +scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev); + +static inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t); + +scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t); +scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t); +scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t); + +scalar_t noise1(scalar_t x); +scalar_t noise2(scalar_t x, scalar_t y); +scalar_t noise3(scalar_t x, scalar_t y, scalar_t z); + +scalar_t fbm1(scalar_t x, int octaves); +scalar_t fbm2(scalar_t x, scalar_t y, int octaves); +scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves); + +scalar_t turbulence1(scalar_t x, int octaves); +scalar_t turbulence2(scalar_t x, scalar_t y, int octaves); +scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves); + +#ifdef __cplusplus +} +#endif /* __cplusplus */ + +#include "vmath.inl" + +#include "vector.h" +#include "matrix.h" +#include "quat.h" +#include "sphvec.h" +#include "ray.h" +#include "geom.h" + +#endif /* VMATH_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vmath.inl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vmath.inl Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,47 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#include + +/** Generates a random number in [0, range) */ +static inline scalar_t frand(scalar_t range) +{ + return range * (scalar_t)rand() / (scalar_t)RAND_MAX; +} + +/** Generates a random vector on the surface of a sphere */ +static inline vec3_t sphrand(scalar_t rad) +{ + scalar_t u = (scalar_t)rand() / RAND_MAX; + scalar_t v = (scalar_t)rand() / RAND_MAX; + + scalar_t theta = 2.0 * M_PI * u; + scalar_t phi = acos(2.0 * v - 1.0); + + vec3_t res; + res.x = rad * cos(theta) * sin(phi); + res.y = rad * sin(theta) * sin(phi); + res.z = rad * cos(phi); + return res; +} + +/** linear interpolation */ +static inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t) +{ + return a + (b - a) * t; +} diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vmath_config.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vmath_config.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,19 @@ +#ifndef VMATH_CONFIG_H_ +#define VMATH_CONFIG_H_ + +#if (__STDC_VERSION__ < 199999) +#if defined(__GNUC__) || defined(_MSC_VER) +#define inline __inline +#else +#define inline + +#ifdef VECTOR_H_ +#warning "compiling vector operations without inline, performance might suffer" +#endif /* VECTOR_H_ */ + +#endif /* gcc/msvc */ +#endif /* not C99 */ + +#define SINGLE_PRECISION_MATH + +#endif /* VMATH_CONFIG_H_ */ diff -r aebcd71cc3cf -r 96de911d05d4 prototype/vmath/vmath_types.h --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/prototype/vmath/vmath_types.h Thu Jun 28 06:05:50 2012 +0300 @@ -0,0 +1,58 @@ +/* +libvmath - a vector math library +Copyright (C) 2004-2011 John Tsiombikas + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published +by the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with this program. If not, see . +*/ + +#ifndef VMATH_TYPES_H_ +#define VMATH_TYPES_H_ + +#include "vmath_config.h" + +#define SMALL_NUMBER 1.e-4 +#define XSMALL_NUMBER 1.e-8 +#define ERROR_MARGIN 1.e-6 + + +#ifdef SINGLE_PRECISION_MATH +typedef float scalar_t; +#else +typedef double scalar_t; +#endif /* floating point precision */ + +/* vectors */ +typedef struct { scalar_t x, y; } vec2_t; +typedef struct { scalar_t x, y, z; } vec3_t; +typedef struct { scalar_t x, y, z, w; } vec4_t; + +/* quaternions */ +typedef vec4_t quat_t; + +/* matrices */ +typedef scalar_t mat3_t[3][3]; +typedef scalar_t mat4_t[4][4]; + + +#ifdef __cplusplus +class Vector2; +class Vector3; +class Vector4; +class Quaternion; +class Matrix3x3; +class Matrix4x4; +class SphVector; +#endif /* __cplusplus */ + +#endif /* VMATH_TYPES_H_ */