1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/libs/vmath/basis.h	Thu Sep 08 08:30:42 2011 +0300
     1.3 @@ -0,0 +1,39 @@
     1.4 +#ifndef VMATH_BASIS_H_
     1.5 +#define VMATH_BASIS_H_
     1.6 +
     1.7 +#include "vector.h"
     1.8 +#include "matrix.h"
     1.9 +
    1.10 +enum {
    1.11 +	LEFT_HANDED,
    1.12 +	RIGHT_HANDED
    1.13 +};
    1.14 +
    1.15 +#ifdef __cplusplus
    1.16 +extern "C" {
    1.17 +#endif	/* __cplusplus */
    1.18 +
    1.19 +void basis_matrix(mat4_t res, vec3_t i, vec3_t j, vec3_t k);
    1.20 +void basis_matrix_dir(mat4_t res, vec3_t dir);
    1.21 +
    1.22 +#ifdef __cplusplus
    1.23 +}	/* extern "C" */
    1.24 +
    1.25 +class Basis {
    1.26 +public:
    1.27 +	Vector3 i, j, k;
    1.28 +
    1.29 +	Basis();
    1.30 +	Basis(const Vector3 &i, const Vector3 &j, const Vector3 &k);
    1.31 +	Basis(const Vector3 &dir, bool left_handed = true);
    1.32 +
    1.33 +	void rotate(scalar_t x, scalar_t y, scalar_t z);
    1.34 +	void rotate(const Vector3 &axis, scalar_t angle);
    1.35 +	void rotate(const Matrix4x4 &mat);
    1.36 +	void rotate(const Quaternion &quat);
    1.37 +
    1.38 +	Matrix3x3 create_rotation_matrix() const;
    1.39 +};
    1.40 +#endif	/* __cplusplus */
    1.41 +
    1.42 +#endif	/* VMATH_BASIS_H_ */

     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/libs/vmath/basis_c.c	Thu Sep 08 08:30:42 2011 +0300
     2.3 @@ -0,0 +1,19 @@
     2.4 +#include "basis.h"
     2.5 +#include "matrix.h"
     2.6 +
     2.7 +void basis_matrix(mat4_t res, vec3_t i, vec3_t j, vec3_t k)
     2.8 +{
     2.9 +	m4_identity(res);
    2.10 +	m4_set_column(res, v4_cons(i.x, i.y, i.z, 1.0), 0);
    2.11 +	m4_set_column(res, v4_cons(j.x, j.y, j.z, 1.0), 1);
    2.12 +	m4_set_column(res, v4_cons(k.x, k.y, k.z, 1.0), 2);
    2.13 +}
    2.14 +
    2.15 +void basis_matrix_dir(mat4_t res, vec3_t dir)
    2.16 +{
    2.17 +	vec3_t k = v3_normalize(dir);
    2.18 +	vec3_t j = {0, 1, 0};
    2.19 +	vec3_t i = v3_cross(j, k);
    2.20 +	j = v3_cross(k, i);
    2.21 +	basis_matrix(res, i, j, k);
    2.22 +}

     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/libs/vmath/geom.c	Thu Sep 08 08:30:42 2011 +0300
     3.3 @@ -0,0 +1,131 @@
     3.4 +#include <math.h>
     3.5 +#include "geom.h"
     3.6 +#include "vector.h"
     3.7 +
     3.8 +plane_t plane_cons(scalar_t nx, scalar_t ny, scalar_t nz, scalar_t d)
     3.9 +{
    3.10 +	plane_t p;
    3.11 +	p.norm.x = nx;
    3.12 +	p.norm.y = ny;
    3.13 +	p.norm.z = nz;
    3.14 +	p.d = d;
    3.15 +	return p;
    3.16 +}
    3.17 +
    3.18 +plane_t plane_poly(vec3_t v0, vec3_t v1, vec3_t v2)
    3.19 +{
    3.20 +	vec3_t a, b, norm;
    3.21 +
    3.22 +	a = v3_sub(v1, v0);
    3.23 +	b = v3_sub(v2, v0);
    3.24 +	norm = v3_cross(a, b);
    3.25 +	norm = v3_normalize(norm);
    3.26 +
    3.27 +	return plane_ptnorm(v0, norm);
    3.28 +}
    3.29 +
    3.30 +plane_t plane_ptnorm(vec3_t pt, vec3_t normal)
    3.31 +{
    3.32 +	plane_t plane;
    3.33 +
    3.34 +	plane.norm = normal;
    3.35 +	plane.d = v3_dot(pt, normal);
    3.36 +
    3.37 +	return plane;
    3.38 +}
    3.39 +
    3.40 +plane_t plane_invert(plane_t p)
    3.41 +{
    3.42 +	p.norm = v3_neg(p.norm);
    3.43 +	p.d = -p.d;
    3.44 +	return p;
    3.45 +}
    3.46 +
    3.47 +scalar_t plane_signed_dist(plane_t plane, vec3_t pt)
    3.48 +{
    3.49 +	vec3_t pp = plane_point(plane);
    3.50 +	vec3_t pptopt = v3_sub(pt, pp);
    3.51 +	return v3_dot(pptopt, plane.norm);
    3.52 +}
    3.53 +
    3.54 +scalar_t plane_dist(plane_t plane, vec3_t pt)
    3.55 +{
    3.56 +	return fabs(plane_signed_dist(plane, pt));
    3.57 +}
    3.58 +
    3.59 +vec3_t plane_point(plane_t plane)
    3.60 +{
    3.61 +	return v3_scale(plane.norm, plane.d);
    3.62 +}
    3.63 +
    3.64 +int plane_ray_intersect(ray_t ray, plane_t plane, scalar_t *pos)
    3.65 +{
    3.66 +	vec3_t pt, orig_to_pt;
    3.67 +	scalar_t ndotdir;
    3.68 +
    3.69 +	pt = plane_point(plane);
    3.70 +	ndotdir = v3_dot(plane.norm, ray.dir);
    3.71 +
    3.72 +	if(fabs(ndotdir) < 1e-7) {
    3.73 +		return 0;
    3.74 +	}
    3.75 +
    3.76 +	if(pos) {
    3.77 +		orig_to_pt = v3_sub(pt, ray.origin);
    3.78 +		*pos = v3_dot(plane.norm, orig_to_pt) / ndotdir;
    3.79 +	}
    3.80 +	return 1;
    3.81 +}
    3.82 +
    3.83 +sphere_t sphere_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t rad)
    3.84 +{
    3.85 +	sphere_t sph;
    3.86 +	sph.pos.x = x;
    3.87 +	sph.pos.y = y;
    3.88 +	sph.pos.z = z;
    3.89 +	sph.rad = rad;
    3.90 +	return sph;
    3.91 +}
    3.92 +
    3.93 +int sphere_ray_intersect(ray_t ray, sphere_t sph, scalar_t *pos)
    3.94 +{
    3.95 +	scalar_t a, b, c, d, sqrt_d, t1, t2, t;
    3.96 +
    3.97 +	a = v3_dot(ray.dir, ray.dir);
    3.98 +	b = 2.0 * ray.dir.x * (ray.origin.x - sph.pos.x) +
    3.99 +		2.0 * ray.dir.y * (ray.origin.y - sph.pos.y) +
   3.100 +		2.0 * ray.dir.z * (ray.origin.z - sph.pos.z);
   3.101 +	c = v3_dot(sph.pos, sph.pos) + v3_dot(ray.origin, ray.origin) +
   3.102 +		2.0 * v3_dot(v3_neg(sph.pos), ray.origin) - sph.rad * sph.rad;
   3.103 +
   3.104 +	d = b * b - 4.0 * a * c;
   3.105 +	if(d < 0.0) {
   3.106 +		return 0;
   3.107 +	}
   3.108 +
   3.109 +	sqrt_d = sqrt(d);
   3.110 +	t1 = (-b + sqrt_d) / (2.0 * a);
   3.111 +	t2 = (-b - sqrt_d) / (2.0 * a);
   3.112 +
   3.113 +	if(t1 < 1e-7 || t1 > 1.0) {
   3.114 +		t1 = t2;
   3.115 +	}
   3.116 +	if(t2 < 1e-7 || t2 > 1.0) {
   3.117 +		t2 = t1;
   3.118 +	}
   3.119 +	t = t1 < t2 ? t1 : t2;
   3.120 +
   3.121 +	if(t < 1e-7 || t > 1.0) {
   3.122 +		return 0;
   3.123 +	}
   3.124 +
   3.125 +	if(pos) {
   3.126 +		*pos = t;
   3.127 +	}
   3.128 +	return 1;
   3.129 +}
   3.130 +
   3.131 +int sphere_sphere_intersect(sphere_t sph1, sphere_t sph2, scalar_t *pos, scalar_t *rad)
   3.132 +{
   3.133 +	return -1;
   3.134 +}

     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/libs/vmath/geom.h	Thu Sep 08 08:30:42 2011 +0300
     4.3 @@ -0,0 +1,55 @@
     4.4 +#ifndef GEOM_H_
     4.5 +#define GEOM_H_
     4.6 +
     4.7 +#include "vector.h"
     4.8 +#include "ray.h"
     4.9 +
    4.10 +typedef struct {
    4.11 +	vec3_t norm;
    4.12 +	scalar_t d;
    4.13 +} plane_t;
    4.14 +
    4.15 +typedef struct {
    4.16 +	vec3_t pos;
    4.17 +	scalar_t rad;
    4.18 +} sphere_t;
    4.19 +
    4.20 +typedef struct {
    4.21 +	vec3_t min, max;
    4.22 +} aabox_t;
    4.23 +
    4.24 +#ifdef __cplusplus
    4.25 +extern "C" {
    4.26 +#endif
    4.27 +
    4.28 +/* planes are good... you need planes, yes you do */
    4.29 +plane_t plane_cons(scalar_t nx, scalar_t ny, scalar_t nz, scalar_t d);
    4.30 +plane_t plane_poly(vec3_t v0, vec3_t v1, vec3_t v2);
    4.31 +plane_t plane_ptnorm(vec3_t pt, vec3_t normal);
    4.32 +
    4.33 +plane_t plane_invert(plane_t p);
    4.34 +
    4.35 +scalar_t plane_signed_dist(plane_t plane, vec3_t pt);
    4.36 +scalar_t plane_dist(plane_t plane, vec3_t pt);
    4.37 +vec3_t plane_point(plane_t plane);
    4.38 +
    4.39 +int plane_ray_intersect(ray_t ray, plane_t plane, scalar_t *pos);
    4.40 +
    4.41 +/* spheres always come in handy */
    4.42 +sphere_t sphere_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t rad);
    4.43 +
    4.44 +int sphere_ray_intersect(ray_t ray, sphere_t sph, scalar_t *pos);
    4.45 +int sphere_sphere_intersect(sphere_t sph1, sphere_t sph2, scalar_t *pos, scalar_t *rad);
    4.46 +
    4.47 +#ifdef __cplusplus
    4.48 +}
    4.49 +
    4.50 +/* TODO
    4.51 +class Plane : public plane_t {
    4.52 +public:
    4.53 +};
    4.54 +*/
    4.55 +
    4.56 +#endif
    4.57 +
    4.58 +#endif	/* GEOM_H_ */

     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/libs/vmath/matrix.h	Thu Sep 08 08:30:42 2011 +0300
     5.3 @@ -0,0 +1,187 @@
     5.4 +#ifndef VMATH_MATRIX_H_
     5.5 +#define VMATH_MATRIX_H_
     5.6 +
     5.7 +#include <stdio.h>
     5.8 +#include "vmath_types.h"
     5.9 +
    5.10 +#ifdef __cplusplus
    5.11 +extern "C" {
    5.12 +#endif	/* __cplusplus */
    5.13 +
    5.14 +/* C matrix 3x3 functions */
    5.15 +static inline void m3_identity(mat3_t m);
    5.16 +static inline void m3_cons(mat3_t m,
    5.17 +		scalar_t m11, scalar_t m12, scalar_t m13,
    5.18 +		scalar_t m21, scalar_t m22, scalar_t m23,
    5.19 +		scalar_t m31, scalar_t m32, scalar_t m33);
    5.20 +static inline void m3_copy(mat3_t dest, mat3_t src);
    5.21 +void m3_to_m4(mat4_t dest, mat3_t src);
    5.22 +
    5.23 +void m3_print(FILE *fp, mat3_t m);
    5.24 +
    5.25 +/* C matrix 4x4 functions */
    5.26 +static inline void m4_identity(mat4_t m);
    5.27 +static inline void m4_cons(mat4_t m,
    5.28 +		scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14,
    5.29 +		scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24,
    5.30 +		scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34,
    5.31 +		scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44);
    5.32 +static inline void m4_copy(mat4_t dest, mat4_t src);
    5.33 +void m4_to_m3(mat3_t dest, mat4_t src);
    5.34 +
    5.35 +static inline void m4_mult(mat4_t res, mat4_t m1, mat4_t m2);
    5.36 +
    5.37 +void m4_translate(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
    5.38 +void m4_rotate(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
    5.39 +void m4_rotate_x(mat4_t m, scalar_t angle);
    5.40 +void m4_rotate_y(mat4_t m, scalar_t angle);
    5.41 +void m4_rotate_z(mat4_t m, scalar_t angle);
    5.42 +void m4_rotate_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
    5.43 +void m4_rotate_quat(mat4_t m, quat_t q);
    5.44 +void m4_scale(mat4_t m, scalar_t x, scalar_t y, scalar_t z);
    5.45 +static inline void m4_set_column(mat4_t m, vec4_t v, int idx);
    5.46 +static inline void m4_set_row(mat4_t m, vec4_t v, int idx);
    5.47 +
    5.48 +void m4_transpose(mat4_t res, mat4_t m);
    5.49 +scalar_t m4_determinant(mat4_t m);
    5.50 +void m4_adjoint(mat4_t res, mat4_t m);
    5.51 +void m4_inverse(mat4_t res, mat4_t m);
    5.52 +
    5.53 +void m4_print(FILE *fp, mat4_t m);
    5.54 +
    5.55 +#ifdef __cplusplus
    5.56 +}
    5.57 +
    5.58 +/* when included from C++ source files, also define the matrix classes */
    5.59 +#include <iostream>
    5.60 +
    5.61 +/** 3x3 matrix */
    5.62 +class Matrix3x3 {
    5.63 +private:
    5.64 +	scalar_t m[3][3];
    5.65 +
    5.66 +public:
    5.67 +
    5.68 +	static Matrix3x3 identity;
    5.69 +
    5.70 +	Matrix3x3();
    5.71 +	Matrix3x3(	scalar_t m11, scalar_t m12, scalar_t m13,
    5.72 +				scalar_t m21, scalar_t m22, scalar_t m23,
    5.73 +				scalar_t m31, scalar_t m32, scalar_t m33);
    5.74 +	Matrix3x3(const mat3_t cmat);
    5.75 +
    5.76 +	Matrix3x3(const Matrix4x4 &mat4x4);
    5.77 +
    5.78 +	/* binary operations matrix (op) matrix */
    5.79 +	friend Matrix3x3 operator +(const Matrix3x3 &m1, const Matrix3x3 &m2);
    5.80 +	friend Matrix3x3 operator -(const Matrix3x3 &m1, const Matrix3x3 &m2);
    5.81 +	friend Matrix3x3 operator *(const Matrix3x3 &m1, const Matrix3x3 &m2);
    5.82 +
    5.83 +	friend void operator +=(Matrix3x3 &m1, const Matrix3x3 &m2);
    5.84 +	friend void operator -=(Matrix3x3 &m1, const Matrix3x3 &m2);
    5.85 +	friend void operator *=(Matrix3x3 &m1, const Matrix3x3 &m2);
    5.86 +
    5.87 +	/* binary operations matrix (op) scalar and scalar (op) matrix */
    5.88 +	friend Matrix3x3 operator *(const Matrix3x3 &mat, scalar_t scalar);
    5.89 +	friend Matrix3x3 operator *(scalar_t scalar, const Matrix3x3 &mat);
    5.90 +
    5.91 +	friend void operator *=(Matrix3x3 &mat, scalar_t scalar);
    5.92 +
    5.93 +	inline scalar_t *operator [](int index);
    5.94 +	inline const scalar_t *operator [](int index) const;
    5.95 +
    5.96 +	inline void reset_identity();
    5.97 +
    5.98 +	void translate(const Vector2 &trans);
    5.99 +	void set_translation(const Vector2 &trans);
   5.100 +
   5.101 +	void rotate(scalar_t angle);						/* 2d rotation */
   5.102 +	void rotate(const Vector3 &euler_angles);			/* 3d rotation with euler angles */
   5.103 +	void rotate(const Vector3 &axis, scalar_t angle);	/* 3d axis/angle rotation */
   5.104 +	void set_rotation(scalar_t angle);
   5.105 +	void set_rotation(const Vector3 &euler_angles);
   5.106 +	void set_rotation(const Vector3 &axis, scalar_t angle);
   5.107 +
   5.108 +	void scale(const Vector3 &scale_vec);
   5.109 +	void set_scaling(const Vector3 &scale_vec);
   5.110 +
   5.111 +	void set_column_vector(const Vector3 &vec, unsigned int col_index);
   5.112 +	void set_row_vector(const Vector3 &vec, unsigned int row_index);
   5.113 +	Vector3 get_column_vector(unsigned int col_index) const;
   5.114 +	Vector3 get_row_vector(unsigned int row_index) const;
   5.115 +
   5.116 +	void transpose();
   5.117 +	Matrix3x3 transposed() const;
   5.118 +	scalar_t determinant() const;
   5.119 +	Matrix3x3 inverse() const;
   5.120 +
   5.121 +	friend std::ostream &operator <<(std::ostream &out, const Matrix3x3 &mat);
   5.122 +};
   5.123 +
   5.124 +/** 4x4 matrix */
   5.125 +class Matrix4x4 {
   5.126 +private:
   5.127 +	scalar_t m[4][4];
   5.128 +
   5.129 +public:
   5.130 +
   5.131 +	static Matrix4x4 identity;
   5.132 +
   5.133 +	Matrix4x4();
   5.134 +	Matrix4x4(	scalar_t m11, scalar_t m12, scalar_t m13, scalar_t m14,
   5.135 +				scalar_t m21, scalar_t m22, scalar_t m23, scalar_t m24,
   5.136 +				scalar_t m31, scalar_t m32, scalar_t m33, scalar_t m34,
   5.137 +				scalar_t m41, scalar_t m42, scalar_t m43, scalar_t m44);
   5.138 +	Matrix4x4(const mat4_t cmat);
   5.139 +
   5.140 +	Matrix4x4(const Matrix3x3 &mat3x3);
   5.141 +
   5.142 +	/* binary operations matrix (op) matrix */
   5.143 +	friend Matrix4x4 operator +(const Matrix4x4 &m1, const Matrix4x4 &m2);
   5.144 +	friend Matrix4x4 operator -(const Matrix4x4 &m1, const Matrix4x4 &m2);
   5.145 +	friend Matrix4x4 operator *(const Matrix4x4 &m1, const Matrix4x4 &m2);
   5.146 +
   5.147 +	friend void operator +=(Matrix4x4 &m1, const Matrix4x4 &m2);
   5.148 +	friend void operator -=(Matrix4x4 &m1, const Matrix4x4 &m2);
   5.149 +	friend inline void operator *=(Matrix4x4 &m1, const Matrix4x4 &m2);
   5.150 +
   5.151 +	/* binary operations matrix (op) scalar and scalar (op) matrix */
   5.152 +	friend Matrix4x4 operator *(const Matrix4x4 &mat, scalar_t scalar);
   5.153 +	friend Matrix4x4 operator *(scalar_t scalar, const Matrix4x4 &mat);
   5.154 +
   5.155 +	friend void operator *=(Matrix4x4 &mat, scalar_t scalar);
   5.156 +
   5.157 +	inline scalar_t *operator [](int index);
   5.158 +	inline const scalar_t *operator [](int index) const;
   5.159 +
   5.160 +	inline void reset_identity();
   5.161 +
   5.162 +	void translate(const Vector3 &trans);
   5.163 +	void set_translation(const Vector3 &trans);
   5.164 +
   5.165 +	void rotate(const Vector3 &euler_angles);			/* 3d rotation with euler angles */
   5.166 +	void rotate(const Vector3 &axis, scalar_t angle);	/* 3d axis/angle rotation */
   5.167 +	void set_rotation(const Vector3 &euler_angles);
   5.168 +	void set_rotation(const Vector3 &axis, scalar_t angle);
   5.169 +
   5.170 +	void scale(const Vector4 &scale_vec);
   5.171 +	void set_scaling(const Vector4 &scale_vec);
   5.172 +
   5.173 +	void set_column_vector(const Vector4 &vec, unsigned int col_index);
   5.174 +	void set_row_vector(const Vector4 &vec, unsigned int row_index);
   5.175 +	Vector4 get_column_vector(unsigned int col_index) const;
   5.176 +	Vector4 get_row_vector(unsigned int row_index) const;
   5.177 +
   5.178 +	void transpose();
   5.179 +	Matrix4x4 transposed() const;
   5.180 +	scalar_t determinant() const;
   5.181 +	Matrix4x4 adjoint() const;
   5.182 +	Matrix4x4 inverse() const;
   5.183 +
   5.184 +	friend std::ostream &operator <<(std::ostream &out, const Matrix4x4 &mat);
   5.185 +};
   5.186 +#endif	/* __cplusplus */
   5.187 +
   5.188 +#include "matrix.inl"
   5.189 +
   5.190 +#endif	/* VMATH_MATRIX_H_ */

     6.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     6.2 +++ b/libs/vmath/matrix_c.c	Thu Sep 08 08:30:42 2011 +0300
     6.3 @@ -0,0 +1,246 @@
     6.4 +#include <stdio.h>
     6.5 +#include "matrix.h"
     6.6 +#include "vector.h"
     6.7 +#include "quat.h"
     6.8 +
     6.9 +void m3_to_m4(mat4_t dest, mat3_t src)
    6.10 +{
    6.11 +	int i, j;
    6.12 +
    6.13 +	memset(dest, 0, sizeof(mat4_t));
    6.14 +	for(i=0; i<3; i++) {
    6.15 +		for(j=0; j<3; j++) {
    6.16 +			dest[i][j] = src[i][j];
    6.17 +		}
    6.18 +	}
    6.19 +	dest[3][3] = 1.0;
    6.20 +}
    6.21 +
    6.22 +void m3_print(FILE *fp, mat3_t m)
    6.23 +{
    6.24 +	int i;
    6.25 +	for(i=0; i<3; i++) {
    6.26 +		fprintf(fp, "[ %12.5f %12.5f %12.5f ]\n", (float)m[i][0], (float)m[i][1], (float)m[i][2]);
    6.27 +	}
    6.28 +}
    6.29 +
    6.30 +/* C matrix 4x4 functions */
    6.31 +void m4_to_m3(mat3_t dest, mat4_t src)
    6.32 +{
    6.33 +	int i, j;
    6.34 +	for(i=0; i<3; i++) {
    6.35 +		for(j=0; j<3; j++) {
    6.36 +			dest[i][j] = src[i][j];
    6.37 +		}
    6.38 +	}
    6.39 +}
    6.40 +
    6.41 +void m4_translate(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
    6.42 +{
    6.43 +	mat4_t tm;
    6.44 +	m4_identity(tm);
    6.45 +	tm[0][3] = x;
    6.46 +	tm[1][3] = y;
    6.47 +	tm[2][3] = z;
    6.48 +	m4_mult(m, m, tm);
    6.49 +}
    6.50 +
    6.51 +void m4_rotate(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
    6.52 +{
    6.53 +	m4_rotate_x(m, x);
    6.54 +	m4_rotate_y(m, y);
    6.55 +	m4_rotate_z(m, z);
    6.56 +}
    6.57 +
    6.58 +void m4_rotate_x(mat4_t m, scalar_t angle)
    6.59 +{
    6.60 +	mat4_t rm;
    6.61 +	m4_identity(rm);
    6.62 +	rm[1][1] = cos(angle); rm[1][2] = -sin(angle);
    6.63 +	rm[2][1] = sin(angle); rm[2][2] = cos(angle);
    6.64 +	m4_mult(m, m, rm);
    6.65 +}
    6.66 +
    6.67 +void m4_rotate_y(mat4_t m, scalar_t angle)
    6.68 +{
    6.69 +	mat4_t rm;
    6.70 +	m4_identity(rm);
    6.71 +	rm[0][0] = cos(angle); rm[0][2] = sin(angle);
    6.72 +	rm[2][0] = -sin(angle); rm[2][2] = cos(angle);
    6.73 +	m4_mult(m, m, rm);
    6.74 +}
    6.75 +
    6.76 +void m4_rotate_z(mat4_t m, scalar_t angle)
    6.77 +{
    6.78 +	mat4_t rm;
    6.79 +	m4_identity(rm);
    6.80 +	rm[0][0] = cos(angle); rm[0][1] = -sin(angle);
    6.81 +	rm[1][0] = sin(angle); rm[1][1] = cos(angle);
    6.82 +	m4_mult(m, m, rm);
    6.83 +}
    6.84 +
    6.85 +void m4_rotate_axis(mat4_t m, scalar_t angle, scalar_t x, scalar_t y, scalar_t z)
    6.86 +{
    6.87 +	mat4_t xform;
    6.88 +	scalar_t sina = sin(angle);
    6.89 +	scalar_t cosa = cos(angle);
    6.90 +	scalar_t one_minus_cosa = 1.0 - cosa;
    6.91 +	scalar_t nxsq = x * x;
    6.92 +	scalar_t nysq = y * y;
    6.93 +	scalar_t nzsq = z * z;
    6.94 +
    6.95 +	m4_identity(xform);
    6.96 +	xform[0][0] = nxsq + (1.0 - nxsq) * cosa;
    6.97 +	xform[0][1] = x * y * one_minus_cosa - z * sina;
    6.98 +	xform[0][2] = x * z * one_minus_cosa + y * sina;
    6.99 +	xform[1][0] = x * y * one_minus_cosa + z * sina;
   6.100 +	xform[1][1] = nysq + (1.0 - nysq) * cosa;
   6.101 +	xform[1][2] = y * z * one_minus_cosa - x * sina;
   6.102 +	xform[2][0] = x * z * one_minus_cosa - y * sina;
   6.103 +	xform[2][1] = y * z * one_minus_cosa + x * sina;
   6.104 +	xform[2][2] = nzsq + (1.0 - nzsq) * cosa;
   6.105 +
   6.106 +	m4_mult(m, m, xform);
   6.107 +}
   6.108 +
   6.109 +void m4_rotate_quat(mat4_t m, quat_t q)
   6.110 +{
   6.111 +	mat4_t rm;
   6.112 +	quat_to_mat4(rm, q);
   6.113 +	m4_mult(m, m, rm);
   6.114 +}
   6.115 +
   6.116 +void m4_scale(mat4_t m, scalar_t x, scalar_t y, scalar_t z)
   6.117 +{
   6.118 +	mat4_t sm;
   6.119 +	m4_identity(sm);
   6.120 +	sm[0][0] = x;
   6.121 +	sm[1][1] = y;
   6.122 +	sm[2][2] = z;
   6.123 +	m4_mult(m, m, sm);
   6.124 +}
   6.125 +
   6.126 +void m4_transpose(mat4_t res, mat4_t m)
   6.127 +{
   6.128 +	int i, j;
   6.129 +	mat4_t tmp;
   6.130 +	m4_copy(tmp, m);
   6.131 +
   6.132 +	for(i=0; i<4; i++) {
   6.133 +		for(j=0; j<4; j++) {
   6.134 +			res[i][j] = tmp[j][i];
   6.135 +		}
   6.136 +	}
   6.137 +}
   6.138 +
   6.139 +scalar_t m4_determinant(mat4_t m)
   6.140 +{
   6.141 +	scalar_t det11 =	(m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
   6.142 +						(m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
   6.143 +						(m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
   6.144 +
   6.145 +	scalar_t det12 =	(m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
   6.146 +						(m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
   6.147 +						(m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
   6.148 +
   6.149 +	scalar_t det13 =	(m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
   6.150 +						(m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
   6.151 +						(m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
   6.152 +
   6.153 +	scalar_t det14 =	(m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
   6.154 +						(m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
   6.155 +						(m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
   6.156 +
   6.157 +	return m[0][0] * det11 - m[0][1] * det12 + m[0][2] * det13 - m[0][3] * det14;
   6.158 +}
   6.159 +
   6.160 +void m4_adjoint(mat4_t res, mat4_t m)
   6.161 +{
   6.162 +	int i, j;
   6.163 +	mat4_t coef;
   6.164 +
   6.165 +	coef[0][0] =	(m[1][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
   6.166 +					(m[1][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
   6.167 +					(m[1][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
   6.168 +	coef[0][1] =	(m[1][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
   6.169 +					(m[1][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
   6.170 +					(m[1][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
   6.171 +	coef[0][2] =	(m[1][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
   6.172 +					(m[1][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
   6.173 +					(m[1][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
   6.174 +	coef[0][3] =	(m[1][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
   6.175 +					(m[1][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
   6.176 +					(m[1][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
   6.177 +
   6.178 +	coef[1][0] =	(m[0][1] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
   6.179 +					(m[0][2] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) +
   6.180 +					(m[0][3] * (m[2][1] * m[3][2] - m[3][1] * m[2][2]));
   6.181 +	coef[1][1] =	(m[0][0] * (m[2][2] * m[3][3] - m[3][2] * m[2][3])) -
   6.182 +					(m[0][2] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
   6.183 +					(m[0][3] * (m[2][0] * m[3][2] - m[3][0] * m[2][2]));
   6.184 +	coef[1][2] =	(m[0][0] * (m[2][1] * m[3][3] - m[3][1] * m[2][3])) -
   6.185 +					(m[0][1] * (m[2][0] * m[3][3] - m[3][0] * m[2][3])) +
   6.186 +					(m[0][3] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
   6.187 +	coef[1][3] =	(m[0][0] * (m[2][1] * m[3][2] - m[3][1] * m[2][2])) -
   6.188 +					(m[0][1] * (m[2][0] * m[3][2] - m[3][0] * m[2][2])) +
   6.189 +					(m[0][2] * (m[2][0] * m[3][1] - m[3][0] * m[2][1]));
   6.190 +
   6.191 +	coef[2][0] =	(m[0][1] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) -
   6.192 +					(m[0][2] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) +
   6.193 +					(m[0][3] * (m[1][1] * m[3][2] - m[3][1] * m[1][2]));
   6.194 +	coef[2][1] =	(m[0][0] * (m[1][2] * m[3][3] - m[3][2] * m[1][3])) -
   6.195 +					(m[0][2] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) +
   6.196 +					(m[0][3] * (m[1][0] * m[3][2] - m[3][0] * m[1][2]));
   6.197 +	coef[2][2] =	(m[0][0] * (m[1][1] * m[3][3] - m[3][1] * m[1][3])) -
   6.198 +					(m[0][1] * (m[1][0] * m[3][3] - m[3][0] * m[1][3])) +
   6.199 +					(m[0][3] * (m[1][0] * m[3][1] - m[3][0] * m[1][1]));
   6.200 +	coef[2][3] =	(m[0][0] * (m[1][1] * m[3][2] - m[3][1] * m[1][2])) -
   6.201 +					(m[0][1] * (m[1][0] * m[3][2] - m[3][0] * m[1][2])) +
   6.202 +					(m[0][2] * (m[1][0] * m[3][1] - m[3][0] * m[1][1]));
   6.203 +
   6.204 +	coef[3][0] =	(m[0][1] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) -
   6.205 +					(m[0][2] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) +
   6.206 +					(m[0][3] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]));
   6.207 +	coef[3][1] =	(m[0][0] * (m[1][2] * m[2][3] - m[2][2] * m[1][3])) -
   6.208 +					(m[0][2] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) +
   6.209 +					(m[0][3] * (m[1][0] * m[2][2] - m[2][0] * m[1][2]));
   6.210 +	coef[3][2] =	(m[0][0] * (m[1][1] * m[2][3] - m[2][1] * m[1][3])) -
   6.211 +					(m[0][1] * (m[1][0] * m[2][3] - m[2][0] * m[1][3])) +
   6.212 +					(m[0][3] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]));
   6.213 +	coef[3][3] =	(m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2])) -
   6.214 +					(m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2])) +
   6.215 +					(m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]));
   6.216 +
   6.217 +	m4_transpose(res, coef);
   6.218 +
   6.219 +	for(i=0; i<4; i++) {
   6.220 +		for(j=0; j<4; j++) {
   6.221 +			res[i][j] = j % 2 ? -res[i][j] : res[i][j];
   6.222 +			if(i % 2) res[i][j] = -res[i][j];
   6.223 +		}
   6.224 +	}
   6.225 +}
   6.226 +
   6.227 +void m4_inverse(mat4_t res, mat4_t m)
   6.228 +{
   6.229 +	int i, j;
   6.230 +	mat4_t adj;
   6.231 +	scalar_t det;
   6.232 +
   6.233 +	m4_adjoint(adj, m);
   6.234 +	det = m4_determinant(m);
   6.235 +	
   6.236 +	for(i=0; i<4; i++) {
   6.237 +		for(j=0; j<4; j++) {
   6.238 +			res[i][j] = adj[i][j] / det;
   6.239 +		}
   6.240 +	}
   6.241 +}
   6.242 +
   6.243 +void m4_print(FILE *fp, mat4_t m)
   6.244 +{
   6.245 +	int i;
   6.246 +	for(i=0; i<4; i++) {
   6.247 +		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]);
   6.248 +	}
   6.249 +}

     7.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     7.2 +++ b/libs/vmath/quat.h	Thu Sep 08 08:30:42 2011 +0300
     7.3 @@ -0,0 +1,101 @@
     7.4 +#ifndef VMATH_QUATERNION_H_
     7.5 +#define VMATH_QUATERNION_H_
     7.6 +
     7.7 +#include <stdio.h>
     7.8 +#include "vmath_types.h"
     7.9 +#include "vector.h"
    7.10 +
    7.11 +#ifdef __cplusplus
    7.12 +extern "C" {
    7.13 +#endif	/* __cplusplus */
    7.14 +
    7.15 +#define quat_cons(s, x, y, z)	v4_cons(x, y, z, s)
    7.16 +#define quat_vec(q)				v3_cons((q).x, (q).y, (q).z)
    7.17 +#define quat_s(q)				((q).w)
    7.18 +#define quat_identity()			quat_cons(1.0, 0.0, 0.0, 0.0)
    7.19 +void quat_print(FILE *fp, quat_t q);
    7.20 +
    7.21 +#define quat_add		v4_add
    7.22 +#define quat_sub		v4_sub
    7.23 +#define quat_neg		v4_neg
    7.24 +
    7.25 +static inline quat_t quat_mul(quat_t q1, quat_t q2);
    7.26 +
    7.27 +static inline quat_t quat_conjugate(quat_t q);
    7.28 +
    7.29 +#define quat_length		v4_length
    7.30 +#define quat_length_sq	v4_length_sq
    7.31 +
    7.32 +#define quat_normalize	v4_normalize
    7.33 +static inline quat_t quat_inverse(quat_t q);
    7.34 +
    7.35 +quat_t quat_rotate(quat_t q, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
    7.36 +quat_t quat_rotate_quat(quat_t q, quat_t rotq);
    7.37 +
    7.38 +static inline void quat_to_mat3(mat3_t res, quat_t q);
    7.39 +static inline void quat_to_mat4(mat4_t res, quat_t q);
    7.40 +
    7.41 +#define quat_lerp quat_slerp
    7.42 +quat_t quat_slerp(quat_t q1, quat_t q2, scalar_t t);
    7.43 +
    7.44 +
    7.45 +#ifdef __cplusplus
    7.46 +}	/* extern "C" */
    7.47 +
    7.48 +#include <iostream>
    7.49 +
    7.50 +/* Quaternion */
    7.51 +class Quaternion {
    7.52 +public:
    7.53 +	scalar_t s;
    7.54 +	Vector3 v;
    7.55 +
    7.56 +	Quaternion();
    7.57 +	Quaternion(scalar_t s, const Vector3 &v);
    7.58 +	Quaternion(scalar_t s, scalar_t x, scalar_t y, scalar_t z);
    7.59 +	Quaternion(const Vector3 &axis, scalar_t angle);
    7.60 +	Quaternion(const quat_t &quat);
    7.61 +
    7.62 +	Quaternion operator +(const Quaternion &quat) const;
    7.63 +	Quaternion operator -(const Quaternion &quat) const;
    7.64 +	Quaternion operator -() const;
    7.65 +	Quaternion operator *(const Quaternion &quat) const;
    7.66 +
    7.67 +	void operator +=(const Quaternion &quat);
    7.68 +	void operator -=(const Quaternion &quat);
    7.69 +	void operator *=(const Quaternion &quat);
    7.70 +
    7.71 +	void reset_identity();
    7.72 +
    7.73 +	Quaternion conjugate() const;
    7.74 +
    7.75 +	scalar_t length() const;
    7.76 +	scalar_t length_sq() const;
    7.77 +
    7.78 +	void normalize();
    7.79 +	Quaternion normalized() const;
    7.80 +
    7.81 +	Quaternion inverse() const;
    7.82 +
    7.83 +	void set_rotation(const Vector3 &axis, scalar_t angle);
    7.84 +	void rotate(const Vector3 &axis, scalar_t angle);
    7.85 +	/* note: this is a totally different operation from the above
    7.86 +	 * this treats the quaternion as signifying direction and rotates
    7.87 +	 * it by a rotation quaternion by rot * q * rot'
    7.88 +	 */
    7.89 +	void rotate(const Quaternion &q);
    7.90 +
    7.91 +	Matrix3x3 get_rotation_matrix() const;
    7.92 +
    7.93 +	friend Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t);
    7.94 +
    7.95 +	friend std::ostream &operator <<(std::ostream &out, const Quaternion &q);
    7.96 +};
    7.97 +
    7.98 +Quaternion slerp(const Quaternion &q1, const Quaternion &q2, scalar_t t);
    7.99 +inline Quaternion lerp(const Quaternion &q1, const Quaternion &q2, scalar_t t);
   7.100 +#endif	/* __cplusplus */
   7.101 +
   7.102 +#include "quat.inl"
   7.103 +
   7.104 +#endif	/* VMATH_QUATERNION_H_ */

     8.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     8.2 +++ b/libs/vmath/quat_c.c	Thu Sep 08 08:30:42 2011 +0300
     8.3 @@ -0,0 +1,42 @@
     8.4 +#include <stdio.h>
     8.5 +#include <math.h>
     8.6 +#include "quat.h"
     8.7 +
     8.8 +void quat_print(FILE *fp, quat_t q)
     8.9 +{
    8.10 +	fprintf(fp, "([ %.4f %.4f %.4f ] %.4f)", q.x, q.y, q.z, q.w);
    8.11 +}
    8.12 +
    8.13 +quat_t quat_rotate(quat_t q, scalar_t angle, scalar_t x, scalar_t y, scalar_t z)
    8.14 +{
    8.15 +	quat_t rq;
    8.16 +	scalar_t half_angle = angle * 0.5;
    8.17 +	scalar_t sin_half = sin(half_angle);
    8.18 +
    8.19 +	rq.w = cos(half_angle);
    8.20 +	rq.x = x * sin_half;
    8.21 +	rq.y = y * sin_half;
    8.22 +	rq.z = z * sin_half;
    8.23 +
    8.24 +	return quat_mul(q, rq);
    8.25 +}
    8.26 +
    8.27 +quat_t quat_rotate_quat(quat_t q, quat_t rotq)
    8.28 +{
    8.29 +	return quat_mul(quat_mul(rotq, q), quat_conjugate(rotq));
    8.30 +}
    8.31 +
    8.32 +quat_t quat_slerp(quat_t q1, quat_t q2, scalar_t t)
    8.33 +{
    8.34 +	quat_t res;
    8.35 +	scalar_t angle = acos(q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z);
    8.36 +	scalar_t a = sin((1.0f - t) * angle);
    8.37 +	scalar_t b = sin(t * angle);
    8.38 +	scalar_t c = sin(angle);
    8.39 +
    8.40 +	res.x = (q1.x * a + q2.x * b) / c;
    8.41 +	res.y = (q1.y * a + q2.y * b) / c;
    8.42 +	res.z = (q1.z * a + q2.z * b) / c;
    8.43 +	res.w = (q1.w * a + q2.w * b) / c;
    8.44 +	return quat_normalize(res);
    8.45 +}

     9.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     9.2 +++ b/libs/vmath/ray.h	Thu Sep 08 08:30:42 2011 +0300
     9.3 @@ -0,0 +1,55 @@
     9.4 +#ifndef VMATH_RAY_H_
     9.5 +#define VMATH_RAY_H_
     9.6 +
     9.7 +#include "matrix.h"
     9.8 +#include "vector.h"
     9.9 +
    9.10 +typedef struct {
    9.11 +	vec3_t origin, dir;
    9.12 +} ray_t;
    9.13 +
    9.14 +#ifdef __cplusplus
    9.15 +extern "C" {
    9.16 +#endif	/* __cplusplus */
    9.17 +
    9.18 +static inline ray_t ray_cons(vec3_t origin, vec3_t dir);
    9.19 +ray_t ray_transform(ray_t r, mat4_t m);
    9.20 +
    9.21 +#ifdef __cplusplus
    9.22 +}	/* __cplusplus */
    9.23 +
    9.24 +#include <stack>
    9.25 +
    9.26 +class Ray {
    9.27 +private:
    9.28 +	std::stack<scalar_t> ior_stack;
    9.29 +
    9.30 +public:
    9.31 +	/* enviornmental index of refraction, normally 1.0 */
    9.32 +	static scalar_t env_ior;
    9.33 +	
    9.34 +	Vector3 origin, dir;
    9.35 +	scalar_t energy;
    9.36 +	int iter;
    9.37 +	scalar_t ior;
    9.38 +	long time;
    9.39 +
    9.40 +	Ray();
    9.41 +	Ray(const Vector3 &origin, const Vector3 &dir);
    9.42 +
    9.43 +	void transform(const Matrix4x4 &xform);
    9.44 +	Ray transformed(const Matrix4x4 &xform) const;
    9.45 +
    9.46 +	void enter(scalar_t new_ior);
    9.47 +	void leave();
    9.48 +
    9.49 +	scalar_t calc_ior(bool entering, scalar_t mat_ior = 1.0) const;
    9.50 +};
    9.51 +
    9.52 +inline Ray reflect_ray(const Ray &inray, const Vector3 &norm);
    9.53 +inline Ray refract_ray(const Ray &inray, const Vector3 &norm, scalar_t ior, bool entering, scalar_t ray_mag = -1.0);
    9.54 +#endif	/* __cplusplus */
    9.55 +
    9.56 +#include "ray.inl"
    9.57 +
    9.58 +#endif	/* VMATH_RAY_H_ */

    10.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    10.2 +++ b/libs/vmath/ray_c.c	Thu Sep 08 08:30:42 2011 +0300
    10.3 @@ -0,0 +1,18 @@
    10.4 +#include "ray.h"
    10.5 +#include "vector.h"
    10.6 +
    10.7 +ray_t ray_transform(ray_t r, mat4_t xform)
    10.8 +{
    10.9 +	mat4_t upper;
   10.10 +	vec3_t dir;
   10.11 +
   10.12 +	m4_copy(upper, xform);
   10.13 +	upper[0][3] = upper[1][3] = upper[2][3] = upper[3][0] = upper[3][1] = upper[3][2] = 0.0;
   10.14 +	upper[3][3] = 1.0;
   10.15 +
   10.16 +	dir = v3_sub(r.dir, r.origin);
   10.17 +	dir = v3_transform(dir, upper);
   10.18 +	r.origin = v3_transform(r.origin, xform);
   10.19 +	r.dir = v3_add(dir, r.origin);
   10.20 +	return r;
   10.21 +}

    11.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    11.2 +++ b/libs/vmath/sphvec.h	Thu Sep 08 08:30:42 2011 +0300
    11.3 @@ -0,0 +1,18 @@
    11.4 +#ifndef VMATH_SPHVEC_H_
    11.5 +#define VMATH_SPHVEC_H_
    11.6 +
    11.7 +#include "vmath_types.h"
    11.8 +
    11.9 +#ifdef __cplusplus
   11.10 +/* Vector in spherical coordinates */
   11.11 +class SphVector {
   11.12 +public:
   11.13 +	scalar_t theta, phi, r;
   11.14 +
   11.15 +	SphVector(scalar_t theta = 0.0, scalar_t phi = 0.0, scalar_t r = 1.0);
   11.16 +	SphVector(const Vector3 &cvec);
   11.17 +	SphVector &operator =(const Vector3 &cvec);
   11.18 +};
   11.19 +#endif	/* __cplusplus */
   11.20 +
   11.21 +#endif	/* VMATH_SPHVEC_H_ */

    12.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    12.2 +++ b/libs/vmath/vector.h	Thu Sep 08 08:30:42 2011 +0300
    12.3 @@ -0,0 +1,274 @@
    12.4 +#ifndef VMATH_VECTOR_H_
    12.5 +#define VMATH_VECTOR_H_
    12.6 +
    12.7 +#include <stdio.h>
    12.8 +#include "vmath_types.h"
    12.9 +
   12.10 +#ifdef __cplusplus
   12.11 +extern "C" {
   12.12 +#endif	/* __cplusplus */
   12.13 +
   12.14 +/* C 2D vector functions */
   12.15 +static inline vec2_t v2_cons(scalar_t x, scalar_t y);
   12.16 +static inline void v2_print(FILE *fp, vec2_t v);
   12.17 +
   12.18 +static inline vec2_t v2_add(vec2_t v1, vec2_t v2);
   12.19 +static inline vec2_t v2_sub(vec2_t v1, vec2_t v2);
   12.20 +static inline vec2_t v2_scale(vec2_t v, scalar_t s);
   12.21 +static inline scalar_t v2_dot(vec2_t v1, vec2_t v2);
   12.22 +static inline scalar_t v2_length(vec2_t v);
   12.23 +static inline scalar_t v2_length_sq(vec2_t v);
   12.24 +static inline vec2_t v2_normalize(vec2_t v);
   12.25 +
   12.26 +static inline vec2_t v2_lerp(vec2_t v1, vec2_t v2, scalar_t t);
   12.27 +
   12.28 +/* C 3D vector functions */
   12.29 +static inline vec3_t v3_cons(scalar_t x, scalar_t y, scalar_t z);
   12.30 +static inline void v3_print(FILE *fp, vec3_t v);
   12.31 +
   12.32 +static inline vec3_t v3_add(vec3_t v1, vec3_t v2);
   12.33 +static inline vec3_t v3_sub(vec3_t v1, vec3_t v2);
   12.34 +static inline vec3_t v3_neg(vec3_t v);
   12.35 +static inline vec3_t v3_mul(vec3_t v1, vec3_t v2);
   12.36 +static inline vec3_t v3_scale(vec3_t v1, scalar_t s);
   12.37 +static inline scalar_t v3_dot(vec3_t v1, vec3_t v2);
   12.38 +static inline vec3_t v3_cross(vec3_t v1, vec3_t v2);
   12.39 +static inline scalar_t v3_length(vec3_t v);
   12.40 +static inline scalar_t v3_length_sq(vec3_t v);
   12.41 +static inline vec3_t v3_normalize(vec3_t v);
   12.42 +static inline vec3_t v3_transform(vec3_t v, mat4_t m);
   12.43 +
   12.44 +static inline vec3_t v3_rotate(vec3_t v, scalar_t x, scalar_t y, scalar_t z);
   12.45 +static inline vec3_t v3_rotate_axis(vec3_t v, scalar_t angle, scalar_t x, scalar_t y, scalar_t z);
   12.46 +static inline vec3_t v3_rotate_quat(vec3_t v, quat_t q);
   12.47 +
   12.48 +static inline vec3_t v3_reflect(vec3_t v, vec3_t n);
   12.49 +
   12.50 +static inline vec3_t v3_lerp(vec3_t v1, vec3_t v2, scalar_t t);
   12.51 +
   12.52 +/* C 4D vector functions */
   12.53 +static inline vec4_t v4_cons(scalar_t x, scalar_t y, scalar_t z, scalar_t w);
   12.54 +static inline void v4_print(FILE *fp, vec4_t v);
   12.55 +
   12.56 +static inline vec4_t v4_add(vec4_t v1, vec4_t v2);
   12.57 +static inline vec4_t v4_sub(vec4_t v1, vec4_t v2);
   12.58 +static inline vec4_t v4_neg(vec4_t v);
   12.59 +static inline vec4_t v4_mul(vec4_t v1, vec4_t v2);
   12.60 +static inline vec4_t v4_scale(vec4_t v, scalar_t s);
   12.61 +static inline scalar_t v4_dot(vec4_t v1, vec4_t v2);
   12.62 +static inline scalar_t v4_length(vec4_t v);
   12.63 +static inline scalar_t v4_length_sq(vec4_t v);
   12.64 +static inline vec4_t v4_normalize(vec4_t v);
   12.65 +static inline vec4_t v4_transform(vec4_t v, mat4_t m);
   12.66 +
   12.67 +#ifdef __cplusplus
   12.68 +}	/* extern "C" */
   12.69 +
   12.70 +/* when included from C++ source files, also define the vector classes */
   12.71 +#include <iostream>
   12.72 +
   12.73 +/** 2D Vector */
   12.74 +class Vector2 {
   12.75 +public:
   12.76 +	scalar_t x, y;
   12.77 +
   12.78 +	explicit Vector2(scalar_t x = 0.0, scalar_t y = 0.0);
   12.79 +	Vector2(const vec2_t &vec);
   12.80 +	Vector2(const Vector3 &vec);
   12.81 +	Vector2(const Vector4 &vec);
   12.82 +
   12.83 +	inline scalar_t &operator [](int elem);
   12.84 +	inline const scalar_t &operator [](int elem) const;
   12.85 +
   12.86 +	inline scalar_t length() const;
   12.87 +	inline scalar_t length_sq() const;
   12.88 +	void normalize();
   12.89 +	Vector2 normalized() const;
   12.90 +
   12.91 +	void transform(const Matrix3x3 &mat);
   12.92 +	Vector2 transformed(const Matrix3x3 &mat) const;
   12.93 +
   12.94 +	void rotate(scalar_t angle);
   12.95 +	Vector2 rotated(scalar_t angle) const;
   12.96 +
   12.97 +	Vector2 reflection(const Vector2 &normal) const;
   12.98 +	Vector2 refraction(const Vector2 &normal, scalar_t src_ior, scalar_t dst_ior) const;
   12.99 +};
  12.100 +
  12.101 +/* unary operations */
  12.102 +inline Vector2 operator -(const Vector2 &vec);
  12.103 +
  12.104 +/* binary vector (op) vector operations */
  12.105 +inline scalar_t dot_product(const Vector2 &v1, const Vector2 &v2);
  12.106 +
  12.107 +inline Vector2 operator +(const Vector2 &v1, const Vector2 &v2);
  12.108 +inline Vector2 operator -(const Vector2 &v1, const Vector2 &v2);
  12.109 +inline Vector2 operator *(const Vector2 &v1, const Vector2 &v2);
  12.110 +inline Vector2 operator /(const Vector2 &v1, const Vector2 &v2);
  12.111 +inline bool operator ==(const Vector2 &v1, const Vector2 &v2);
  12.112 +
  12.113 +inline void operator +=(Vector2 &v1, const Vector2 &v2);
  12.114 +inline void operator -=(Vector2 &v1, const Vector2 &v2);
  12.115 +inline void operator *=(Vector2 &v1, const Vector2 &v2);
  12.116 +inline void operator /=(Vector2 &v1, const Vector2 &v2);
  12.117 +
  12.118 +/* binary vector (op) scalar and scalar (op) vector operations */
  12.119 +inline Vector2 operator +(const Vector2 &vec, scalar_t scalar);
  12.120 +inline Vector2 operator +(scalar_t scalar, const Vector2 &vec);
  12.121 +inline Vector2 operator -(const Vector2 &vec, scalar_t scalar);
  12.122 +inline Vector2 operator *(const Vector2 &vec, scalar_t scalar);
  12.123 +inline Vector2 operator *(scalar_t scalar, const Vector2 &vec);
  12.124 +inline Vector2 operator /(const Vector2 &vec, scalar_t scalar);
  12.125 +
  12.126 +inline void operator +=(Vector2 &vec, scalar_t scalar);
  12.127 +inline void operator -=(Vector2 &vec, scalar_t scalar);
  12.128 +inline void operator *=(Vector2 &vec, scalar_t scalar);
  12.129 +inline void operator /=(Vector2 &vec, scalar_t scalar);
  12.130 +
  12.131 +std::ostream &operator <<(std::ostream &out, const Vector2 &vec);
  12.132 +
  12.133 +inline Vector2 lerp(const Vector2 &a, const Vector2 &b, scalar_t t);
  12.134 +inline Vector2 catmull_rom_spline(const Vector2 &v0, const Vector2 &v1,
  12.135 +		const Vector2 &v2, const Vector2 &v3, scalar_t t);
  12.136 +
  12.137 +/* 3D Vector */
  12.138 +class Vector3 {
  12.139 +public:
  12.140 +	scalar_t x, y, z;
  12.141 +
  12.142 +	explicit Vector3(scalar_t x = 0.0, scalar_t y = 0.0, scalar_t z = 0.0);
  12.143 +	Vector3(const vec3_t &vec);
  12.144 +	Vector3(const Vector2 &vec);
  12.145 +	Vector3(const Vector4 &vec);
  12.146 +	Vector3(const SphVector &sph);
  12.147 +
  12.148 +	Vector3 &operator =(const SphVector &sph);
  12.149 +
  12.150 +	inline scalar_t &operator [](int elem);
  12.151 +	inline const scalar_t &operator [](int elem) const;
  12.152 +
  12.153 +	inline scalar_t length() const;
  12.154 +	inline scalar_t length_sq() const;
  12.155 +	void normalize();
  12.156 +	Vector3 normalized() const;
  12.157 +
  12.158 +	void transform(const Matrix3x3 &mat);
  12.159 +	Vector3 transformed(const Matrix3x3 &mat) const;
  12.160 +	void transform(const Matrix4x4 &mat);
  12.161 +	Vector3 transformed(const Matrix4x4 &mat) const;
  12.162 +	void transform(const Quaternion &quat);
  12.163 +	Vector3 transformed(const Quaternion &quat) const;
  12.164 +
  12.165 +	void rotate(const Vector3 &euler);
  12.166 +	Vector3 rotated(const Vector3 &euler) const;
  12.167 +
  12.168 +	Vector3 reflection(const Vector3 &normal) const;
  12.169 +	Vector3 refraction(const Vector3 &normal, scalar_t src_ior, scalar_t dst_ior) const;
  12.170 +	Vector3 refraction(const Vector3 &normal, scalar_t ior) const;
  12.171 +};
  12.172 +
  12.173 +/* unary operations */
  12.174 +inline Vector3 operator -(const Vector3 &vec);
  12.175 +
  12.176 +/* binary vector (op) vector operations */
  12.177 +inline scalar_t dot_product(const Vector3 &v1, const Vector3 &v2);
  12.178 +inline Vector3 cross_product(const Vector3 &v1, const Vector3 &v2);
  12.179 +
  12.180 +inline Vector3 operator +(const Vector3 &v1, const Vector3 &v2);
  12.181 +inline Vector3 operator -(const Vector3 &v1, const Vector3 &v2);
  12.182 +inline Vector3 operator *(const Vector3 &v1, const Vector3 &v2);
  12.183 +inline Vector3 operator /(const Vector3 &v1, const Vector3 &v2);
  12.184 +inline bool operator ==(const Vector3 &v1, const Vector3 &v2);
  12.185 +
  12.186 +inline void operator +=(Vector3 &v1, const Vector3 &v2);
  12.187 +inline void operator -=(Vector3 &v1, const Vector3 &v2);
  12.188 +inline void operator *=(Vector3 &v1, const Vector3 &v2);
  12.189 +inline void operator /=(Vector3 &v1, const Vector3 &v2);
  12.190 +
  12.191 +/* binary vector (op) scalar and scalar (op) vector operations */
  12.192 +inline Vector3 operator +(const Vector3 &vec, scalar_t scalar);
  12.193 +inline Vector3 operator +(scalar_t scalar, const Vector3 &vec);
  12.194 +inline Vector3 operator -(const Vector3 &vec, scalar_t scalar);
  12.195 +inline Vector3 operator *(const Vector3 &vec, scalar_t scalar);
  12.196 +inline Vector3 operator *(scalar_t scalar, const Vector3 &vec);
  12.197 +inline Vector3 operator /(const Vector3 &vec, scalar_t scalar);
  12.198 +
  12.199 +inline void operator +=(Vector3 &vec, scalar_t scalar);
  12.200 +inline void operator -=(Vector3 &vec, scalar_t scalar);
  12.201 +inline void operator *=(Vector3 &vec, scalar_t scalar);
  12.202 +inline void operator /=(Vector3 &vec, scalar_t scalar);
  12.203 +
  12.204 +std::ostream &operator <<(std::ostream &out, const Vector3 &vec);
  12.205 +
  12.206 +inline Vector3 lerp(const Vector3 &a, const Vector3 &b, scalar_t t);
  12.207 +inline Vector3 catmull_rom_spline(const Vector3 &v0, const Vector3 &v1,
  12.208 +		const Vector3 &v2, const Vector3 &v3, scalar_t t);
  12.209 +
  12.210 +/* 4D Vector */
  12.211 +class Vector4 {
  12.212 +public:
  12.213 +	scalar_t x, y, z, w;
  12.214 +
  12.215 +	explicit Vector4(scalar_t x = 0.0, scalar_t y = 0.0, scalar_t z = 0.0, scalar_t w = 0.0);
  12.216 +	Vector4(const vec4_t &vec);
  12.217 +	Vector4(const Vector2 &vec);
  12.218 +	Vector4(const Vector3 &vec);
  12.219 +
  12.220 +	inline scalar_t &operator [](int elem);
  12.221 +	inline const scalar_t &operator [](int elem) const;
  12.222 +
  12.223 +	inline scalar_t length() const;
  12.224 +	inline scalar_t length_sq() const;
  12.225 +	void normalize();
  12.226 +	Vector4 normalized() const;
  12.227 +
  12.228 +	void transform(const Matrix4x4 &mat);
  12.229 +	Vector4 transformed(const Matrix4x4 &mat) const;
  12.230 +
  12.231 +	Vector4 reflection(const Vector4 &normal) const;
  12.232 +	Vector4 refraction(const Vector4 &normal, scalar_t src_ior, scalar_t dst_ior) const;
  12.233 +};
  12.234 +
  12.235 +
  12.236 +/* unary operations */
  12.237 +inline Vector4 operator -(const Vector4 &vec);
  12.238 +
  12.239 +/* binary vector (op) vector operations */
  12.240 +inline scalar_t dot_product(const Vector4 &v1, const Vector4 &v2);
  12.241 +inline Vector4 cross_product(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3);
  12.242 +
  12.243 +inline Vector4 operator +(const Vector4 &v1, const Vector4 &v2);
  12.244 +inline Vector4 operator -(const Vector4 &v1, const Vector4 &v2);
  12.245 +inline Vector4 operator *(const Vector4 &v1, const Vector4 &v2);
  12.246 +inline Vector4 operator /(const Vector4 &v1, const Vector4 &v2);
  12.247 +inline bool operator ==(const Vector4 &v1, const Vector4 &v2);
  12.248 +
  12.249 +inline void operator +=(Vector4 &v1, const Vector4 &v2);
  12.250 +inline void operator -=(Vector4 &v1, const Vector4 &v2);
  12.251 +inline void operator *=(Vector4 &v1, const Vector4 &v2);
  12.252 +inline void operator /=(Vector4 &v1, const Vector4 &v2);
  12.253 +
  12.254 +/* binary vector (op) scalar and scalar (op) vector operations */
  12.255 +inline Vector4 operator +(const Vector4 &vec, scalar_t scalar);
  12.256 +inline Vector4 operator +(scalar_t scalar, const Vector4 &vec);
  12.257 +inline Vector4 operator -(const Vector4 &vec, scalar_t scalar);
  12.258 +inline Vector4 operator *(const Vector4 &vec, scalar_t scalar);
  12.259 +inline Vector4 operator *(scalar_t scalar, const Vector4 &vec);
  12.260 +inline Vector4 operator /(const Vector4 &vec, scalar_t scalar);
  12.261 +
  12.262 +inline void operator +=(Vector4 &vec, scalar_t scalar);
  12.263 +inline void operator -=(Vector4 &vec, scalar_t scalar);
  12.264 +inline void operator *=(Vector4 &vec, scalar_t scalar);
  12.265 +inline void operator /=(Vector4 &vec, scalar_t scalar);
  12.266 +
  12.267 +std::ostream &operator <<(std::ostream &out, const Vector4 &vec);
  12.268 +
  12.269 +inline Vector4 lerp(const Vector4 &v0, const Vector4 &v1, scalar_t t);
  12.270 +inline Vector4 catmull_rom_spline(const Vector4 &v0, const Vector4 &v1,
  12.271 +		const Vector4 &v2, const Vector4 &v3, scalar_t t);
  12.272 +
  12.273 +#endif	/* __cplusplus */
  12.274 +
  12.275 +#include "vector.inl"
  12.276 +
  12.277 +#endif	/* VMATH_VECTOR_H_ */

    13.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    13.2 +++ b/libs/vmath/vmath.c	Thu Sep 08 08:30:42 2011 +0300
    13.3 @@ -0,0 +1,317 @@
    13.4 +#include <stdlib.h>
    13.5 +#include <math.h>
    13.6 +#include "vmath.h"
    13.7 +
    13.8 +/** Numerical calculation of integrals using simpson's rule */
    13.9 +scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples)
   13.10 +{
   13.11 +	int i;
   13.12 +	scalar_t h = (high - low) / (scalar_t)samples;
   13.13 +	scalar_t sum = 0.0;
   13.14 +
   13.15 +	for(i=0; i<samples+1; i++) {
   13.16 +		scalar_t y = f((scalar_t)i * h + low);
   13.17 +		sum += ((!i || i == samples) ? y : ((i % 2) ? 4.0 * y : 2.0 * y)) * (h / 3.0);
   13.18 +	}
   13.19 +	return sum;
   13.20 +}
   13.21 +
   13.22 +/** Gaussuan function */
   13.23 +scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev)
   13.24 +{
   13.25 +	scalar_t exponent = -SQ(x - mean) / (2.0 * SQ(sdev));
   13.26 +	return 1.0 - -pow(M_E, exponent) / (sdev * sqrt(TWO_PI));
   13.27 +}
   13.28 +
   13.29 +
   13.30 +/** b-spline approximation */
   13.31 +scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
   13.32 +{
   13.33 +	vec4_t tmp;
   13.34 +	scalar_t tsq = t * t;
   13.35 +
   13.36 +	static mat4_t bspline_mat = {
   13.37 +		{-1,  3, -3,  1},
   13.38 +		{3, -6,  3,  0},
   13.39 +		{-3,  0,  3,  0},
   13.40 +		{1,  4,  1,  0}
   13.41 +	};
   13.42 +
   13.43 +	tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), bspline_mat), 1.0 / 6.0);
   13.44 +	return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);
   13.45 +}
   13.46 +
   13.47 +/** Catmull-rom spline interpolation */
   13.48 +scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t) {
   13.49 +	vec4_t tmp;
   13.50 +	scalar_t tsq = t * t;
   13.51 +
   13.52 +	static mat4_t crspline_mat = {
   13.53 +		{-1,  3, -3,  1},
   13.54 +		{2, -5,  4, -1},
   13.55 +		{-1,  0,  1,  0},
   13.56 +		{0,  2,  0,  0}
   13.57 +	};
   13.58 +
   13.59 +	tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), crspline_mat), 0.5);
   13.60 +	return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);
   13.61 +}
   13.62 +
   13.63 +/** Bezier interpolation */
   13.64 +scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)
   13.65 +{
   13.66 +	scalar_t omt, omt3, t3, f;
   13.67 +	t3 = t * t * t;
   13.68 +	omt = 1.0f - t;
   13.69 +	omt3 = omt * omt * omt;
   13.70 +	f = 3 * t * omt;
   13.71 +
   13.72 +	return (a * omt3) + (b * f * omt) + (c * f * t) + (d * t3);
   13.73 +}
   13.74 +
   13.75 +/* ---- Ken Perlin's implementation of noise ---- */
   13.76 +
   13.77 +#define B	0x100
   13.78 +#define BM	0xff
   13.79 +#define N	0x1000
   13.80 +#define NP	12   /* 2^N */
   13.81 +#define NM	0xfff
   13.82 +
   13.83 +#define s_curve(t) (t * t * (3.0f - 2.0f * t))
   13.84 +
   13.85 +#define setup(elem, b0, b1, r0, r1) \
   13.86 +	do {							\
   13.87 +		scalar_t t = elem + N;		\
   13.88 +		b0 = ((int)t) & BM;			\
   13.89 +		b1 = (b0 + 1) & BM;			\
   13.90 +		r0 = t - (int)t;			\
   13.91 +		r1 = r0 - 1.0f;				\
   13.92 +	} while(0)
   13.93 +
   13.94 +
   13.95 +static int perm[B + B + 2];			/* permuted index from g_n onto themselves */
   13.96 +static vec3_t grad3[B + B + 2];		/* 3D random gradients */
   13.97 +static vec2_t grad2[B + B + 2];		/* 2D random gradients */
   13.98 +static scalar_t grad1[B + B + 2];	/* 1D random ... slopes */
   13.99 +static int tables_valid;
  13.100 +
  13.101 +static void init_noise()
  13.102 +{
  13.103 +	int i;
  13.104 +
  13.105 +	/* calculate random gradients */
  13.106 +	for(i=0; i<B; i++) {
  13.107 +		perm[i] = i;	/* .. and initialize permutation mapping to identity */
  13.108 +
  13.109 +		grad1[i] = (scalar_t)((rand() % (B + B)) - B) / B;
  13.110 +
  13.111 +		grad2[i].x = (scalar_t)((rand() % (B + B)) - B) / B;
  13.112 +		grad2[i].y = (scalar_t)((rand() % (B + B)) - B) / B;
  13.113 +		grad2[i] = v2_normalize(grad2[i]);
  13.114 +
  13.115 +		grad3[i].x = (scalar_t)((rand() % (B + B)) - B) / B;
  13.116 +		grad3[i].y = (scalar_t)((rand() % (B + B)) - B) / B;
  13.117 +		grad3[i].z = (scalar_t)((rand() % (B + B)) - B) / B;
  13.118 +		grad3[i] = v3_normalize(grad3[i]);
  13.119 +	}
  13.120 +
  13.121 +	/* permute indices by swapping them randomly */
  13.122 +	for(i=0; i<B; i++) {
  13.123 +		int rand_idx = rand() % B;
  13.124 +
  13.125 +		int tmp = perm[i];
  13.126 +		perm[i] = perm[rand_idx];
  13.127 +		perm[rand_idx] = tmp;
  13.128 +	}
  13.129 +
  13.130 +	/* fill up the rest of the arrays by duplicating the existing gradients */
  13.131 +	/* and permutations */
  13.132 +	for(i=0; i<B+2; i++) {
  13.133 +		perm[B + i] = perm[i];
  13.134 +		grad1[B + i] = grad1[i];
  13.135 +		grad2[B + i] = grad2[i];
  13.136 +		grad3[B + i] = grad3[i];
  13.137 +	}
  13.138 +}
  13.139 +
  13.140 +scalar_t noise1(scalar_t x)
  13.141 +{
  13.142 +	int bx0, bx1;
  13.143 +	scalar_t rx0, rx1, sx, u, v;
  13.144 +
  13.145 +	if(!tables_valid) {
  13.146 +		init_noise();
  13.147 +		tables_valid = 1;
  13.148 +	}
  13.149 +
  13.150 +	setup(x, bx0, bx1, rx0, rx1);
  13.151 +	sx = s_curve(rx0);
  13.152 +	u = rx0 * grad1[perm[bx0]];
  13.153 +	v = rx1 * grad1[perm[bx1]];
  13.154 +
  13.155 +	return lerp(u, v, sx);
  13.156 +}
  13.157 +
  13.158 +scalar_t noise2(scalar_t x, scalar_t y)
  13.159 +{
  13.160 +	int i, j, b00, b10, b01, b11;
  13.161 +	int bx0, bx1, by0, by1;
  13.162 +	scalar_t rx0, rx1, ry0, ry1;
  13.163 +	scalar_t sx, sy, u, v, a, b;
  13.164 +
  13.165 +	if(!tables_valid) {
  13.166 +		init_noise();
  13.167 +		tables_valid = 1;
  13.168 +	}
  13.169 +
  13.170 +	setup(x, bx0, bx1, rx0, rx1);
  13.171 +	setup(y, by0, by1, ry0, ry1);
  13.172 +
  13.173 +	i = perm[bx0];
  13.174 +	j = perm[bx1];
  13.175 +
  13.176 +	b00 = perm[i + by0];
  13.177 +	b10 = perm[j + by0];
  13.178 +	b01 = perm[i + by1];
  13.179 +	b11 = perm[j + by1];
  13.180 +
  13.181 +	/* calculate hermite inteprolating factors */
  13.182 +	sx = s_curve(rx0);
  13.183 +	sy = s_curve(ry0);
  13.184 +
  13.185 +	/* interpolate along the left edge */
  13.186 +	u = v2_dot(grad2[b00], v2_cons(rx0, ry0));
  13.187 +	v = v2_dot(grad2[b10], v2_cons(rx1, ry0));
  13.188 +	a = lerp(u, v, sx);
  13.189 +
  13.190 +	/* interpolate along the right edge */
  13.191 +	u = v2_dot(grad2[b01], v2_cons(rx0, ry1));
  13.192 +	v = v2_dot(grad2[b11], v2_cons(rx1, ry1));
  13.193 +	b = lerp(u, v, sx);
  13.194 +
  13.195 +	/* interpolate between them */
  13.196 +	return lerp(a, b, sy);
  13.197 +}
  13.198 +
  13.199 +scalar_t noise3(scalar_t x, scalar_t y, scalar_t z)
  13.200 +{
  13.201 +	int i, j;
  13.202 +	int bx0, bx1, by0, by1, bz0, bz1;
  13.203 +	int b00, b10, b01, b11;
  13.204 +	scalar_t rx0, rx1, ry0, ry1, rz0, rz1;
  13.205 +	scalar_t sx, sy, sz;
  13.206 +	scalar_t u, v, a, b, c, d;
  13.207 +
  13.208 +	if(!tables_valid) {
  13.209 +		init_noise();
  13.210 +		tables_valid = 1;
  13.211 +	}
  13.212 +
  13.213 +	setup(x, bx0, bx1, rx0, rx1);
  13.214 +	setup(y, by0, by1, ry0, ry1);
  13.215 +	setup(z, bz0, bz1, rz0, rz1);
  13.216 +
  13.217 +	i = perm[bx0];
  13.218 +	j = perm[bx1];
  13.219 +
  13.220 +	b00 = perm[i + by0];
  13.221 +	b10 = perm[j + by0];
  13.222 +	b01 = perm[i + by1];
  13.223 +	b11 = perm[j + by1];
  13.224 +
  13.225 +	/* calculate hermite interpolating factors */
  13.226 +	sx = s_curve(rx0);
  13.227 +	sy = s_curve(ry0);
  13.228 +	sz = s_curve(rz0);
  13.229 +
  13.230 +	/* interpolate along the top slice of the cell */
  13.231 +	u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz0));
  13.232 +	v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz0));
  13.233 +	a = lerp(u, v, sx);
  13.234 +
  13.235 +	u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz0));
  13.236 +	v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz0));
  13.237 +	b = lerp(u, v, sx);
  13.238 +
  13.239 +	c = lerp(a, b, sy);
  13.240 +
  13.241 +	/* interpolate along the bottom slice of the cell */
  13.242 +	u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz1));
  13.243 +	v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz1));
  13.244 +	a = lerp(u, v, sx);
  13.245 +
  13.246 +	u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz1));
  13.247 +	v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz1));
  13.248 +	b = lerp(u, v, sx);
  13.249 +
  13.250 +	d = lerp(a, b, sy);
  13.251 +
  13.252 +	/* interpolate between slices */
  13.253 +	return lerp(c, d, sz);
  13.254 +}
  13.255 +
  13.256 +scalar_t fbm1(scalar_t x, int octaves)
  13.257 +{
  13.258 +	int i;
  13.259 +	scalar_t res = 0.0f, freq = 1.0f;
  13.260 +	for(i=0; i<octaves; i++) {
  13.261 +		res += noise1(x * freq) / freq;
  13.262 +		freq *= 2.0f;
  13.263 +	}
  13.264 +	return res;
  13.265 +}
  13.266 +
  13.267 +scalar_t fbm2(scalar_t x, scalar_t y, int octaves)
  13.268 +{
  13.269 +	int i;
  13.270 +	scalar_t res = 0.0f, freq = 1.0f;
  13.271 +	for(i=0; i<octaves; i++) {
  13.272 +		res += noise2(x * freq, y * freq) / freq;
  13.273 +		freq *= 2.0f;
  13.274 +	}
  13.275 +	return res;
  13.276 +}
  13.277 +
  13.278 +scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves)
  13.279 +{
  13.280 +	int i;
  13.281 +	scalar_t res = 0.0f, freq = 1.0f;
  13.282 +	for(i=0; i<octaves; i++) {
  13.283 +		res += noise3(x * freq, y * freq, z * freq) / freq;
  13.284 +		freq *= 2.0f;
  13.285 +	}
  13.286 +	return res;
  13.287 +}
  13.288 +
  13.289 +scalar_t turbulence1(scalar_t x, int octaves)
  13.290 +{
  13.291 +	int i;
  13.292 +	scalar_t res = 0.0f, freq = 1.0f;
  13.293 +	for(i=0; i<octaves; i++) {
  13.294 +		res += fabs(noise1(x * freq) / freq);
  13.295 +		freq *= 2.0f;
  13.296 +	}
  13.297 +	return res;
  13.298 +}
  13.299 +
  13.300 +scalar_t turbulence2(scalar_t x, scalar_t y, int octaves)
  13.301 +{
  13.302 +	int i;
  13.303 +	scalar_t res = 0.0f, freq = 1.0f;
  13.304 +	for(i=0; i<octaves; i++) {
  13.305 +		res += fabs(noise2(x * freq, y * freq) / freq);
  13.306 +		freq *= 2.0f;
  13.307 +	}
  13.308 +	return res;
  13.309 +}
  13.310 +
  13.311 +scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves)
  13.312 +{
  13.313 +	int i;
  13.314 +	scalar_t res = 0.0f, freq = 1.0f;
  13.315 +	for(i=0; i<octaves; i++) {
  13.316 +		res += fabs(noise3(x * freq, y * freq, z * freq) / freq);
  13.317 +		freq *= 2.0f;
  13.318 +	}
  13.319 +	return res;
  13.320 +}

    14.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    14.2 +++ b/libs/vmath/vmath.h	Thu Sep 08 08:30:42 2011 +0300
    14.3 @@ -0,0 +1,74 @@
    14.4 +#ifndef VMATH_H_
    14.5 +#define VMATH_H_
    14.6 +
    14.7 +#include <math.h>
    14.8 +#include "vmath_types.h"
    14.9 +
   14.10 +#ifndef M_PI
   14.11 +#define M_PI	PI
   14.12 +#endif
   14.13 +
   14.14 +#ifndef M_E
   14.15 +#define M_E				2.718281828459045
   14.16 +#endif
   14.17 +
   14.18 +#define PI				3.141592653589793
   14.19 +#define HALF_PI			1.570796326794897
   14.20 +#define QUARTER_PI		0.785398163397448
   14.21 +#define TWO_PI			6.283185307179586
   14.22 +
   14.23 +
   14.24 +#define RAD_TO_DEG(a) ((((scalar_t)a) * 360.0) / TWO_PI)
   14.25 +#define DEG_TO_RAD(a) (((scalar_t)a) * (PI / 180.0))
   14.26 +
   14.27 +#define SQ(x) ((x) * (x))
   14.28 +
   14.29 +#define MIN(a, b)	((a) < (b) ? (a) : (b))
   14.30 +#define MAX(a, b)	((a) > (b) ? (a) : (b))
   14.31 +
   14.32 +#ifndef __GNUC__
   14.33 +#define round(x)	((x) >= 0 ? (x) + 0.5 : (x) - 0.5)
   14.34 +#endif
   14.35 +
   14.36 +#ifdef __cplusplus
   14.37 +extern "C" {
   14.38 +#endif	/* __cplusplus */
   14.39 +
   14.40 +static inline scalar_t frand(scalar_t range);
   14.41 +static inline vec3_t sphrand(scalar_t rad);
   14.42 +
   14.43 +scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples);
   14.44 +scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev);
   14.45 +
   14.46 +static inline scalar_t lerp(scalar_t a, scalar_t b, scalar_t t);
   14.47 +
   14.48 +scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
   14.49 +scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
   14.50 +scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t);
   14.51 +
   14.52 +scalar_t noise1(scalar_t x);
   14.53 +scalar_t noise2(scalar_t x, scalar_t y);
   14.54 +scalar_t noise3(scalar_t x, scalar_t y, scalar_t z);
   14.55 +
   14.56 +scalar_t fbm1(scalar_t x, int octaves);
   14.57 +scalar_t fbm2(scalar_t x, scalar_t y, int octaves);
   14.58 +scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves);
   14.59 +
   14.60 +scalar_t turbulence1(scalar_t x, int octaves);
   14.61 +scalar_t turbulence2(scalar_t x, scalar_t y, int octaves);
   14.62 +scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves);
   14.63 +
   14.64 +#ifdef __cplusplus
   14.65 +}
   14.66 +#endif	/* __cplusplus */
   14.67 +
   14.68 +#include "vmath.inl"
   14.69 +
   14.70 +#include "vector.h"
   14.71 +#include "matrix.h"
   14.72 +#include "quat.h"
   14.73 +#include "sphvec.h"
   14.74 +#include "ray.h"
   14.75 +#include "geom.h"
   14.76 +
   14.77 +#endif	/* VMATH_H_ */

    15.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    15.2 +++ b/libs/vmath/vmath_config.h	Thu Sep 08 08:30:42 2011 +0300
    15.3 @@ -0,0 +1,19 @@
    15.4 +#ifndef VMATH_CONFIG_H_
    15.5 +#define VMATH_CONFIG_H_
    15.6 +
    15.7 +#if (__STDC_VERSION__ < 199999)
    15.8 +#if defined(__GNUC__) || defined(_MSC_VER)
    15.9 +#define inline __inline
   15.10 +#else
   15.11 +#define inline
   15.12 +
   15.13 +#ifdef VECTOR_H_
   15.14 +#warning "compiling vector operations without inline, performance might suffer"
   15.15 +#endif	/* VECTOR_H_ */
   15.16 +
   15.17 +#endif	/* gcc/msvc */
   15.18 +#endif	/* not C99 */
   15.19 +
   15.20 +#define SINGLE_PRECISION_MATH
   15.21 +
   15.22 +#endif	/* VMATH_CONFIG_H_ */

    16.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
    16.2 +++ b/libs/vmath/vmath_types.h	Thu Sep 08 08:30:42 2011 +0300
    16.3 @@ -0,0 +1,40 @@
    16.4 +#ifndef VMATH_TYPES_H_
    16.5 +#define VMATH_TYPES_H_
    16.6 +
    16.7 +#include "vmath_config.h"
    16.8 +
    16.9 +#define SMALL_NUMBER	1.e-4
   16.10 +#define XSMALL_NUMBER	1.e-8
   16.11 +#define ERROR_MARGIN	1.e-6
   16.12 +
   16.13 +
   16.14 +#ifdef SINGLE_PRECISION_MATH
   16.15 +typedef float scalar_t;
   16.16 +#else
   16.17 +typedef double scalar_t;
   16.18 +#endif	/* floating point precision */
   16.19 +
   16.20 +/* vectors */
   16.21 +typedef struct { scalar_t x, y; } vec2_t;
   16.22 +typedef struct { scalar_t x, y, z; } vec3_t;
   16.23 +typedef struct { scalar_t x, y, z, w; } vec4_t;
   16.24 +
   16.25 +/* quaternions */
   16.26 +typedef vec4_t quat_t;
   16.27 +
   16.28 +/* matrices */
   16.29 +typedef scalar_t mat3_t[3][3];
   16.30 +typedef scalar_t mat4_t[4][4];
   16.31 +
   16.32 +
   16.33 +#ifdef __cplusplus
   16.34 +class Vector2;
   16.35 +class Vector3;
   16.36 +class Vector4;
   16.37 +class Quaternion;
   16.38 +class Matrix3x3;
   16.39 +class Matrix4x4;
   16.40 +class SphVector;
   16.41 +#endif	/* __cplusplus */
   16.42 +
   16.43 +#endif	/* VMATH_TYPES_H_ */