nuclear@1: #include <math.h>
nuclear@1: #include "vmath.h"
nuclear@1: 
nuclear@1: static float bernstein(int i, float x);
nuclear@1: 
nuclear@1: 
nuclear@1: struct vec3 v3_add(struct vec3 a, struct vec3 b)
nuclear@1: {
nuclear@1: 	a.x += b.x;
nuclear@1: 	a.y += b.y;
nuclear@1: 	a.z += b.z;
nuclear@1: 	return a;
nuclear@1: }
nuclear@1: 
nuclear@1: struct vec3 v3_sub(struct vec3 a, struct vec3 b)
nuclear@1: {
nuclear@1: 	a.x -= b.x;
nuclear@1: 	a.y -= b.y;
nuclear@1: 	a.z -= b.z;
nuclear@1: 	return a;
nuclear@1: }
nuclear@1: 
nuclear@1: struct vec3 v3_cross(struct vec3 a, struct vec3 b)
nuclear@1: {
nuclear@1: 	struct vec3 res;
nuclear@1: 	res.x = a.y * b.z - a.z * b.y;
nuclear@1: 	res.y = a.z * b.x - a.x * b.z;
nuclear@1: 	res.z = a.x * b.y - a.y * b.x;
nuclear@1: 	return res;
nuclear@1: }
nuclear@1: 
nuclear@1: struct vec3 v3_normalize(struct vec3 v)
nuclear@1: {
nuclear@1: 	float len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
nuclear@1: 	v.x /= len;
nuclear@1: 	v.y /= len;
nuclear@1: 	v.z /= len;
nuclear@1: 	return v;
nuclear@1: }
nuclear@1: 
nuclear@1: struct vec3 bezier_patch(struct vec3 *cp, float u, float v)
nuclear@1: {
nuclear@1: 	int i, j;
nuclear@1: 	struct vec3 res = {0, 0, 0};
nuclear@1: 
nuclear@1: 	for(j=0; j<4; j++) {
nuclear@1: 		for(i=0; i<4; i++) {
nuclear@1: 			float bu = bernstein(i, u);
nuclear@1: 			float bv = bernstein(j, v);
nuclear@1: 
nuclear@1: 			res.x += cp->x * bu * bv;
nuclear@1: 			res.y += cp->y * bu * bv;
nuclear@1: 			res.z += cp->z * bu * bv;
nuclear@1: 
nuclear@1: 			cp++;
nuclear@1: 		}
nuclear@1: 	}
nuclear@1: 	return res;
nuclear@1: }
nuclear@1: 
nuclear@1: #define DT	0.001
nuclear@1: 
nuclear@1: struct vec3 bezier_patch_norm(struct vec3 *cp, float u, float v)
nuclear@1: {
nuclear@1: 	struct vec3 tang, bitan, norm;
nuclear@1: 
nuclear@1: 	tang = v3_sub(bezier_patch(cp, u + DT, v), bezier_patch(cp, u - DT, v));
nuclear@1: 	bitan = v3_sub(bezier_patch(cp, u, v + DT), bezier_patch(cp, u, v - DT));
nuclear@1: 	norm = v3_cross(tang, bitan);
nuclear@1: 
nuclear@1: 	return v3_normalize(norm);
nuclear@1: }
nuclear@1: 
nuclear@1: 
nuclear@1: 
nuclear@1: static float bernstein(int i, float x)
nuclear@1: {
nuclear@1: 	float invx = 1.0 - x;
nuclear@1: 
nuclear@1: 	switch(i) {
nuclear@1: 	case 0:
nuclear@1: 		return invx * invx * invx;
nuclear@1: 	case 1:
nuclear@1: 		return 3 * x * invx * invx;
nuclear@1: 	case 2:
nuclear@1: 		return 3 * x * x * invx;
nuclear@1: 	case 3:
nuclear@1: 		return x * x * x;
nuclear@1: 	default:
nuclear@1: 		break;
nuclear@1: 	}
nuclear@1: 	return 0;
nuclear@1: }