nuclear@0: #include <stdlib.h>
nuclear@0: #include <math.h>
nuclear@0: #include "noise.h"
nuclear@0: 
nuclear@0: typedef struct vec2 {
nuclear@0: 	float x, y;
nuclear@0: } vec2_t;
nuclear@0: 
nuclear@0: typedef struct vec3 {
nuclear@0: 	float x, y, z;
nuclear@0: } vec3_t;
nuclear@0: 
nuclear@0: /* ---- Ken Perlin's implementation of noise ---- */
nuclear@0: 
nuclear@0: #define B	0x100
nuclear@0: #define BM	0xff
nuclear@0: #define N	0x1000
nuclear@0: #define NP	12   /* 2^N */
nuclear@0: #define NM	0xfff
nuclear@0: 
nuclear@0: #define lerp(a, b, t)	((a) + ((b) - (a)) * t)
nuclear@0: #define s_curve(t) (t * t * (3.0f - 2.0f * t))
nuclear@0: 
nuclear@0: #define setup(elem, b0, b1, r0, r1) \
nuclear@0: 	do {							\
nuclear@0: 		float t = elem + N;		\
nuclear@0: 		b0 = ((int)t) & BM;			\
nuclear@0: 		b1 = (b0 + 1) & BM;			\
nuclear@0: 		r0 = t - (int)t;			\
nuclear@0: 		r1 = r0 - 1.0f;				\
nuclear@0: 	} while(0)
nuclear@0: 
nuclear@0: 
nuclear@0: static int perm[B + B + 2];			/* permuted index from g_n onto themselves */
nuclear@0: static vec3_t grad3[B + B + 2];		/* 3D random gradients */
nuclear@0: static vec2_t grad2[B + B + 2];		/* 2D random gradients */
nuclear@0: static float grad1[B + B + 2];	/* 1D random ... slopes */
nuclear@0: static int tables_valid;
nuclear@0: 
nuclear@0: static vec2_t v2_normalize(vec2_t v)
nuclear@0: {
nuclear@0: 	vec2_t res;
nuclear@0: 	float len = sqrt(v.x * v.x + v.y * v.y);
nuclear@0: 	if(len != 0.0f) {
nuclear@0: 		res.x = v.x / len;
nuclear@0: 		res.y = v.y / len;
nuclear@0: 	} else {
nuclear@0: 		res = v;
nuclear@0: 	}
nuclear@0: 	return res;
nuclear@0: }
nuclear@0: 
nuclear@0: static vec3_t v3_normalize(vec3_t v)
nuclear@0: {
nuclear@0: 	vec3_t res;
nuclear@0: 	float len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
nuclear@0: 	if(len != 0.0f) {
nuclear@0: 		res.x = v.x / len;
nuclear@0: 		res.y = v.y / len;
nuclear@0: 		res.z = v.z / len;
nuclear@0: 	} else {
nuclear@0: 		res = v;
nuclear@0: 	}
nuclear@0: 	return res;
nuclear@0: }
nuclear@0: 
nuclear@0: static void init_noise()
nuclear@0: {
nuclear@0: 	int i;
nuclear@0: 
nuclear@0: 	/* calculate random gradients */
nuclear@0: 	for(i=0; i<B; i++) {
nuclear@0: 		perm[i] = i;	/* .. and initialize permutation mapping to identity */
nuclear@0: 
nuclear@0: 		grad1[i] = (float)((rand() % (B + B)) - B) / B;
nuclear@0: 
nuclear@0: 		grad2[i].x = (float)((rand() % (B + B)) - B) / B;
nuclear@0: 		grad2[i].y = (float)((rand() % (B + B)) - B) / B;
nuclear@0: 		grad2[i] = v2_normalize(grad2[i]);
nuclear@0: 
nuclear@0: 		grad3[i].x = (float)((rand() % (B + B)) - B) / B;
nuclear@0: 		grad3[i].y = (float)((rand() % (B + B)) - B) / B;
nuclear@0: 		grad3[i].z = (float)((rand() % (B + B)) - B) / B;
nuclear@0: 		grad3[i] = v3_normalize(grad3[i]);
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	/* permute indices by swapping them randomly */
nuclear@0: 	for(i=0; i<B; i++) {
nuclear@0: 		int rand_idx = rand() % B;
nuclear@0: 
nuclear@0: 		int tmp = perm[i];
nuclear@0: 		perm[i] = perm[rand_idx];
nuclear@0: 		perm[rand_idx] = tmp;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	/* fill up the rest of the arrays by duplicating the existing gradients */
nuclear@0: 	/* and permutations */
nuclear@0: 	for(i=0; i<B+2; i++) {
nuclear@0: 		perm[B + i] = perm[i];
nuclear@0: 		grad1[B + i] = grad1[i];
nuclear@0: 		grad2[B + i] = grad2[i];
nuclear@0: 		grad3[B + i] = grad3[i];
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: 
nuclear@0: float noise1(float x)
nuclear@0: {
nuclear@0: 	int bx0, bx1;
nuclear@0: 	float rx0, rx1, sx, u, v;
nuclear@0: 
nuclear@0: 	if(!tables_valid) {
nuclear@0: 		init_noise();
nuclear@0: 		tables_valid = 1;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	setup(x, bx0, bx1, rx0, rx1);
nuclear@0: 	sx = s_curve(rx0);
nuclear@0: 	u = rx0 * grad1[perm[bx0]];
nuclear@0: 	v = rx1 * grad1[perm[bx1]];
nuclear@0: 
nuclear@0: 	return lerp(u, v, sx);
nuclear@0: }
nuclear@0: 
nuclear@0: float fbm1(float x, int octaves)
nuclear@0: {
nuclear@0: 	int i;
nuclear@0: 	float res = 0.0f, freq = 1.0f;
nuclear@0: 	for(i=0; i<octaves; i++) {
nuclear@0: 		res += noise1(x * freq) / freq;
nuclear@0: 		freq *= 2.0f;
nuclear@0: 	}
nuclear@0: 	return res;
nuclear@0: }
nuclear@0: 
nuclear@0: float turbulence1(float x, int octaves)
nuclear@0: {
nuclear@0: 	int i;
nuclear@0: 	float res = 0.0f, freq = 1.0f;
nuclear@0: 	for(i=0; i<octaves; i++) {
nuclear@0: 		res += fabs(noise1(x * freq) / freq);
nuclear@0: 		freq *= 2.0f;
nuclear@0: 	}
nuclear@0: 	return res;
nuclear@0: }
nuclear@0: