|
1 /***************************************************************** |
|
2 |
|
3 Implementation of the fractional Brownian motion algorithm. These |
|
4 functions were originally the work of F. Kenton Musgrave. |
|
5 For documentation of the different functions please refer to the |
|
6 book: |
|
7 "Texturing and modeling: a procedural approach" |
|
8 by David S. Ebert et. al. |
|
9 |
|
10 ******************************************************************/ |
|
11 |
|
12 #if defined (_MSC_VER) |
|
13 #include <qglobal.h> |
|
14 #endif |
|
15 |
|
16 #include <time.h> |
|
17 #include <stdlib.h> |
|
18 #include "fbm.h" |
|
19 |
|
20 #if defined(Q_CC_MSVC) |
|
21 #pragma warning(disable:4244) |
|
22 #endif |
|
23 |
|
24 /* Definitions used by the noise2() functions */ |
|
25 |
|
26 //#define B 0x100 |
|
27 //#define BM 0xff |
|
28 #define B 0x20 |
|
29 #define BM 0x1f |
|
30 |
|
31 #define N 0x1000 |
|
32 #define NP 12 /* 2^N */ |
|
33 #define NM 0xfff |
|
34 |
|
35 static int p[B + B + 2]; |
|
36 static float g3[B + B + 2][3]; |
|
37 static float g2[B + B + 2][2]; |
|
38 static float g1[B + B + 2]; |
|
39 static int start = 1; |
|
40 |
|
41 static void init(void); |
|
42 |
|
43 #define s_curve(t) ( t * t * (3. - 2. * t) ) |
|
44 |
|
45 #define lerp(t, a, b) ( a + t * (b - a) ) |
|
46 |
|
47 #define setup(i,b0,b1,r0,r1)\ |
|
48 t = vec[i] + N;\ |
|
49 b0 = ((int)t) & BM;\ |
|
50 b1 = (b0+1) & BM;\ |
|
51 r0 = t - (int)t;\ |
|
52 r1 = r0 - 1.; |
|
53 #define at3(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] ) |
|
54 |
|
55 /* Fractional Brownian Motion function */ |
|
56 |
|
57 double fBm( Vector point, double H, double lacunarity, double octaves, |
|
58 int init ) |
|
59 { |
|
60 |
|
61 double value, frequency, remainder; |
|
62 int i; |
|
63 static double exponent_array[10]; |
|
64 float vec[3]; |
|
65 |
|
66 /* precompute and store spectral weights */ |
|
67 if ( init ) { |
|
68 start = 1; |
|
69 srand( time(0) ); |
|
70 /* seize required memory for exponent_array */ |
|
71 frequency = 1.0; |
|
72 for (i=0; i<=octaves; i++) { |
|
73 /* compute weight for each frequency */ |
|
74 exponent_array[i] = pow( frequency, -H ); |
|
75 frequency *= lacunarity; |
|
76 } |
|
77 } |
|
78 |
|
79 value = 0.0; /* initialize vars to proper values */ |
|
80 frequency = 1.0; |
|
81 vec[0]=point.x; |
|
82 vec[1]=point.y; |
|
83 vec[2]=point.z; |
|
84 |
|
85 |
|
86 /* inner loop of spectral construction */ |
|
87 for (i=0; i<octaves; i++) { |
|
88 /* value += noise3( vec ) * exponent_array[i];*/ |
|
89 value += noise3( vec ) * exponent_array[i]; |
|
90 vec[0] *= lacunarity; |
|
91 vec[1] *= lacunarity; |
|
92 vec[2] *= lacunarity; |
|
93 } /* for */ |
|
94 |
|
95 remainder = octaves - (int)octaves; |
|
96 if ( remainder ) /* add in ``octaves'' remainder */ |
|
97 /* ``i'' and spatial freq. are preset in loop above */ |
|
98 value += remainder * noise3( vec ) * exponent_array[i]; |
|
99 |
|
100 return( value ); |
|
101 |
|
102 } /* fBm() */ |
|
103 |
|
104 |
|
105 float noise3(float vec[3]) |
|
106 { |
|
107 int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11; |
|
108 float rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v; |
|
109 register int i, j; |
|
110 |
|
111 if (start) { |
|
112 start = 0; |
|
113 init(); |
|
114 } |
|
115 |
|
116 setup(0, bx0,bx1, rx0,rx1); |
|
117 setup(1, by0,by1, ry0,ry1); |
|
118 setup(2, bz0,bz1, rz0,rz1); |
|
119 |
|
120 i = p[ bx0 ]; |
|
121 j = p[ bx1 ]; |
|
122 |
|
123 b00 = p[ i + by0 ]; |
|
124 b10 = p[ j + by0 ]; |
|
125 b01 = p[ i + by1 ]; |
|
126 b11 = p[ j + by1 ]; |
|
127 |
|
128 t = s_curve(rx0); |
|
129 sy = s_curve(ry0); |
|
130 sz = s_curve(rz0); |
|
131 |
|
132 |
|
133 q = g3[ b00 + bz0 ] ; u = at3(rx0,ry0,rz0); |
|
134 q = g3[ b10 + bz0 ] ; v = at3(rx1,ry0,rz0); |
|
135 a = lerp(t, u, v); |
|
136 |
|
137 q = g3[ b01 + bz0 ] ; u = at3(rx0,ry1,rz0); |
|
138 q = g3[ b11 + bz0 ] ; v = at3(rx1,ry1,rz0); |
|
139 b = lerp(t, u, v); |
|
140 |
|
141 c = lerp(sy, a, b); |
|
142 |
|
143 q = g3[ b00 + bz1 ] ; u = at3(rx0,ry0,rz1); |
|
144 q = g3[ b10 + bz1 ] ; v = at3(rx1,ry0,rz1); |
|
145 a = lerp(t, u, v); |
|
146 |
|
147 q = g3[ b01 + bz1 ] ; u = at3(rx0,ry1,rz1); |
|
148 q = g3[ b11 + bz1 ] ; v = at3(rx1,ry1,rz1); |
|
149 b = lerp(t, u, v); |
|
150 |
|
151 d = lerp(sy, a, b); |
|
152 |
|
153 return lerp(sz, c, d); |
|
154 } |
|
155 |
|
156 static void normalize2(float v[2]) |
|
157 { |
|
158 float s; |
|
159 |
|
160 s = sqrt(v[0] * v[0] + v[1] * v[1]); |
|
161 v[0] = v[0] / s; |
|
162 v[1] = v[1] / s; |
|
163 } |
|
164 |
|
165 static void normalize3(float v[3]) |
|
166 { |
|
167 float s; |
|
168 |
|
169 s = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); |
|
170 v[0] = v[0] / s; |
|
171 v[1] = v[1] / s; |
|
172 v[2] = v[2] / s; |
|
173 } |
|
174 |
|
175 static void init(void) |
|
176 { |
|
177 int i, j, k; |
|
178 |
|
179 for (i = 0 ; i < B ; i++) { |
|
180 p[i] = i; |
|
181 |
|
182 g1[i] = (float)((rand() % (B + B)) - B) / B; |
|
183 |
|
184 for (j = 0 ; j < 2 ; j++) |
|
185 g2[i][j] = (float)((rand() % (B + B)) - B) / B; |
|
186 normalize2(g2[i]); |
|
187 |
|
188 for (j = 0 ; j < 3 ; j++) |
|
189 g3[i][j] = (float)((rand() % (B + B)) - B) / B; |
|
190 normalize3(g3[i]); |
|
191 } |
|
192 |
|
193 while (--i) { |
|
194 k = p[i]; |
|
195 p[i] = p[j = rand() % B]; |
|
196 p[j] = k; |
|
197 } |
|
198 |
|
199 for (i = 0 ; i < B + 2 ; i++) { |
|
200 p[B + i] = p[i]; |
|
201 g1[B + i] = g1[i]; |
|
202 for (j = 0 ; j < 2 ; j++) |
|
203 g2[B + i][j] = g2[i][j]; |
|
204 for (j = 0 ; j < 3 ; j++) |
|
205 g3[B + i][j] = g3[i][j]; |
|
206 } |
|
207 } |