metasurf: 52664d3451ad src/metasurf.c

metasurf

view src/metasurf.c @ 3:52664d3451ad

added volume rendering example

author	John Tsiombikas <nuclear@member.fsf.org>
date	Tue, 25 Oct 2011 13:30:03 +0300
parents	7aa4627e492b
children	2c575855f707

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 vec3 vbuf[3];

42 int nverts;

43 };

45 static int msurf_init(struct metasurface *ms);

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

47 #ifdef USE_MTETRA

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

49 #endif

50 #ifdef USE_MCUBES

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

52 #endif

55 struct metasurface *msurf_create(void)

56 {

57 struct metasurface *ms;

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

60 return 0;

61 }

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

63 free(ms);

64 }

65 return ms;

66 }

68 void msurf_free(struct metasurface *ms)

69 {

70 free(ms);

71 }

73 static int msurf_init(struct metasurface *ms)

74 {

75 ms->thres = 0.0;

76 ms->eval = 0;

77 ms->vertex = 0;

78 ms->normal = 0;

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

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

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

82 ms->nverts = 0;

84 return 0;

85 }

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

88 {

89 ms->eval = func;

90 }

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

93 {

94 ms->vertex = func;

95 }

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

98 {

99 ms->normal = func;

100 }

101

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

103 {

104 ms->min[0] = xmin;

105 ms->min[1] = ymin;

106 ms->min[2] = zmin;

107 ms->max[0] = xmax;

108 ms->max[1] = ymax;

109 ms->max[2] = zmax;

110 }

111

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

113 {

114 ms->res[0] = xres;

115 ms->res[1] = yres;

116 ms->res[2] = zres;

117 }

118

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

120 {

121 ms->thres = thres;

122 }

123

124

125 int msurf_polygonize(struct metasurface *ms)

126 {

127 int i, j, k;

128 vec3 pos, delta;

129

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

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

132 return -1;

133 }

134

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

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

137 }

138

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

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

141

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

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

144

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

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

147

148 process_cell(ms, pos, delta);

149

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

151 }

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

153 }

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

155 }

156 return 0;

157 }

158

159

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

161 {

162 int i;

163 vec3 p[8];

164 float val[8];

165

166 #ifdef USE_MTETRA

167 static int tetra[][4] = {

168 {0, 2, 3, 7},

169 {0, 2, 6, 7},

170 {0, 4, 6, 7},

171 {0, 6, 1, 2},

172 {0, 6, 1, 4},

173 {5, 6, 1, 4}

174 };

175 #endif

176

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

178 {0.0f, 0.0f, 0.0f},

179 {1.0f, 0.0f, 0.0f},

180 {1.0f, 1.0f, 0.0f},

181 {0.0f, 1.0f, 0.0f},

182 {0.0f, 0.0f, 1.0f},

183 {1.0f, 0.0f, 1.0f},

184 {1.0f, 1.0f, 1.0f},

185 {0.0f, 1.0f, 1.0f}

186 };

187

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

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

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

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

192

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

194 }

195

196 #ifdef USE_MTETRA

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

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

199 }

200 #endif

201 #ifdef USE_MCUBES

202 process_cube(ms, p, val);

203 #endif

204 }

205

206

207 /* ---- marching cubes implementation ---- */

208 #ifdef USE_MCUBES

209

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

211

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

213 {

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

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

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

217 };

218 int i, j;

219 vec3 vert[12];

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

221

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

223 return;

224 }

225

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

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

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

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

230

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

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

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

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

235 }

236 }

237

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

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

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

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

242 }

243 }

244 }

245

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

247 {

248 unsigned int i, res = 0;

249

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

251 if(val[i] > thres) {

252 res |= 1 << i;

253 }

254 }

255 return res;

256 }

257 #endif /* USE_MCUBES */

258

259

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

261 #ifdef USE_MTETRA

262

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

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

265

266

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

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

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

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

271 {

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

273

274 switch(code) {

275 case 1:

276 case INV(1):

277 EDGE(0, 1);

278 EDGE(0, 2);

279 EDGE(0, 3);

280 break;

281

282 case 2:

283 case INV(2):

284 EDGE(1, 0);

285 EDGE(1, 3);

286 EDGE(1, 2);

287 break;

288

289 case 3:

290 case INV(3):

291 EDGE(0, 3);

292 EDGE(0, 2);

293 EDGE(1, 3);

294

295 EDGE(1, 3);

296 EDGE(1, 2);

297 EDGE(0, 2);

298 break;

299

300 case 4:

301 case INV(4):

302 EDGE(2, 0);

303 EDGE(2, 1);

304 EDGE(2, 3);

305 break;

306

307 case 5:

308 case INV(5):

309 EDGE(0, 1);

310 EDGE(2, 3);

311 EDGE(0, 3);

312

313 EDGE(0, 1);

314 EDGE(1, 2);

315 EDGE(2, 3);

316 break;

317

318 case 6:

319 case INV(6):

320 EDGE(0, 1);

321 EDGE(1, 3);

322 EDGE(2, 3);

323

324 EDGE(0, 1);

325 EDGE(0, 2);

326 EDGE(2, 3);

327 break;

328

329 case 7:

330 case INV(7):

331 EDGE(3, 0);

332 EDGE(3, 2);

333 EDGE(3, 1);

334 break;

335

336 default:

337 break; /* cases 0 and 15 */

338 }

339 }

340

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

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

343 {

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

345 }

346

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

348 {

349 int i;

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

351

352 vec3 p;

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

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

355 }

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

357

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

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

360 }

361

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

363 ms->nverts = 0;

364

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

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

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

368 }

369 }*/

370 }

371

372 #endif /* USE_MTETRA */