GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
gvl_calc2.c
Go to the documentation of this file.
1/*!
2 \file lib/ogsf/gvl_calc2.c
3
4 \brief OGSF library - loading and manipulating volumes, MarchingCubes 33 Algorithm (lower level functions)
5
6 GRASS OpenGL gsurf OGSF Library
7
8 Based on implementation of MarchingCubes 33 Algorithm by
9 Thomas Lewiner, thomas.lewiner@polytechnique.org, Math Dept, PUC-Rio
10
11 (C) 1999-2008 by the GRASS Development Team
12
13 This program is free software under the
14 GNU General Public License (>=v2).
15 Read the file COPYING that comes with GRASS
16 for details.
17
18 \author Tomas Paudits, 2004
19 \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
20 */
21
22#include <float.h>
23
24#include <grass/gis.h>
25#include <grass/ogsf.h>
26
27#include "mc33_table.h"
28
29unsigned char m_case, m_config, m_subconfig;
30
31/*!
32 \brief ADD
33
34 \param face
35 \param v
36
37 \return
38 */
39int mc33_test_face(char face, float *v)
40{
41 float A, B, C, D;
42
43 switch (face) {
44 case -1:
45 case 1:
46 A = v[0];
47 B = v[4];
48 C = v[5];
49 D = v[1];
50 break;
51
52 case -2:
53 case 2:
54 A = v[1];
55 B = v[5];
56 C = v[6];
57 D = v[2];
58 break;
59
60 case -3:
61 case 3:
62 A = v[2];
63 B = v[6];
64 C = v[7];
65 D = v[3];
66 break;
67
68 case -4:
69 case 4:
70 A = v[3];
71 B = v[7];
72 C = v[4];
73 D = v[0];
74 break;
75
76 case -5:
77 case 5:
78 A = v[0];
79 B = v[3];
80 C = v[2];
81 D = v[1];
82 break;
83
84 case -6:
85 case 6:
86 A = v[4];
87 B = v[7];
88 C = v[6];
89 D = v[5];
90 break;
91
92 default:
93 fprintf(stderr, "Invalid face code %d\n", face);
94 A = B = C = D = 0;
95 };
96
97 return face * A * (A * C - B * D) >= 0;
98}
99
100/*!
101 \brief ADD
102
103 \param s
104 \param v
105
106 \return
107 */
108int mc33_test_interior(char s, float *v)
109{
110 float t, At = 0, Bt = 0, Ct = 0, Dt = 0, a, b;
111 char test = 0;
112 char edge = -1;
113
114 switch (m_case) {
115 case 4:
116 case 10:
117 a = (v[4] - v[0]) * (v[6] - v[2]) - (v[7] - v[3]) * (v[5] - v[1]);
118 b = v[2] * (v[4] - v[0]) + v[0] * (v[6] - v[2])
119 - v[1] * (v[7] - v[3]) - v[3] * (v[5] - v[1]);
120 t = -b / (2 * a);
121
122 if (t < 0 || t > 1)
123 return s > 0;
124
125 At = v[0] + (v[4] - v[0]) * t;
126 Bt = v[3] + (v[7] - v[3]) * t;
127 Ct = v[2] + (v[6] - v[2]) * t;
128 Dt = v[1] + (v[5] - v[1]) * t;
129 break;
130
131 case 6:
132 case 7:
133 case 12:
134 case 13:
135 switch (m_case) {
136 case 6:
138 break;
139 case 7:
141 break;
142 case 12:
144 break;
145 case 13:
146 edge =
148 m_config * 4].polys[2];
149 break;
150 }
151
152 switch (edge) {
153 case 0:
154 t = v[0] / (v[0] - v[1]);
155 At = 0;
156 Bt = v[3] + (v[2] - v[3]) * t;
157 Ct = v[7] + (v[6] - v[7]) * t;
158 Dt = v[4] + (v[5] - v[4]) * t;
159 break;
160 case 1:
161 t = v[1] / (v[1] - v[2]);
162 At = 0;
163 Bt = v[0] + (v[3] - v[0]) * t;
164 Ct = v[4] + (v[7] - v[4]) * t;
165 Dt = v[5] + (v[6] - v[5]) * t;
166 break;
167 case 2:
168 t = v[2] / (v[2] - v[3]);
169 At = 0;
170 Bt = v[1] + (v[0] - v[1]) * t;
171 Ct = v[5] + (v[4] - v[5]) * t;
172 Dt = v[6] + (v[7] - v[6]) * t;
173 break;
174 case 3:
175 t = v[3] / (v[3] - v[0]);
176 At = 0;
177 Bt = v[2] + (v[1] - v[2]) * t;
178 Ct = v[6] + (v[5] - v[6]) * t;
179 Dt = v[7] + (v[4] - v[7]) * t;
180 break;
181 case 4:
182 t = v[4] / (v[4] - v[5]);
183 At = 0;
184 Bt = v[7] + (v[6] - v[7]) * t;
185 Ct = v[3] + (v[2] - v[3]) * t;
186 Dt = v[0] + (v[1] - v[0]) * t;
187 break;
188 case 5:
189 t = v[5] / (v[5] - v[6]);
190 At = 0;
191 Bt = v[4] + (v[7] - v[4]) * t;
192 Ct = v[0] + (v[3] - v[0]) * t;
193 Dt = v[1] + (v[2] - v[1]) * t;
194 break;
195 case 6:
196 t = v[6] / (v[6] - v[7]);
197 At = 0;
198 Bt = v[5] + (v[4] - v[5]) * t;
199 Ct = v[1] + (v[0] - v[1]) * t;
200 Dt = v[2] + (v[3] - v[2]) * t;
201 break;
202 case 7:
203 t = v[7] / (v[7] - v[4]);
204 At = 0;
205 Bt = v[6] + (v[5] - v[6]) * t;
206 Ct = v[2] + (v[1] - v[2]) * t;
207 Dt = v[3] + (v[0] - v[3]) * t;
208 break;
209 case 8:
210 t = v[0] / (v[0] - v[4]);
211 At = 0;
212 Bt = v[3] + (v[7] - v[3]) * t;
213 Ct = v[2] + (v[6] - v[2]) * t;
214 Dt = v[1] + (v[5] - v[1]) * t;
215 break;
216 case 9:
217 t = v[1] / (v[1] - v[5]);
218 At = 0;
219 Bt = v[0] + (v[4] - v[0]) * t;
220 Ct = v[3] + (v[7] - v[3]) * t;
221 Dt = v[2] + (v[6] - v[2]) * t;
222 break;
223 case 10:
224 t = v[2] / (v[2] - v[6]);
225 At = 0;
226 Bt = v[1] + (v[5] - v[1]) * t;
227 Ct = v[0] + (v[4] - v[0]) * t;
228 Dt = v[3] + (v[7] - v[3]) * t;
229 break;
230 case 11:
231 t = v[3] / (v[3] - v[7]);
232 At = 0;
233 Bt = v[2] + (v[6] - v[2]) * t;
234 Ct = v[1] + (v[5] - v[1]) * t;
235 Dt = v[0] + (v[4] - v[0]) * t;
236 break;
237 default:
238 fprintf(stderr, "Invalid edge %d\n", edge);
239 break;
240 }
241 break;
242
243 default:
244 fprintf(stderr, "Invalid ambiguous case %d\n", m_case);
245 break;
246 }
247
248 if (At >= 0)
249 test++;
250 if (Bt >= 0)
251 test += 2;
252 if (Ct >= 0)
253 test += 4;
254 if (Dt >= 0)
255 test += 8;
256
257 switch (test) {
258 case 0:
259 return s > 0;
260 case 1:
261 return s > 0;
262 case 2:
263 return s > 0;
264 case 3:
265 return s > 0;
266 case 4:
267 return s > 0;
268 case 5:
269 if (At * Ct < Bt * Dt)
270 return s > 0;
271 break;
272 case 6:
273 return s > 0;
274 case 7:
275 return s < 0;
276 case 8:
277 return s > 0;
278 case 9:
279 return s > 0;
280 case 10:
281 if (At * Ct >= Bt * Dt)
282 return s > 0;
283 break;
284 case 11:
285 return s < 0;
286 case 12:
287 return s > 0;
288 case 13:
289 return s < 0;
290 case 14:
291 return s < 0;
292 case 15:
293 return s < 0;
294 }
295
296 return s < 0;
297}
298
299/*!
300 \brief ADD
301
302 \param c_ndx
303 \param v
304
305 \return
306 */
307int mc33_process_cube(int c_ndx, float *v)
308{
309 m_case = cases[c_ndx][0];
310 m_config = cases[c_ndx][1];
311 m_subconfig = 0;
312
313 switch (m_case) {
314 case 0:
315 return -1;
316
317 case 1:
318 return OFFSET_T1 + m_config;
319
320 case 2:
321 return OFFSET_T2 + m_config;
322
323 case 3:
324 if (mc33_test_face(test[OFFSET_TEST3 + m_config][0], v))
325 return OFFSET_T3_2 + m_config; /* 3.2 */
326 else
327 return OFFSET_T3_1 + m_config; /* 3.1 */
328
329 case 4:
330 if (mc33_test_interior(test[OFFSET_TEST4 + m_config][0], v))
331 return OFFSET_T4_1 + m_config; /* 4.1 */
332 else
333 return OFFSET_T4_2 + m_config; /* 4.2 */
334
335 case 5:
336 return OFFSET_T5 + m_config;
337
338 case 6:
339 if (mc33_test_face(test[OFFSET_TEST6 + m_config][0], v))
340 return OFFSET_T6_2 + m_config; /* 6.2 */
341 else {
342 if (mc33_test_interior(test[OFFSET_TEST6 + m_config][1], v))
343 return OFFSET_T6_1_1 + m_config; /* 6.1.1 */
344 else
345 return OFFSET_T6_1_2 + m_config; /* 6.1.2 */
346 }
347
348 case 7:
349 if (mc33_test_face(test[OFFSET_TEST7 + m_config][0], v))
350 m_subconfig += 1;
351 if (mc33_test_face(test[OFFSET_TEST7 + m_config][1], v))
352 m_subconfig += 2;
353 if (mc33_test_face(test[OFFSET_TEST7 + m_config][2], v))
354 m_subconfig += 4;
355
356 switch (subconfig7[m_subconfig]) {
357 case 0:
358 if (mc33_test_interior(test[OFFSET_TEST7 + m_config][3], v))
359 return OFFSET_T7_4_2 + m_config; /* 7.4.2 */
360 else
361 return OFFSET_T7_4_1 + m_config; /* 7.4.1 */
362 case 1:
363 return OFFSET_T7_3_S1 + m_config; /* 7.3 */
364 case 2:
365 return OFFSET_T7_3_S2 + m_config; /* 7.3 */
366 case 3:
367 return OFFSET_T7_3_S3 + m_config; /* 7.3 */
368 case 4:
369 return OFFSET_T7_2_S1 + m_config; /* 7.2 */
370 case 5:
371 return OFFSET_T7_2_S2 + m_config; /* 7.2 */
372 case 6:
373 return OFFSET_T7_2_S3 + m_config; /* 7.2 */
374 case 7:
375 return OFFSET_T7_1 + m_config; /* 7.1 */
376 };
377
378 case 8:
379 return OFFSET_T8 + m_config;
380
381 case 9:
382 return OFFSET_T9 + m_config;
383
384 case 10:
385 if (mc33_test_face(test[OFFSET_TEST10 + m_config][0], v)) {
386 if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v))
387 return OFFSET_T10_1_1_S2 + m_config; /* 10.1.1 */
388 else {
389 return OFFSET_T10_2_S1 + m_config; /* 10.2 */
390 }
391 }
392 else {
393 if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) {
394 return OFFSET_T10_2_S2 + m_config; /* 10.2 */
395 }
396 else {
397 if (mc33_test_interior(test[OFFSET_TEST10 + m_config][2], v))
398 return OFFSET_T10_1_1_S1 + m_config; /* 10.1.1 */
399 else
400 return OFFSET_T10_1_2 + m_config; /* 10.1.2 */
401 }
402 }
403
404 case 11:
405 return OFFSET_T11 + m_config;
406
407 case 12:
408 if (mc33_test_face(test[OFFSET_TEST12 + m_config][0], v)) {
409 if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v))
410 return OFFSET_T12_1_1_S2 + m_config; /* 12.1.1 */
411 else {
412 return OFFSET_T12_2_S1 + m_config; /* 12.2 */
413 }
414 }
415 else {
416 if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) {
417 return OFFSET_T12_2_S2 + m_config; /* 12.2 */
418 }
419 else {
420 if (mc33_test_interior(test[OFFSET_TEST12 + m_config][2], v))
421 return OFFSET_T12_1_1_S1 + m_config; /* 12.1.1 */
422 else
423 return OFFSET_T12_1_2 + m_config; /* 12.1.2 */
424 }
425 }
426
427 case 13:
428 if (mc33_test_face(test[OFFSET_TEST13 + m_config][0], v))
429 m_subconfig += 1;
430 if (mc33_test_face(test[OFFSET_TEST13 + m_config][1], v))
431 m_subconfig += 2;
432 if (mc33_test_face(test[OFFSET_TEST13 + m_config][2], v))
433 m_subconfig += 4;
434 if (mc33_test_face(test[OFFSET_TEST13 + m_config][3], v))
435 m_subconfig += 8;
436 if (mc33_test_face(test[OFFSET_TEST13 + m_config][4], v))
437 m_subconfig += 16;
438 if (mc33_test_face(test[OFFSET_TEST13 + m_config][5], v))
439 m_subconfig += 32;
440
441 switch (subconfig13[m_subconfig]) {
442 case 0: /* 13.1 */
443 return OFFSET_T13_1_S1 + m_config;
444
445 case 1: /* 13.2 */
446 return OFFSET_T13_2_S1 + 0 + m_config * 6;
447 case 2: /* 13.2 */
448 return OFFSET_T13_2_S1 + 1 + m_config * 6;
449 case 3: /* 13.2 */
450 return OFFSET_T13_2_S1 + 2 + m_config * 6;
451 case 4: /* 13.2 */
452 return OFFSET_T13_2_S1 + 3 + m_config * 6;
453 case 5: /* 13.2 */
454 return OFFSET_T13_2_S1 + 4 + m_config * 6;
455 case 6: /* 13.2 */
456 return OFFSET_T13_2_S1 + 5 + m_config * 6;
457
458 case 7: /* 13.3 */
459 return OFFSET_T13_3_S1 + 0 + m_config * 12;
460 case 8: /* 13.3 */
461 return OFFSET_T13_3_S1 + 1 + m_config * 12;
462 case 9: /* 13.3 */
463 return OFFSET_T13_3_S1 + 2 + m_config * 12;
464 case 10: /* 13.3 */
465 return OFFSET_T13_3_S1 + 3 + m_config * 12;
466 case 11: /* 13.3 */
467 return OFFSET_T13_3_S1 + 4 + m_config * 12;
468 case 12: /* 13.3 */
469 return OFFSET_T13_3_S1 + 5 + m_config * 12;
470 case 13: /* 13.3 */
471 return OFFSET_T13_3_S1 + 6 + m_config * 12;
472 case 14: /* 13.3 */
473 return OFFSET_T13_3_S1 + 7 + m_config * 12;
474 case 15: /* 13.3 */
475 return OFFSET_T13_3_S1 + 8 + m_config * 12;
476 case 16: /* 13.3 */
477 return OFFSET_T13_3_S1 + 9 + m_config * 12;
478 case 17: /* 13.3 */
479 return OFFSET_T13_3_S1 + 10 + m_config * 12;
480 case 18: /* 13.3 */
481 return OFFSET_T13_3_S1 + 11 + m_config * 12;
482
483 case 19: /* 13.4 */
484 return OFFSET_T13_4 + 0 + m_config * 4;
485 case 20: /* 13.4 */
486 return OFFSET_T13_4 + 1 + m_config * 4;
487 case 21: /* 13.4 */
488 return OFFSET_T13_4 + 2 + m_config * 4;
489 case 22: /* 13.4 */
490 return OFFSET_T13_4 + 3 + m_config * 4;
491
492 case 23: /* 13.5 */
493 m_subconfig = 0;
494 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
495 return OFFSET_T13_5_1 + 0 + m_config * 4;
496 else
497 return OFFSET_T13_5_2 + 0 + m_config * 4;
498 case 24: /* 13.5 */
499 m_subconfig = 1;
500 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
501 return OFFSET_T13_5_1 + 1 + m_config * 4;
502 else
503 return OFFSET_T13_5_2 + 1 + m_config * 4;
504 case 25: /* 13.5 */
505 m_subconfig = 2;
506 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
507 return OFFSET_T13_5_1 + 2 + m_config * 4;
508 else
509 return OFFSET_T13_5_2 + 2 + m_config * 4;
510 case 26: /* 13.5 */
511 m_subconfig = 3;
512 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
513 return OFFSET_T13_5_1 + 3 + m_config * 4;
514 else
515 return OFFSET_T13_5_2 + 3 + m_config * 4;
516
517 case 27: /* 13.3 */
518 return OFFSET_T13_3_S2 + 0 + m_config * 12;
519 case 28: /* 13.3 */
520 return OFFSET_T13_3_S2 + 1 + m_config * 12;
521 case 29: /* 13.3 */
522 return OFFSET_T13_3_S2 + 2 + m_config * 12;
523 case 30: /* 13.3 */
524 return OFFSET_T13_3_S2 + 3 + m_config * 12;
525 case 31: /* 13.3 */
526 return OFFSET_T13_3_S2 + 4 + m_config * 12;
527 case 32: /* 13.3 */
528 return OFFSET_T13_3_S2 + 5 + m_config * 12;
529 case 33: /* 13.3 */
530 return OFFSET_T13_3_S2 + 6 + m_config * 12;
531 case 34: /* 13.3 */
532 return OFFSET_T13_3_S2 + 7 + m_config * 12;
533 case 35: /* 13.3 */
534 return OFFSET_T13_3_S2 + 8 + m_config * 12;
535 case 36: /* 13.3 */
536 return OFFSET_T13_3_S2 + 9 + m_config * 12;
537 case 37: /* 13.3 */
538 return OFFSET_T13_3_S2 + 10 + m_config * 12;
539 case 38: /* 13.3 */
540 return OFFSET_T13_3_S2 + 11 + m_config * 12;
541
542 case 39: /* 13.2 */
543 return OFFSET_T13_2_S2 + 0 + m_config * 6;
544 case 40: /* 13.2 */
545 return OFFSET_T13_2_S2 + 1 + m_config * 6;
546 case 41: /* 13.2 */
547 return OFFSET_T13_2_S2 + 2 + m_config * 6;
548 case 42: /* 13.2 */
549 return OFFSET_T13_2_S2 + 3 + m_config * 6;
550 case 43: /* 13.2 */
551 return OFFSET_T13_2_S2 + 4 + m_config * 6;
552 case 44: /* 13.2 */
553 return OFFSET_T13_2_S2 + 5 + m_config * 6;
554
555 case 45: /* 13.1 */
556 return OFFSET_T13_1_S2 + m_config;
557
558 default:
559 fprintf(stderr, "Marching Cubes: Impossible case 13?\n");
560 }
561
562 case 14:
563 return OFFSET_T14 + m_config;
564 }
565
566 return -1;
567}
CELL_ENTRY cell_table[256]
Definition: cell_table.c:3
double b
double t
int mc33_test_face(char face, float *v)
ADD.
Definition: gvl_calc2.c:39
int mc33_process_cube(int c_ndx, float *v)
ADD.
Definition: gvl_calc2.c:307
unsigned char m_config
Definition: gvl_calc2.c:29
int mc33_test_interior(char s, float *v)
ADD.
Definition: gvl_calc2.c:108
unsigned char m_case
Definition: gvl_calc2.c:29
unsigned char m_subconfig
Definition: gvl_calc2.c:29
#define D
Definition: intersect.c:74
OGSF library -.
#define OFFSET_T6_1_1
Definition: mc33_table.h:303
#define OFFSET_T12_1_1_S1
Definition: mc33_table.h:323
#define OFFSET_T5
Definition: mc33_table.h:302
#define OFFSET_T10_1_2
Definition: mc33_table.h:319
#define OFFSET_T13_5_2
Definition: mc33_table.h:336
#define OFFSET_T13_2_S1
Definition: mc33_table.h:330
#define OFFSET_T3_2
Definition: mc33_table.h:299
#define OFFSET_T7_3_S1
Definition: mc33_table.h:310
#define OFFSET_T10_2_S2
Definition: mc33_table.h:321
#define OFFSET_TEST3
Definition: mc33_table.h:3259
#define OFFSET_T10_2_S1
Definition: mc33_table.h:320
#define OFFSET_T4_1
Definition: mc33_table.h:300
#define OFFSET_T4_2
Definition: mc33_table.h:301
#define OFFSET_T13_5_1
Definition: mc33_table.h:335
#define OFFSET_T8
Definition: mc33_table.h:315
#define OFFSET_T9
Definition: mc33_table.h:316
#define OFFSET_T13_4
Definition: mc33_table.h:334
#define OFFSET_TEST13
Definition: mc33_table.h:3265
#define OFFSET_T13_1_S1
Definition: mc33_table.h:328
#define OFFSET_T13_2_S2
Definition: mc33_table.h:331
#define OFFSET_T7_4_1
Definition: mc33_table.h:313
#define OFFSET_T10_1_1_S1
Definition: mc33_table.h:317
#define OFFSET_TEST4
Definition: mc33_table.h:3260
#define OFFSET_T10_1_1_S2
Definition: mc33_table.h:318
#define OFFSET_T1
Definition: mc33_table.h:296
#define OFFSET_T12_1_2
Definition: mc33_table.h:325
#define OFFSET_T12_1_1_S2
Definition: mc33_table.h:324
#define OFFSET_T2
Definition: mc33_table.h:297
#define OFFSET_TEST7
Definition: mc33_table.h:3262
#define OFFSET_T13_1_S2
Definition: mc33_table.h:329
#define OFFSET_T3_1
Definition: mc33_table.h:298
#define OFFSET_T7_3_S3
Definition: mc33_table.h:312
#define OFFSET_T7_4_2
Definition: mc33_table.h:314
#define OFFSET_T12_2_S2
Definition: mc33_table.h:327
#define OFFSET_T7_2_S2
Definition: mc33_table.h:308
#define OFFSET_T11
Definition: mc33_table.h:322
#define OFFSET_T13_3_S1
Definition: mc33_table.h:332
#define OFFSET_TEST6
Definition: mc33_table.h:3261
#define OFFSET_TEST10
Definition: mc33_table.h:3263
#define OFFSET_T12_2_S1
Definition: mc33_table.h:326
#define OFFSET_T7_3_S2
Definition: mc33_table.h:311
#define OFFSET_TEST12
Definition: mc33_table.h:3264
#define OFFSET_T13_3_S2
Definition: mc33_table.h:333
#define OFFSET_T7_2_S1
Definition: mc33_table.h:307
#define OFFSET_T7_2_S3
Definition: mc33_table.h:309
#define OFFSET_T6_2
Definition: mc33_table.h:305
#define OFFSET_T7_1
Definition: mc33_table.h:306
#define OFFSET_T6_1_2
Definition: mc33_table.h:304
#define OFFSET_T14
Definition: mc33_table.h:337
int polys[30]
Definition: viz.h:78