1 /* Coherent noise function over 1, 2 or 3 dimensions */
2 /* (copyright Ken Perlin) */
3 
4 #include "config.h"
5 #include <stdlib.h>
6 #include <stdio.h>
7 #include "perlin.h"
8 
9 /* random is not portable to all platforms */
10 #define random() g_rand_int (gr)
11 
12 static int p[B + B + 2];
13 static double g3[B + B + 2][3];
14 static double g2[B + B + 2][2];
15 static double g1[B + B + 2];
16 
17 double
noise1(double arg)18 noise1 (double arg)
19 {
20   int       bx0, bx1;
21   double    rx0, rx1, sx, t, u, v, vec[1];
22 
23   vec[0] = arg;
24 
25   setup (0, bx0, bx1, rx0, rx1);
26 
27   sx = s_curve (rx0);
28   u = rx0 * g1[p[bx0]];
29   v = rx1 * g1[p[bx1]];
30 
31   return (lerp (sx, u, v));
32 }
33 
34 double
noise2(double vec[2])35 noise2 (double vec[2])
36 {
37   int       bx0, bx1, by0, by1, b00, b10, b01, b11;
38   double    rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
39   int       i, j;
40 
41   setup (0, bx0, bx1, rx0, rx1);
42   setup (1, by0, by1, ry0, ry1);
43 
44   i = p[bx0];
45   j = p[bx1];
46 
47   b00 = p[i + by0];
48   b10 = p[j + by0];
49   b01 = p[i + by1];
50   b11 = p[j + by1];
51 
52   sx = s_curve (rx0);
53   sy = s_curve (ry0);
54 
55   q = g2[b00];
56   u = at2 (rx0, ry0);
57   q = g2[b10];
58   v = at2 (rx1, ry0);
59   a = lerp (sx, u, v);
60 
61   q = g2[b01];
62   u = at2 (rx0, ry1);
63   q = g2[b11];
64   v = at2 (rx1, ry1);
65   b = lerp (sx, u, v);
66 
67   return lerp (sy, a, b);
68 }
69 
70 double
noise3(double vec[3])71 noise3 (double vec[3])
72 {
73   int       bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
74   double    rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
75   int       i, j;
76 
77   setup (0, bx0, bx1, rx0, rx1);
78   setup (1, by0, by1, ry0, ry1);
79   setup (2, bz0, bz1, rz0, rz1);
80 
81   i = p[bx0];
82   j = p[bx1];
83 
84   b00 = p[i + by0];
85   b10 = p[j + by0];
86   b01 = p[i + by1];
87   b11 = p[j + by1];
88 
89   t = s_curve (rx0);
90   sy = s_curve (ry0);
91   sz = s_curve (rz0);
92 
93   q = g3[b00 + bz0];
94   u = at3 (rx0, ry0, rz0);
95   q = g3[b10 + bz0];
96   v = at3 (rx1, ry0, rz0);
97   a = lerp (t, u, v);
98 
99   q = g3[b01 + bz0];
100   u = at3 (rx0, ry1, rz0);
101   q = g3[b11 + bz0];
102   v = at3 (rx1, ry1, rz0);
103   b = lerp (t, u, v);
104 
105   c = lerp (sy, a, b);
106 
107   q = g3[b00 + bz1];
108   u = at3 (rx0, ry0, rz1);
109   q = g3[b10 + bz1];
110   v = at3 (rx1, ry0, rz1);
111   a = lerp (t, u, v);
112 
113   q = g3[b01 + bz1];
114   u = at3 (rx0, ry1, rz1);
115   q = g3[b11 + bz1];
116   v = at3 (rx1, ry1, rz1);
117   b = lerp (t, u, v);
118 
119   d = lerp (sy, a, b);
120 
121   return lerp (sz, c, d);
122 }
123 
124 void
normalize2(double v[2])125 normalize2 (double v[2])
126 {
127   double    s;
128 
129   s = sqrt (v[0] * v[0] + v[1] * v[1]);
130   v[0] = v[0] / s;
131   v[1] = v[1] / s;
132 }
133 
134 void
normalize3(double v[3])135 normalize3 (double v[3])
136 {
137   double    s;
138 
139   s = sqrt (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
140   v[0] = v[0] / s;
141   v[1] = v[1] / s;
142   v[2] = v[2] / s;
143 }
144 
145 void
perlin_init(void)146 perlin_init (void)
147 {
148   int       i, j, k;
149   static int initialized = 0;
150   if (initialized)
151     return;
152   /* this is racy - but we call it once when creating our perlin noise op */
153   GRand *gr = g_rand_new_with_seed (1234567890);
154 
155   for (i = 0; i < B; i++)
156     {
157       p[i] = i;
158       g1[i] = (double) ((random () % (B + B)) - B) / B;
159 
160       for (j = 0; j < 2; j++)
161         g2[i][j] = (double) ((random () % (B + B)) - B) / B;
162       normalize2 (g2[i]);
163 
164       for (j = 0; j < 3; j++)
165         g3[i][j] = (double) ((random () % (B + B)) - B) / B;
166       normalize3 (g3[i]);
167     }
168 
169   while (--i)
170     {
171       k = p[i];
172       p[i] = p[j = random () % B];
173       p[j] = k;
174     }
175 
176   for (i = 0; i < B + 2; i++)
177     {
178       p[B + i] = p[i];
179       g1[B + i] = g1[i];
180       for (j = 0; j < 2; j++)
181         g2[B + i][j] = g2[i][j];
182       for (j = 0; j < 3; j++)
183         g3[B + i][j] = g3[i][j];
184     }
185   initialized = 1;
186   g_rand_free (gr);
187 }
188 
189 /* --- My harmonic summing functions - PDB --------------------------*/
190 
191 /*
192    In what follows "alpha" is the weight when the sum is formed.
193    Typically it is 2, As this approaches 1 the function is noisier.
194    "beta" is the harmonic scaling/spacing, typically 2.
195 */
196 
197 double
PerlinNoise1D(double x,double alpha,double beta,int n)198 PerlinNoise1D (double x, double alpha, double beta, int n)
199 {
200   int       i;
201   double    val, sum = 0;
202   double    p, scale = 1;
203 
204   p = x;
205   for (i = 0; i < n; i++)
206     {
207       val = noise1 (p);
208       sum += val / scale;
209       scale *= alpha;
210       p *= beta;
211     }
212   return (sum);
213 }
214 
215 double
PerlinNoise2D(double x,double y,double alpha,double beta,int n)216 PerlinNoise2D (double x, double y, double alpha, double beta, int n)
217 {
218   int       i;
219   double    val, sum = 0;
220   double    p[2], scale = 1;
221 
222   p[0] = x;
223   p[1] = y;
224   for (i = 0; i < n; i++)
225     {
226       val = noise2 (p);
227       sum += val / scale;
228       scale *= alpha;
229       p[0] *= beta;
230       p[1] *= beta;
231     }
232   return (sum);
233 }
234 
235 double
PerlinNoise3D(double x,double y,double z,double alpha,double beta,int n)236 PerlinNoise3D (double x, double y, double z, double alpha, double beta, int n)
237 {
238   int       i;
239   double    val, sum = 0;
240   double    p[3], scale = 1.0;
241 
242   if (z < 0.0000)
243     return PerlinNoise2D (x, y, alpha, beta, n);
244   p[0] = x;
245   p[1] = y;
246   p[2] = z;
247 
248   for (i = 0; i < n; i++)
249     {
250       val = noise3 (p);
251       sum += val / scale;
252       scale *= alpha;
253       p[0] *= beta;
254       p[1] *= beta;
255       p[2] *= beta;
256     }
257   return (sum);
258 }
259