nuclear@0: #include <stdlib.h>
nuclear@0: #include "pnoise.h"
nuclear@0: #include "vmath/vmath.h"
nuclear@0: 
nuclear@0: #define B	0x100
nuclear@0: #define BM	0xff
nuclear@0: #define N	0x1000
nuclear@0: #define NP	12
nuclear@0: #define NM	0xfff
nuclear@0: 
nuclear@0: #define s_curve(t)	((t) * (t) * (3.0f - 2.0f * (t)))
nuclear@0: #define setup(elem, b0, b1, r0, r1) \
nuclear@0: 	do {							\
nuclear@0: 		scalar_t 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: #define setup_p(elem, b0, b1, r0, r1, p) \
nuclear@0: 	do {							\
nuclear@0: 		scalar_t t = elem + N;		\
nuclear@0: 		b0 = (((int)t) & BM) % p;	\
nuclear@0: 		b1 = ((b0 + 1) & BM) % p;	\
nuclear@0: 		r0 = t - (int)t;			\
nuclear@0: 		r1 = r0 - 1.0f;				\
nuclear@0: 	} while(0)
nuclear@0: 
nuclear@0: static int perm[B + B + 2];
nuclear@0: static vec2_t grad2[B + B + 2];
nuclear@0: static bool tables_valid;
nuclear@0: 
nuclear@0: static void init_noise()
nuclear@0: {
nuclear@0: 	for(int i=0; i<B; i++) {
nuclear@0: 		perm[i] = i;
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: 
nuclear@0: 	for(int 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: 	for(int i=0; i<B+2; i++) {
nuclear@0: 		perm[B + i] = perm[i];
nuclear@0: 		grad2[B + i] = grad2[i];
nuclear@0: 	}
nuclear@0: }
nuclear@0: 
nuclear@0: #define lerp(a, b, t)	((a) + ((b) - (a)) * t)
nuclear@0: 
nuclear@0: float dbg_noise2(float x, float y)
nuclear@0: {
nuclear@0: 	if(!tables_valid) {
nuclear@0: 		init_noise();
nuclear@0: 		tables_valid = true;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	int bx0, bx1, by0, by1;
nuclear@0: 	float rx0, rx1, ry0, ry1;
nuclear@0: 	setup(x, bx0, bx1, rx0, rx1);
nuclear@0: 	setup(y, by0, by1, ry0, ry1);
nuclear@0: 
nuclear@0: 	int i = perm[bx0];
nuclear@0: 	int j = perm[bx1];
nuclear@0: 
nuclear@0: 	int b00 = perm[i + by0];
nuclear@0: 	int b10 = perm[j + by0];
nuclear@0: 	int b01 = perm[i + by1];
nuclear@0: 	int b11 = perm[j + by1];
nuclear@0: 
nuclear@0: 	float sx = s_curve(rx0);
nuclear@0: 	float sy = s_curve(ry0);
nuclear@0: 
nuclear@0: 	vec2_t g00 = grad2[b00];
nuclear@0: 	vec2_t g10 = grad2[b10];
nuclear@0: 	vec2_t g01 = grad2[b01];
nuclear@0: 	vec2_t g11 = grad2[b11];
nuclear@0: 
nuclear@0: 	float u = g00.x * rx0 + g00.y * ry0;
nuclear@0: 	float v = g10.x * rx1 + g10.y * ry0;
nuclear@0: 	float a = lerp(u, v, sx);
nuclear@0: 
nuclear@0: 	u = g01.x * rx0 + g01.y * ry1;
nuclear@0: 	v = g11.x * rx1 + g11.y * ry1;
nuclear@0: 	float b = lerp(u, v, sx);
nuclear@0: 
nuclear@0: 	return lerp(a, b, sy);
nuclear@0: }
nuclear@0: 
nuclear@0: float pnoise2(float x, float y, int period)
nuclear@0: {
nuclear@0: 	if(!tables_valid) {
nuclear@0: 		init_noise();
nuclear@0: 		tables_valid = true;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	int bx0, bx1, by0, by1;
nuclear@0: 	float rx0, rx1, ry0, ry1;
nuclear@0: 	setup_p(x, bx0, bx1, rx0, rx1, period);
nuclear@0: 	setup_p(y, by0, by1, ry0, ry1, period);
nuclear@0: 
nuclear@0: 	int i = perm[bx0];
nuclear@0: 	int j = perm[bx1];
nuclear@0: 
nuclear@0: 	int b00 = perm[i + by0];
nuclear@0: 	int b10 = perm[j + by0];
nuclear@0: 	int b01 = perm[i + by1];
nuclear@0: 	int b11 = perm[j + by1];
nuclear@0: 
nuclear@0: 	float sx = s_curve(rx0);
nuclear@0: 	float sy = s_curve(ry0);
nuclear@0: 
nuclear@0: 	vec2_t g00 = grad2[b00];
nuclear@0: 	vec2_t g10 = grad2[b10];
nuclear@0: 	vec2_t g01 = grad2[b01];
nuclear@0: 	vec2_t g11 = grad2[b11];
nuclear@0: 
nuclear@0: 	float u = g00.x * rx0 + g00.y * ry0;
nuclear@0: 	float v = g10.x * rx1 + g10.y * ry0;
nuclear@0: 	float a = lerp(u, v, sx);
nuclear@0: 
nuclear@0: 	u = g01.x * rx0 + g01.y * ry1;
nuclear@0: 	v = g11.x * rx1 + g11.y * ry1;
nuclear@0: 	float b = lerp(u, v, sx);
nuclear@0: 
nuclear@0: 	return lerp(a, b, sy);
nuclear@0: }
nuclear@0: 
nuclear@0: float pfbm2(float x, float y, int octaves, int period)
nuclear@0: {
nuclear@0: 	float res = 0.0;
nuclear@0: 	float s = 1.0;
nuclear@0: 
nuclear@0: 	for(int i=0; i<octaves; i++) {
nuclear@0: 		float v = pnoise2(x * s, y * s, period);
nuclear@0: 		res += v / s;
nuclear@0: 		s *= 2.0;
nuclear@0: 		period *= 2;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	return res;
nuclear@0: }
nuclear@0: 
nuclear@0: float pturbulence2(float x, float y, int octaves, int period)
nuclear@0: {
nuclear@0: 	float res = 0.0;
nuclear@0: 	float s = 1.0;
nuclear@0: 
nuclear@0: 	for(int i=0; i<octaves; i++) {
nuclear@0: 		float v = fabs(pnoise2(x * s, y * s, period));
nuclear@0: 		res += v / s;
nuclear@0: 		s *= 2.0;
nuclear@0: 		period *= 2;
nuclear@0: 	}
nuclear@0: 
nuclear@0: 	return res;
nuclear@0: }