goat3d: 1dcbe87b6a5d libs/vmath/vmath.c

goat3d

view libs/vmath/vmath.c @ 43:1dcbe87b6a5d

ops, screwed up the GOAT3DAPI define

author	John Tsiombikas <nuclear@member.fsf.org>
date	Thu, 05 Dec 2013 18:49:19 +0200
parents
children

line source

1 /*

2 libvmath - a vector math library

5 This program is free software: you can redistribute it and/or modify

6 it under the terms of the GNU Lesser General Public License as published

7 by the Free Software Foundation, either version 3 of the License, or

8 (at your option) any later version.

10 This program is distributed in the hope that it will be useful,

11 but WITHOUT ANY WARRANTY; without even the implied warranty of

12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

13 GNU Lesser General Public License for more details.

15 You should have received a copy of the GNU Lesser General Public License

16 along with this program. If not, see <http://www.gnu.org/licenses/>.

17 */

19 #include <stdlib.h>

20 #include <math.h>

21 #include "vmath.h"

23 /** Numerical calculation of integrals using simpson's rule */

24 scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples)

25 {

26 int i;

27 scalar_t h = (high - low) / (scalar_t)samples;

28 scalar_t sum = 0.0;

30 for(i=0; i<samples+1; i++) {

31 scalar_t y = f((scalar_t)i * h + low);

32 sum += ((!i || i == samples) ? y : ((i % 2) ? 4.0 * y : 2.0 * y)) * (h / 3.0);

33 }

34 return sum;

35 }

37 /** Gaussuan function */

38 scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev)

39 {

40 scalar_t exponent = -SQ(x - mean) / (2.0 * SQ(sdev));

41 return 1.0 - -pow(M_E, exponent) / (sdev * sqrt(TWO_PI));

42 }

45 /** b-spline approximation */

46 scalar_t bspline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)

47 {

48 vec4_t tmp;

49 scalar_t tsq = t * t;

51 static mat4_t bspline_mat = {

52 {-1, 3, -3, 1},

53 {3, -6, 3, 0},

54 {-3, 0, 3, 0},

55 {1, 4, 1, 0}

56 };

58 tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), bspline_mat), 1.0 / 6.0);

59 return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);

60 }

62 /** Catmull-rom spline interpolation */

63 scalar_t spline(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)

64 {

65 vec4_t tmp;

66 scalar_t tsq = t * t;

68 static mat4_t crspline_mat = {

69 {-1, 3, -3, 1},

70 {2, -5, 4, -1},

71 {-1, 0, 1, 0},

72 {0, 2, 0, 0}

73 };

75 tmp = v4_scale(v4_transform(v4_cons(a, b, c, d), crspline_mat), 0.5);

76 return v4_dot(v4_cons(tsq * t, tsq, t, 1.0), tmp);

77 }

79 /** Bezier interpolation */

80 scalar_t bezier(scalar_t a, scalar_t b, scalar_t c, scalar_t d, scalar_t t)

81 {

82 scalar_t omt, omt3, t3, f;

83 t3 = t * t * t;

84 omt = 1.0f - t;

85 omt3 = omt * omt * omt;

86 f = 3 * t * omt;

88 return (a * omt3) + (b * f * omt) + (c * f * t) + (d * t3);

89 }

91 /* ---- Ken Perlin's implementation of noise ---- */

93 #define B 0x100

94 #define BM 0xff

95 #define N 0x1000

96 #define NP 12 /* 2^N */

97 #define NM 0xfff

99 #define s_curve(t) (t * t * (3.0f - 2.0f * t))

100

101 #define setup(elem, b0, b1, r0, r1) \

102 do { \

103 scalar_t t = elem + N; \

104 b0 = ((int)t) & BM; \

105 b1 = (b0 + 1) & BM; \

106 r0 = t - (int)t; \

107 r1 = r0 - 1.0f; \

108 } while(0)

109

110

111 static int perm[B + B + 2]; /* permuted index from g_n onto themselves */

