1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/bezmath.c	Sun Oct 27 06:31:18 2013 +0200
     1.3 @@ -0,0 +1,93 @@
     1.4 +#include <math.h>
     1.5 +#include "bezmath.h"
     1.6 +
     1.7 +static float bernstein(int i, float x);
     1.8 +
     1.9 +
    1.10 +struct vec3 v3_add(struct vec3 a, struct vec3 b)
    1.11 +{
    1.12 +	a.x += b.x;
    1.13 +	a.y += b.y;
    1.14 +	a.z += b.z;
    1.15 +	return a;
    1.16 +}
    1.17 +
    1.18 +struct vec3 v3_sub(struct vec3 a, struct vec3 b)
    1.19 +{
    1.20 +	a.x -= b.x;
    1.21 +	a.y -= b.y;
    1.22 +	a.z -= b.z;
    1.23 +	return a;
    1.24 +}
    1.25 +
    1.26 +struct vec3 v3_cross(struct vec3 a, struct vec3 b)
    1.27 +{
    1.28 +	struct vec3 res;
    1.29 +	res.x = a.y * b.z - a.z * b.y;
    1.30 +	res.y = a.z * b.x - a.x * b.z;
    1.31 +	res.z = a.x * b.y - a.y * b.x;
    1.32 +	return res;
    1.33 +}
    1.34 +
    1.35 +struct vec3 v3_normalize(struct vec3 v)
    1.36 +{
    1.37 +	float len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
    1.38 +	v.x /= len;
    1.39 +	v.y /= len;
    1.40 +	v.z /= len;
    1.41 +	return v;
    1.42 +}
    1.43 +
    1.44 +struct vec3 bezier_patch(struct vec3 *cp, float u, float v)
    1.45 +{
    1.46 +	int i, j;
    1.47 +	struct vec3 res = {0, 0, 0};
    1.48 +
    1.49 +	for(j=0; j<4; j++) {
    1.50 +		for(i=0; i<4; i++) {
    1.51 +			float bu = bernstein(i, u);
    1.52 +			float bv = bernstein(j, v);
    1.53 +
    1.54 +			res.x += cp->x * bu * bv;
    1.55 +			res.y += cp->y * bu * bv;
    1.56 +			res.z += cp->z * bu * bv;
    1.57 +
    1.58 +			cp++;
    1.59 +		}
    1.60 +	}
    1.61 +	return res;
    1.62 +}
    1.63 +
    1.64 +#define DT	0.001
    1.65 +
    1.66 +struct vec3 bezier_patch_norm(struct vec3 *cp, float u, float v)
    1.67 +{
    1.68 +	struct vec3 tang, bitan, norm;
    1.69 +
    1.70 +	tang = v3_sub(bezier_patch(cp, u + DT, v), bezier_patch(cp, u - DT, v));
    1.71 +	bitan = v3_sub(bezier_patch(cp, u, v + DT), bezier_patch(cp, u, v - DT));
    1.72 +	norm = v3_cross(tang, bitan);
    1.73 +
    1.74 +	return v3_normalize(norm);
    1.75 +}
    1.76 +
    1.77 +
    1.78 +
    1.79 +static float bernstein(int i, float x)
    1.80 +{
    1.81 +	float invx = 1.0 - x;
    1.82 +
    1.83 +	switch(i) {
    1.84 +	case 0:
    1.85 +		return invx * invx * invx;
    1.86 +	case 1:
    1.87 +		return 3 * x * invx * invx;
    1.88 +	case 2:
    1.89 +		return 3 * x * x * invx;
    1.90 +	case 3:
    1.91 +		return x * x * x;
    1.92 +	default:
    1.93 +		break;
    1.94 +	}
    1.95 +	return 0;
    1.96 +}