dbf-udg: e09cbb2e9d4f libs/metasurf/metasurf.c

dbf-udg

view libs/metasurf/metasurf.c @ 5:e09cbb2e9d4f

metaballs

author	John Tsiombikas <nuclear@member.fsf.org>
date	Mon, 18 Feb 2013 03:46:52 +0200
parents
children

line source

1 /*

2 metasurf - a library for implicit surface polygonization

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 by

7 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 */

18 #include <stdio.h>

19 #include <stdlib.h>

20 #include "metasurf.h"

21 #include "mcubes.h"

23 #undef USE_MTETRA

24 #define USE_MCUBES

26 #if (defined(USE_MTETRA) && defined(USE_MCUBES)) || (!defined(USE_MTETRA) && !defined(USE_MCUBES))

27 #error "pick either USE_MTETRA or USE_MCUBES, not both..."

28 #endif

30 typedef float vec3[3];

32 struct metasurface {

33 vec3 min, max;

34 int res[3];

35 float thres;

37 msurf_eval_func_t eval;

38 msurf_vertex_func_t vertex;

39 msurf_normal_func_t normal;

41 float dx, dy, dz;

42 int flip;

44 vec3 vbuf[3];

45 int nverts;

46 };

48 static int msurf_init(struct metasurface *ms);

49 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz);

50 #ifdef USE_MTETRA

51 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val);

52 #endif

53 #ifdef USE_MCUBES

54 static void process_cube(struct metasurface *ms, vec3 *pos, float *val);

55 #endif

58 struct metasurface *msurf_create(void)

59 {

60 struct metasurface *ms;

62 if(!(ms = malloc(sizeof *ms))) {

63 return 0;

64 }

65 if(msurf_init(ms) == -1) {

66 free(ms);

67 }

68 return ms;

69 }

71 void msurf_free(struct metasurface *ms)

72 {

73 free(ms);

74 }

76 static int msurf_init(struct metasurface *ms)

77 {

78 ms->thres = 0.0;

79 ms->eval = 0;

80 ms->vertex = 0;

81 ms->normal = 0;

82 ms->min[0] = ms->min[1] = ms->min[2] = -1.0;

83 ms->max[0] = ms->max[1] = ms->max[2] = 1.0;

84 ms->res[0] = ms->res[1] = ms->res[2] = 40;

85 ms->nverts = 0;

87 ms->dx = ms->dy = ms->dz = 0.001;

88 ms->flip = 0;

90 return 0;

91 }

93 void msurf_inside(struct metasurface *ms, int inside)

94 {

95 switch(inside) {

96 case MSURF_GREATER:

97 ms->flip = 0;

98 break;

100 case MSURF_LESS:

101 ms->flip = 1;

102 break;

103

104 default:

105 fprintf(stderr, "msurf_inside expects MSURF_GREATER or MSURF_LESS\n");

106 }

107 }

108

109 void msurf_eval_func(struct metasurface *ms, msurf_eval_func_t func)

110 {

111 ms->eval = func;

112 }

113

114 void msurf_vertex_func(struct metasurface *ms, msurf_vertex_func_t func)

115 {

116 ms->vertex = func;

117 }

118

119 void msurf_normal_func(struct metasurface *ms, msurf_normal_func_t func)

120 {

121 ms->normal = func;

122 }

123

124 void msurf_bounds(struct metasurface *ms, float xmin, float ymin, float zmin, float xmax, float ymax, float zmax)

125 {

126 ms->min[0] = xmin;

127 ms->min[1] = ymin;

128 ms->min[2] = zmin;

129 ms->max[0] = xmax;

130 ms->max[1] = ymax;

131 ms->max[2] = zmax;

132 }

133

134 void msurf_resolution(struct metasurface *ms, int xres, int yres, int zres)

135 {

136 ms->res[0] = xres;

137 ms->res[1] = yres;

138 ms->res[2] = zres;

139 }

140

141 void msurf_threshold(struct metasurface *ms, float thres)

142 {

143 ms->thres = thres;

144 }

145

146

147 int msurf_polygonize(struct metasurface *ms)