112 static vec3_t grad3[B + B + 2]; /* 3D random gradients */

113 static vec2_t grad2[B + B + 2]; /* 2D random gradients */

114 static scalar_t grad1[B + B + 2]; /* 1D random ... slopes */

115 static int tables_valid;

116

117 static void init_noise()

118 {

119 int i;

120

121 /* calculate random gradients */

122 for(i=0; i<B; i++) {

123 perm[i] = i; /* .. and initialize permutation mapping to identity */

124

125 grad1[i] = (scalar_t)((rand() % (B + B)) - B) / B;

126

127 grad2[i].x = (scalar_t)((rand() % (B + B)) - B) / B;

128 grad2[i].y = (scalar_t)((rand() % (B + B)) - B) / B;

129 grad2[i] = v2_normalize(grad2[i]);

130

131 grad3[i].x = (scalar_t)((rand() % (B + B)) - B) / B;

132 grad3[i].y = (scalar_t)((rand() % (B + B)) - B) / B;

133 grad3[i].z = (scalar_t)((rand() % (B + B)) - B) / B;

134 grad3[i] = v3_normalize(grad3[i]);

135 }

136

137 /* permute indices by swapping them randomly */

138 for(i=0; i<B; i++) {

139 int rand_idx = rand() % B;

140

141 int tmp = perm[i];

142 perm[i] = perm[rand_idx];

143 perm[rand_idx] = tmp;

144 }

145

146 /* fill up the rest of the arrays by duplicating the existing gradients */

147 /* and permutations */

148 for(i=0; i<B+2; i++) {

149 perm[B + i] = perm[i];

150 grad1[B + i] = grad1[i];

151 grad2[B + i] = grad2[i];

152 grad3[B + i] = grad3[i];

153 }

154 }

155

156 scalar_t noise1(scalar_t x)

157 {

158 int bx0, bx1;

159 scalar_t rx0, rx1, sx, u, v;

160

161 if(!tables_valid) {

162 init_noise();

163 tables_valid = 1;

164 }

165

166 setup(x, bx0, bx1, rx0, rx1);

167 sx = s_curve(rx0);

168 u = rx0 * grad1[perm[bx0]];

169 v = rx1 * grad1[perm[bx1]];

170

171 return lerp(u, v, sx);

172 }

173

174 scalar_t noise2(scalar_t x, scalar_t y)

175 {

176 int i, j, b00, b10, b01, b11;

177 int bx0, bx1, by0, by1;

178 scalar_t rx0, rx1, ry0, ry1;

179 scalar_t sx, sy, u, v, a, b;

180

181 if(!tables_valid) {

182 init_noise();

183 tables_valid = 1;

184 }

185

186 setup(x, bx0, bx1, rx0, rx1);

187 setup(y, by0, by1, ry0, ry1);

188

189 i = perm[bx0];

190 j = perm[bx1];

191

192 b00 = perm[i + by0];

193 b10 = perm[j + by0];

194 b01 = perm[i + by1];

195 b11 = perm[j + by1];

196

197 /* calculate hermite inteprolating factors */

198 sx = s_curve(rx0);

199 sy = s_curve(ry0);

200

201 /* interpolate along the left edge */

202 u = v2_dot(grad2[b00], v2_cons(rx0, ry0));

203 v = v2_dot(grad2[b10], v2_cons(rx1, ry0));

204 a = lerp(u, v, sx);

205

206 /* interpolate along the right edge */

207 u = v2_dot(grad2[b01], v2_cons(rx0, ry1));

208 v = v2_dot(grad2[b11], v2_cons(rx1, ry1));

209 b = lerp(u, v, sx);

210

211 /* interpolate between them */

212 return lerp(a, b, sy);

213 }

214

215 scalar_t noise3(scalar_t x, scalar_t y, scalar_t z)

