1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/noise.c	Thu Jan 15 14:59:38 2015 +0200
     1.3 @@ -0,0 +1,146 @@
     1.4 +#include <stdlib.h>
     1.5 +#include <math.h>
     1.6 +#include "noise.h"
     1.7 +
     1.8 +typedef struct vec2 {
     1.9 +	float x, y;
    1.10 +} vec2_t;
    1.11 +
    1.12 +typedef struct vec3 {
    1.13 +	float x, y, z;
    1.14 +} vec3_t;
    1.15 +
    1.16 +/* ---- Ken Perlin's implementation of noise ---- */
    1.17 +
    1.18 +#define B	0x100
    1.19 +#define BM	0xff
    1.20 +#define N	0x1000
    1.21 +#define NP	12   /* 2^N */
    1.22 +#define NM	0xfff
    1.23 +
    1.24 +#define lerp(a, b, t)	((a) + ((b) - (a)) * t)
    1.25 +#define s_curve(t) (t * t * (3.0f - 2.0f * t))
    1.26 +
    1.27 +#define setup(elem, b0, b1, r0, r1) \
    1.28 +	do {							\
    1.29 +		float t = elem + N;		\
    1.30 +		b0 = ((int)t) & BM;			\
    1.31 +		b1 = (b0 + 1) & BM;			\
    1.32 +		r0 = t - (int)t;			\
    1.33 +		r1 = r0 - 1.0f;				\
    1.34 +	} while(0)
    1.35 +
    1.36 +
    1.37 +static int perm[B + B + 2];			/* permuted index from g_n onto themselves */
    1.38 +static vec3_t grad3[B + B + 2];		/* 3D random gradients */
    1.39 +static vec2_t grad2[B + B + 2];		/* 2D random gradients */
    1.40 +static float grad1[B + B + 2];	/* 1D random ... slopes */
    1.41 +static int tables_valid;
    1.42 +
    1.43 +static vec2_t v2_normalize(vec2_t v)
    1.44 +{
    1.45 +	vec2_t res;
    1.46 +	float len = sqrt(v.x * v.x + v.y * v.y);
    1.47 +	if(len != 0.0f) {
    1.48 +		res.x = v.x / len;
    1.49 +		res.y = v.y / len;
    1.50 +	} else {
    1.51 +		res = v;
    1.52 +	}
    1.53 +	return res;
    1.54 +}
    1.55 +
    1.56 +static vec3_t v3_normalize(vec3_t v)
    1.57 +{
    1.58 +	vec3_t res;
    1.59 +	float len = sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
    1.60 +	if(len != 0.0f) {
    1.61 +		res.x = v.x / len;
    1.62 +		res.y = v.y / len;
    1.63 +		res.z = v.z / len;
    1.64 +	} else {
    1.65 +		res = v;
    1.66 +	}
    1.67 +	return res;
    1.68 +}
    1.69 +
    1.70 +static void init_noise()
    1.71 +{
    1.72 +	int i;
    1.73 +
    1.74 +	/* calculate random gradients */
    1.75 +	for(i=0; i<B; i++) {
    1.76 +		perm[i] = i;	/* .. and initialize permutation mapping to identity */
    1.77 +
    1.78 +		grad1[i] = (float)((rand() % (B + B)) - B) / B;
    1.79 +
    1.80 +		grad2[i].x = (float)((rand() % (B + B)) - B) / B;
    1.81 +		grad2[i].y = (float)((rand() % (B + B)) - B) / B;
    1.82 +		grad2[i] = v2_normalize(grad2[i]);
    1.83 +
    1.84 +		grad3[i].x = (float)((rand() % (B + B)) - B) / B;
    1.85 +		grad3[i].y = (float)((rand() % (B + B)) - B) / B;
    1.86 +		grad3[i].z = (float)((rand() % (B + B)) - B) / B;
    1.87 +		grad3[i] = v3_normalize(grad3[i]);
    1.88 +	}
    1.89 +
    1.90 +	/* permute indices by swapping them randomly */
    1.91 +	for(i=0; i<B; i++) {
    1.92 +		int rand_idx = rand() % B;
    1.93 +
    1.94 +		int tmp = perm[i];
    1.95 +		perm[i] = perm[rand_idx];
    1.96 +		perm[rand_idx] = tmp;
    1.97 +	}
    1.98 +
    1.99 +	/* fill up the rest of the arrays by duplicating the existing gradients */
   1.100 +	/* and permutations */
   1.101 +	for(i=0; i<B+2; i++) {
   1.102 +		perm[B + i] = perm[i];
   1.103 +		grad1[B + i] = grad1[i];
   1.104 +		grad2[B + i] = grad2[i];
   1.105 +		grad3[B + i] = grad3[i];
   1.106 +	}
   1.107 +}
   1.108 +
   1.109 +
   1.110 +float noise1(float x)
   1.111 +{
   1.112 +	int bx0, bx1;
   1.113 +	float rx0, rx1, sx, u, v;
   1.114 +
   1.115 +	if(!tables_valid) {
   1.116 +		init_noise();
   1.117 +		tables_valid = 1;
   1.118 +	}
   1.119 +
   1.120 +	setup(x, bx0, bx1, rx0, rx1);
   1.121 +	sx = s_curve(rx0);
   1.122 +	u = rx0 * grad1[perm[bx0]];
   1.123 +	v = rx1 * grad1[perm[bx1]];
   1.124 +
   1.125 +	return lerp(u, v, sx);
   1.126 +}
   1.127 +
   1.128 +float fbm1(float x, int octaves)
   1.129 +{
   1.130 +	int i;
   1.131 +	float res = 0.0f, freq = 1.0f;
   1.132 +	for(i=0; i<octaves; i++) {
   1.133 +		res += noise1(x * freq) / freq;
   1.134 +		freq *= 2.0f;
   1.135 +	}
   1.136 +	return res;
   1.137 +}
   1.138 +
   1.139 +float turbulence1(float x, int octaves)
   1.140 +{
   1.141 +	int i;
   1.142 +	float res = 0.0f, freq = 1.0f;
   1.143 +	for(i=0; i<octaves; i++) {
   1.144 +		res += fabs(noise1(x * freq) / freq);
   1.145 +		freq *= 2.0f;
   1.146 +	}
   1.147 +	return res;
   1.148 +}
   1.149 +