istereo2: 643f4ab609a4 libs/vmath/vmath.c

istereo2

view libs/vmath/vmath.c @ 35:643f4ab609a4

added readme and license

author	John Tsiombikas <nuclear@member.fsf.org>
date	Sat, 31 Oct 2015 05:45:35 +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 vec4_t tmp;

65 scalar_t tsq = t * t;

67 static mat4_t crspline_mat = {

68 {-1, 3, -3, 1},

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

70 {-1, 0, 1, 0},

71 {0, 2, 0, 0}

72 };

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

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

76 }

78 /** Bezier interpolation */

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

80 {

81 scalar_t omt, omt3, t3, f;

82 t3 = t * t * t;

83 omt = 1.0f - t;

84 omt3 = omt * omt * omt;

85 f = 3 * t * omt;

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

88 }

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

92 #define B 0x100

93 #define BM 0xff

94 #define N 0x1000

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

96 #define NM 0xfff

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

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

101 do { \

102 scalar_t t = elem + N; \

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

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

105 r0 = t - (int)t; \

106 r1 = r0 - 1.0f; \

107 } while(0)

108

109

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

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

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

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

114 static int tables_valid;

115

116 static void init_noise()

117 {

118 int i;

119

120 /* calculate random gradients */

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

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

123

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

125

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

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

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

129

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

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

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

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

134 }

135

136 /* permute indices by swapping them randomly */

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

138 int rand_idx = rand() % B;

139

140 int tmp = perm[i];

141 perm[i] = perm[rand_idx];

142 perm[rand_idx] = tmp;

143 }

144

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

146 /* and permutations */

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

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

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

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

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

152 }

153 }

154

155 scalar_t noise1(scalar_t x)

156 {

157 int bx0, bx1;

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

159

160 if(!tables_valid) {

161 init_noise();

162 tables_valid = 1;

163 }

164

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

166 sx = s_curve(rx0);

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

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

169

170 return lerp(u, v, sx);

171 }

172

173 scalar_t noise2(scalar_t x, scalar_t y)

174 {

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

176 int bx0, bx1, by0, by1;

177 scalar_t rx0, rx1, ry0, ry1;

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

179

180 if(!tables_valid) {

181 init_noise();

182 tables_valid = 1;

183 }

184

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

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

187

188 i = perm[bx0];

189 j = perm[bx1];

190

191 b00 = perm[i + by0];

192 b10 = perm[j + by0];

193 b01 = perm[i + by1];

194 b11 = perm[j + by1];

195

196 /* calculate hermite inteprolating factors */

197 sx = s_curve(rx0);

198 sy = s_curve(ry0);

199

200 /* interpolate along the left edge */

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

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

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

204

205 /* interpolate along the right edge */

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

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

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

209

210 /* interpolate between them */

211 return lerp(a, b, sy);

212 }

213

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

215 {

216 int i, j;

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

218 int b00, b10, b01, b11;

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

220 scalar_t sx, sy, sz;

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

222

223 if(!tables_valid) {

224 init_noise();

225 tables_valid = 1;

226 }

227

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

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

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

231

232 i = perm[bx0];

233 j = perm[bx1];

234

235 b00 = perm[i + by0];

236 b10 = perm[j + by0];

237 b01 = perm[i + by1];

238 b11 = perm[j + by1];

239

240 /* calculate hermite interpolating factors */

241 sx = s_curve(rx0);

242 sy = s_curve(ry0);

243 sz = s_curve(rz0);

244

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

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

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

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

249

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

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

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

253

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

255

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

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

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

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

260

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

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

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

264

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

266

267 /* interpolate between slices */

268 return lerp(c, d, sz);

269 }

270

271 scalar_t fbm1(scalar_t x, int octaves)

272 {

273 int i;

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

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

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

277 freq *= 2.0f;

278 }

279 return res;

280 }

281

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

283 {

284 int i;

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

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

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

288 freq *= 2.0f;

289 }

290 return res;

291 }

292

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

294 {

295 int i;

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

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

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

299 freq *= 2.0f;

300 }

301 return res;

302 }

303

304 scalar_t turbulence1(scalar_t x, int octaves)

305 {

306 int i;

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

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

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

310 freq *= 2.0f;

311 }

312 return res;

313 }

314

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

316 {

317 int i;

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

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

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

321 freq *= 2.0f;

322 }

323 return res;

324 }

325

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

327 {

328 int i;

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

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

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

332 freq *= 2.0f;

333 }

334 return res;

335 }