148 {

149 int i, j, k;

150 vec3 pos, delta;

151

152 if(!ms->eval || !ms->vertex) {

153 fprintf(stderr, "you need to set eval and vertex callbacks before calling msurf_polygonize\n");

154 return -1;

155 }

156

157 for(i=0; i<3; i++) {

158 delta[i] = (ms->max[i] - ms->min[i]) / (float)ms->res[i];

159 }

160

161 pos[0] = ms->min[0];

162 for(i=0; i<ms->res[0] - 1; i++) {

163

164 pos[1] = ms->min[1];

165 for(j=0; j<ms->res[1] - 1; j++) {

166

167 pos[2] = ms->min[2];

168 for(k=0; k<ms->res[2] - 1; k++) {

169

170 process_cell(ms, pos, delta);

171

172 pos[2] += delta[2];

173 }

174 pos[1] += delta[1];

175 }

176 pos[0] += delta[0];

177 }

178 return 0;

179 }

180

181

182 static void process_cell(struct metasurface *ms, vec3 pos, vec3 sz)

183 {

184 int i;

185 vec3 p[8];

186 float val[8];

187

188 #ifdef USE_MTETRA

189 static int tetra[][4] = {

190 {0, 2, 3, 7},

191 {0, 2, 6, 7},

192 {0, 4, 6, 7},

193 {0, 6, 1, 2},

194 {0, 6, 1, 4},

195 {5, 6, 1, 4}

196 };

197 #endif

198

199 static const float offs[][3] = {

200 {0.0f, 0.0f, 0.0f},

201 {1.0f, 0.0f, 0.0f},

202 {1.0f, 1.0f, 0.0f},

203 {0.0f, 1.0f, 0.0f},

204 {0.0f, 0.0f, 1.0f},

205 {1.0f, 0.0f, 1.0f},

206 {1.0f, 1.0f, 1.0f},

207 {0.0f, 1.0f, 1.0f}

208 };

209

210 for(i=0; i<8; i++) {

211 p[i][0] = pos[0] + sz[0] * offs[i][2];

212 p[i][1] = pos[1] + sz[1] * offs[i][1];

213 p[i][2] = pos[2] + sz[2] * offs[i][0];

214

215 val[i] = ms->eval(p[i][0], p[i][1], p[i][2]);

216 }

217

218 #ifdef USE_MTETRA

219 for(i=0; i<6; i++) {

220 process_tetra(ms, tetra[i], p, val);

221 }

222 #endif

223 #ifdef USE_MCUBES

224 process_cube(ms, p, val);

225 #endif

226 }

227

228

229 /* ---- marching cubes implementation ---- */

230 #ifdef USE_MCUBES

231

232 static unsigned int mc_bitcode(float *val, float thres);

233

234 static void process_cube(struct metasurface *ms, vec3 *pos, float *val)

235 {

236 static const int pidx[12][2] = {

237 {0, 1}, {1, 2}, {2, 3}, {3, 0}, {4, 5}, {5, 6},

238 {6, 7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}

239 };

240 int i, j;

241 vec3 vert[12];

242 unsigned int code = mc_bitcode(val, ms->thres);

243

244 if(ms->flip) {

245 code = ~code & 0xff;

246 }

247

248 if(mc_edge_table[code] == 0) {

249 return;

250 }

251

252 for(i=0; i<12; i++) {

253 if(mc_edge_table[code] & (1 << i)) {

254 int p0 = pidx[i][0];

255 int p1 = pidx[i][1];

256

257 float t = (ms->thres - val[p0]) / (val[p1] - val[p0]);

258 vert[i][0] = pos[p0][0] + (pos[p1][0] - pos[p0][0]) * t;

259 vert[i][1] = pos[p0][1] + (pos[p1][1] - pos[p0][1]) * t;

260 vert[i][2] = pos[p0][2] + (pos[p1][2] - pos[p0][2]) * t;

261 }

262 }

263

264 for(i=0; mc_tri_table[code][i] != -1; i+=3) {

265 for(j=0; j<3; j++) {

266 float *v = vert[mc_tri_table[code][i + j]];

267

268 if(ms->normal) {

269 float dfdx, dfdy, dfdz;

270 dfdx = ms->eval(v[0] - ms->dx, v[1], v[2]) - ms->eval(v[0] + ms->dx, v[1], v[2]);

271 dfdy = ms->eval(v[0], v[1] - ms->dy, v[2]) - ms->eval(v[0], v[1] + ms->dy, v[2]);

272 dfdz = ms->eval(v[0], v[1], v[2] - ms->dz) - ms->eval(v[0], v[1], v[2] + ms->dz);

273

274 if(ms->flip) {

275 dfdx = -dfdx;

276 dfdy = -dfdy;

277 dfdz = -dfdz;

278 }

279 ms->normal(dfdx, dfdy, dfdz);

280 }

281

282 ms->vertex(v[0], v[1], v[2]);

283 }

284 }

285 }

