| 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( seed: time(timer: 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( x: frequency, y: -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 | 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(x: 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(x: 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(v: g2[i]); |
| 187 | |
| 188 | for (j = 0 ; j < 3 ; j++) |
| 189 | g3[i][j] = (float)((rand() % (B + B)) - B) / B; |
| 190 | normalize3(v: 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 | } |
| 208 | |