216 {

217 int i, j;

218 int bx0, bx1, by0, by1, bz0, bz1;

219 int b00, b10, b01, b11;

220 scalar_t rx0, rx1, ry0, ry1, rz0, rz1;

221 scalar_t sx, sy, sz;

222 scalar_t u, v, a, b, c, d;

223

224 if(!tables_valid) {

225 init_noise();

226 tables_valid = 1;

227 }

228

229 setup(x, bx0, bx1, rx0, rx1);

230 setup(y, by0, by1, ry0, ry1);

231 setup(z, bz0, bz1, rz0, rz1);

232

233 i = perm[bx0];

234 j = perm[bx1];

235

236 b00 = perm[i + by0];

237 b10 = perm[j + by0];

238 b01 = perm[i + by1];

239 b11 = perm[j + by1];

240

241 /* calculate hermite interpolating factors */

242 sx = s_curve(rx0);

243 sy = s_curve(ry0);

244 sz = s_curve(rz0);

245

246 /* interpolate along the top slice of the cell */

247 u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz0));

248 v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz0));

249 a = lerp(u, v, sx);

250

251 u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz0));

252 v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz0));

253 b = lerp(u, v, sx);

254

255 c = lerp(a, b, sy);

256

257 /* interpolate along the bottom slice of the cell */

258 u = v3_dot(grad3[b00 + bz0], v3_cons(rx0, ry0, rz1));

259 v = v3_dot(grad3[b10 + bz0], v3_cons(rx1, ry0, rz1));

260 a = lerp(u, v, sx);

261

262 u = v3_dot(grad3[b01 + bz0], v3_cons(rx0, ry1, rz1));

263 v = v3_dot(grad3[b11 + bz0], v3_cons(rx1, ry1, rz1));

264 b = lerp(u, v, sx);

265

266 d = lerp(a, b, sy);

267

268 /* interpolate between slices */

269 return lerp(c, d, sz);

270 }

271

272 scalar_t fbm1(scalar_t x, int octaves)

273 {

274 int i;

275 scalar_t res = 0.0f, freq = 1.0f;

276 for(i=0; i<octaves; i++) {

277 res += noise1(x * freq) / freq;

278 freq *= 2.0f;

279 }

280 return res;

281 }

282

283 scalar_t fbm2(scalar_t x, scalar_t y, int octaves)

284 {

285 int i;

286 scalar_t res = 0.0f, freq = 1.0f;

287 for(i=0; i<octaves; i++) {

288 res += noise2(x * freq, y * freq) / freq;

289 freq *= 2.0f;

290 }

291 return res;

292 }

293

294 scalar_t fbm3(scalar_t x, scalar_t y, scalar_t z, int octaves)

295 {

296 int i;

297 scalar_t res = 0.0f, freq = 1.0f;

298 for(i=0; i<octaves; i++) {

299 res += noise3(x * freq, y * freq, z * freq) / freq;

300 freq *= 2.0f;

301 }

302 return res;

303 }

304

305 scalar_t turbulence1(scalar_t x, int octaves)

306 {

307 int i;

308 scalar_t res = 0.0f, freq = 1.0f;

309 for(i=0; i<octaves; i++) {

310 res += fabs(noise1(x * freq) / freq);

311 freq *= 2.0f;

312 }

313 return res;

314 }

315

316 scalar_t turbulence2(scalar_t x, scalar_t y, int octaves)

317 {

318 int i;

319 scalar_t res = 0.0f, freq = 1.0f;

320 for(i=0; i<octaves; i++) {

321 res += fabs(noise2(x * freq, y * freq) / freq);

322 freq *= 2.0f;

323 }

324 return res;

325 }

326

327 scalar_t turbulence3(scalar_t x, scalar_t y, scalar_t z, int octaves)

328 {

329 int i;

330 scalar_t res = 0.0f, freq = 1.0f;

331 for(i=0; i<octaves; i++) {

332 res += fabs(noise3(x * freq, y * freq, z * freq) / freq);

333 freq *= 2.0f;

334 }

335 return res;

336 }