286

287 static unsigned int mc_bitcode(float *val, float thres)

288 {

289 unsigned int i, res = 0;

290

291 for(i=0; i<8; i++) {

292 if(val[i] > thres) {

293 res |= 1 << i;

294 }

295 }

296 return res;

297 }

298 #endif /* USE_MCUBES */

299

300

301 /* ---- marching tetrahedra implementation (incomplete) ---- */

302 #ifdef USE_MTETRA

303

304 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres);

305 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)

306

307

308 #define REVBIT(x) ((x) & 8)

309 #define INV(x) (~(x) & 0xf)

310 #define EDGE(a, b) emmit(ms, val[idx[a]], val[idx[b]], pos[idx[a]], pos[idx[b]], REVBIT(code))

311 static void process_tetra(struct metasurface *ms, int *idx, vec3 *pos, float *val)

312 {

313 unsigned int code = mt_bitcode(val[idx[0]], val[idx[1]], val[idx[2]], val[idx[3]], ms->thres);

314

315 switch(code) {

316 case 1:

317 case INV(1):

318 EDGE(0, 1);

319 EDGE(0, 2);

320 EDGE(0, 3);

321 break;

322

323 case 2:

324 case INV(2):

325 EDGE(1, 0);

326 EDGE(1, 3);

327 EDGE(1, 2);

328 break;

329

330 case 3:

331 case INV(3):

332 EDGE(0, 3);

333 EDGE(0, 2);

334 EDGE(1, 3);

335

336 EDGE(1, 3);

337 EDGE(1, 2);

338 EDGE(0, 2);

339 break;

340

341 case 4:

342 case INV(4):

343 EDGE(2, 0);

344 EDGE(2, 1);

345 EDGE(2, 3);

346 break;

347

348 case 5:

349 case INV(5):

350 EDGE(0, 1);

351 EDGE(2, 3);

352 EDGE(0, 3);

353

354 EDGE(0, 1);

355 EDGE(1, 2);

356 EDGE(2, 3);

357 break;

358

359 case 6:

360 case INV(6):

361 EDGE(0, 1);

362 EDGE(1, 3);

363 EDGE(2, 3);

364

365 EDGE(0, 1);

366 EDGE(0, 2);

367 EDGE(2, 3);

368 break;

369

370 case 7:

371 case INV(7):

372 EDGE(3, 0);

373 EDGE(3, 2);

374 EDGE(3, 1);

375 break;

376

377 default:

378 break; /* cases 0 and 15 */

379 }

380 }

381

382 #define BIT(i) ((v##i > thres) ? (1 << i) : 0)

383 static unsigned int mt_bitcode(float v0, float v1, float v2, float v3, float thres)

384 {

385 return BIT(0) | BIT(1) | BIT(2) | BIT(3);

386 }

387

388 static void emmit(struct metasurface *ms, float v0, float v1, vec3 p0, vec3 p1, int rev)

389 {

390 int i;

391 float t = (ms->thres - v0) / (v1 - v0);

392

393 vec3 p;

394 for(i=0; i<3; i++) {

395 p[i] = p0[i] + (p1[i] - p0[i]) * t;

396 }

397 ms->vertex(p[0], p[1], p[2]);

398

399 /*for(i=0; i<3; i++) {

400 ms->vbuf[ms->nverts][i] = p0[i] + (p1[i] - p0[i]) * t;

401 }

402

403 if(++ms->nverts >= 3) {

404 ms->nverts = 0;

405

406 for(i=0; i<3; i++) {

407 int idx = rev ? (2 - i) : i;

408 ms->vertex(ms->vbuf[idx][0], ms->vbuf[idx][1], ms->vbuf[idx][2]);

409 }

410 }*/

411 }

412

413 #endif /* USE_MTETRA */