GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
5 Algorithm (lower level functions)
6
7 GRASS OpenGL gsurf OGSF Library
8
9 Based on implementation of MarchingCubes 33 Algorithm by
10 Thomas Lewiner, thomas.lewiner@polytechnique.org, Math Dept, PUC-Rio
11
12 (C) 1999-2008 by the GRASS Development Team
13
14 This program is free software under the
15 GNU General Public License (>=v2).
16 Read the file COPYING that comes with GRASS
17 for details.
18
19 \author Tomas Paudits, 2004
20 \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
21 */
22
23#include <float.h>
24
25#include <grass/gis.h>
26#include <grass/ogsf.h>
27
28#include "mc33_table.h"
29
30unsigned char m_case, m_config, m_subconfig;
31
32/*!
33 \brief ADD
34
35 \param face
36 \param v
37
38 \return
39 */
40int mc33_test_face(char face, float *v)
41{
42 float A, B, C, D;
43
44 switch (face) {
45 case -1:
46 case 1:
47 A = v[0];
48 B = v[4];
49 C = v[5];
50 D = v[1];
51 break;
52
53 case -2:
54 case 2:
55 A = v[1];
56 B = v[5];
57 C = v[6];
58 D = v[2];
59 break;
60
61 case -3:
62 case 3:
63 A = v[2];
64 B = v[6];
65 C = v[7];
66 D = v[3];
67 break;
68
69 case -4:
70 case 4:
71 A = v[3];
72 B = v[7];
73 C = v[4];
74 D = v[0];
75 break;
76
77 case -5:
78 case 5:
79 A = v[0];
80 B = v[3];
81 C = v[2];
82 D = v[1];
83 break;
84
85 case -6:
86 case 6:
87 A = v[4];
88 B = v[7];
89 C = v[6];
90 D = v[5];
91 break;
92
93 default:
94 fprintf(stderr, "Invalid face code %d\n", face);
95 A = B = C = D = 0;
96 };
97
98 return face * A * (A * C - B * D) >= 0;
99}
100
101/*!
102 \brief ADD
103
104 \param s
105 \param v
106
107 \return
108 */
109int mc33_test_interior(char s, float *v)
110{
111 float t, At = 0, Bt = 0, Ct = 0, Dt = 0, a, b;
112 char test = 0;
113 char edge = -1;
114
115 switch (m_case) {
116 case 4:
117 case 10:
118 a = (v[4] - v[0]) * (v[6] - v[2]) - (v[7] - v[3]) * (v[5] - v[1]);
119 b = v[2] * (v[4] - v[0]) + v[0] * (v[6] - v[2]) - v[1] * (v[7] - v[3]) -
120 v[3] * (v[5] - v[1]);
121 t = -b / (2 * a);
122
123 if (t < 0 || t > 1)
124 return s > 0;
125
126 At = v[0] + (v[4] - v[0]) * t;
127 Bt = v[3] + (v[7] - v[3]) * t;
128 Ct = v[2] + (v[6] - v[2]) * t;
129 Dt = v[1] + (v[5] - v[1]) * t;
130 break;
131
132 case 6:
133 case 7:
134 case 12:
135 case 13:
136 switch (m_case) {
137 case 6:
138 edge = cell_table[OFFSET_T6_1_1 + m_config].polys[0];
139 break;
140 case 7:
141 edge = cell_table[OFFSET_T7_4_1 + m_config].polys[13];
142 break;
143 case 12:
144 edge = cell_table[OFFSET_T12_2_S1 + m_config].polys[14];
145 break;
146 case 13:
148 .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 break; /* will not reach this as previous switch is exhaustive */
378
379 case 8:
380 return OFFSET_T8 + m_config;
381
382 case 9:
383 return OFFSET_T9 + m_config;
384
385 case 10:
386 if (mc33_test_face(test[OFFSET_TEST10 + m_config][0], v)) {
387 if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v))
388 return OFFSET_T10_1_1_S2 + m_config; /* 10.1.1 */
389 else {
390 return OFFSET_T10_2_S1 + m_config; /* 10.2 */
391 }
392 }
393 else {
394 if (mc33_test_face(test[OFFSET_TEST10 + m_config][1], v)) {
395 return OFFSET_T10_2_S2 + m_config; /* 10.2 */
396 }
397 else {
398 if (mc33_test_interior(test[OFFSET_TEST10 + m_config][2], v))
399 return OFFSET_T10_1_1_S1 + m_config; /* 10.1.1 */
400 else
401 return OFFSET_T10_1_2 + m_config; /* 10.1.2 */
402 }
403 }
404
405 case 11:
406 return OFFSET_T11 + m_config;
407
408 case 12:
409 if (mc33_test_face(test[OFFSET_TEST12 + m_config][0], v)) {
410 if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v))
411 return OFFSET_T12_1_1_S2 + m_config; /* 12.1.1 */
412 else {
413 return OFFSET_T12_2_S1 + m_config; /* 12.2 */
414 }
415 }
416 else {
417 if (mc33_test_face(test[OFFSET_TEST12 + m_config][1], v)) {
418 return OFFSET_T12_2_S2 + m_config; /* 12.2 */
419 }
420 else {
421 if (mc33_test_interior(test[OFFSET_TEST12 + m_config][2], v))
422 return OFFSET_T12_1_1_S1 + m_config; /* 12.1.1 */
423 else
424 return OFFSET_T12_1_2 + m_config; /* 12.1.2 */
425 }
426 }
427
428 case 13:
429 if (mc33_test_face(test[OFFSET_TEST13 + m_config][0], v))
430 m_subconfig += 1;
431 if (mc33_test_face(test[OFFSET_TEST13 + m_config][1], v))
432 m_subconfig += 2;
433 if (mc33_test_face(test[OFFSET_TEST13 + m_config][2], v))
434 m_subconfig += 4;
435 if (mc33_test_face(test[OFFSET_TEST13 + m_config][3], v))
436 m_subconfig += 8;
437 if (mc33_test_face(test[OFFSET_TEST13 + m_config][4], v))
438 m_subconfig += 16;
439 if (mc33_test_face(test[OFFSET_TEST13 + m_config][5], v))
440 m_subconfig += 32;
441
442 switch (subconfig13[m_subconfig]) {
443 case 0: /* 13.1 */
444 return OFFSET_T13_1_S1 + m_config;
445
446 case 1: /* 13.2 */
447 return OFFSET_T13_2_S1 + 0 + m_config * 6;
448 case 2: /* 13.2 */
449 return OFFSET_T13_2_S1 + 1 + m_config * 6;
450 case 3: /* 13.2 */
451 return OFFSET_T13_2_S1 + 2 + m_config * 6;
452 case 4: /* 13.2 */
453 return OFFSET_T13_2_S1 + 3 + m_config * 6;
454 case 5: /* 13.2 */
455 return OFFSET_T13_2_S1 + 4 + m_config * 6;
456 case 6: /* 13.2 */
457 return OFFSET_T13_2_S1 + 5 + m_config * 6;
458
459 case 7: /* 13.3 */
460 return OFFSET_T13_3_S1 + 0 + m_config * 12;
461 case 8: /* 13.3 */
462 return OFFSET_T13_3_S1 + 1 + m_config * 12;
463 case 9: /* 13.3 */
464 return OFFSET_T13_3_S1 + 2 + m_config * 12;
465 case 10: /* 13.3 */
466 return OFFSET_T13_3_S1 + 3 + m_config * 12;
467 case 11: /* 13.3 */
468 return OFFSET_T13_3_S1 + 4 + m_config * 12;
469 case 12: /* 13.3 */
470 return OFFSET_T13_3_S1 + 5 + m_config * 12;
471 case 13: /* 13.3 */
472 return OFFSET_T13_3_S1 + 6 + m_config * 12;
473 case 14: /* 13.3 */
474 return OFFSET_T13_3_S1 + 7 + m_config * 12;
475 case 15: /* 13.3 */
476 return OFFSET_T13_3_S1 + 8 + m_config * 12;
477 case 16: /* 13.3 */
478 return OFFSET_T13_3_S1 + 9 + m_config * 12;
479 case 17: /* 13.3 */
480 return OFFSET_T13_3_S1 + 10 + m_config * 12;
481 case 18: /* 13.3 */
482 return OFFSET_T13_3_S1 + 11 + m_config * 12;
483
484 case 19: /* 13.4 */
485 return OFFSET_T13_4 + 0 + m_config * 4;
486 case 20: /* 13.4 */
487 return OFFSET_T13_4 + 1 + m_config * 4;
488 case 21: /* 13.4 */
489 return OFFSET_T13_4 + 2 + m_config * 4;
490 case 22: /* 13.4 */
491 return OFFSET_T13_4 + 3 + m_config * 4;
492
493 case 23: /* 13.5 */
494 m_subconfig = 0;
495 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
496 return OFFSET_T13_5_1 + 0 + m_config * 4;
497 else
498 return OFFSET_T13_5_2 + 0 + m_config * 4;
499 case 24: /* 13.5 */
500 m_subconfig = 1;
501 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
502 return OFFSET_T13_5_1 + 1 + m_config * 4;
503 else
504 return OFFSET_T13_5_2 + 1 + m_config * 4;
505 case 25: /* 13.5 */
506 m_subconfig = 2;
507 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
508 return OFFSET_T13_5_1 + 2 + m_config * 4;
509 else
510 return OFFSET_T13_5_2 + 2 + m_config * 4;
511 case 26: /* 13.5 */
512 m_subconfig = 3;
513 if (mc33_test_interior(test[OFFSET_TEST13 + m_config][6], v))
514 return OFFSET_T13_5_1 + 3 + m_config * 4;
515 else
516 return OFFSET_T13_5_2 + 3 + m_config * 4;
517
518 case 27: /* 13.3 */
519 return OFFSET_T13_3_S2 + 0 + m_config * 12;
520 case 28: /* 13.3 */
521 return OFFSET_T13_3_S2 + 1 + m_config * 12;
522 case 29: /* 13.3 */
523 return OFFSET_T13_3_S2 + 2 + m_config * 12;
524 case 30: /* 13.3 */
525 return OFFSET_T13_3_S2 + 3 + m_config * 12;
526 case 31: /* 13.3 */
527 return OFFSET_T13_3_S2 + 4 + m_config * 12;
528 case 32: /* 13.3 */
529 return OFFSET_T13_3_S2 + 5 + m_config * 12;
530 case 33: /* 13.3 */
531 return OFFSET_T13_3_S2 + 6 + m_config * 12;
532 case 34: /* 13.3 */
533 return OFFSET_T13_3_S2 + 7 + m_config * 12;
534 case 35: /* 13.3 */
535 return OFFSET_T13_3_S2 + 8 + m_config * 12;
536 case 36: /* 13.3 */
537 return OFFSET_T13_3_S2 + 9 + m_config * 12;
538 case 37: /* 13.3 */
539 return OFFSET_T13_3_S2 + 10 + m_config * 12;
540 case 38: /* 13.3 */
541 return OFFSET_T13_3_S2 + 11 + m_config * 12;
542
543 case 39: /* 13.2 */
544 return OFFSET_T13_2_S2 + 0 + m_config * 6;
545 case 40: /* 13.2 */
546 return OFFSET_T13_2_S2 + 1 + m_config * 6;
547 case 41: /* 13.2 */
548 return OFFSET_T13_2_S2 + 2 + m_config * 6;
549 case 42: /* 13.2 */
550 return OFFSET_T13_2_S2 + 3 + m_config * 6;
551 case 43: /* 13.2 */
552 return OFFSET_T13_2_S2 + 4 + m_config * 6;
553 case 44: /* 13.2 */
554 return OFFSET_T13_2_S2 + 5 + m_config * 6;
555
556 case 45: /* 13.1 */
557 return OFFSET_T13_1_S2 + m_config;
558
559 default:
560 fprintf(stderr, "Marching Cubes: Impossible case 13?\n");
561 }
562 break; /* will not reach this as previous switch is exhaustive */
563
564 case 14:
565 return OFFSET_T14 + m_config;
566 }
567
568 return -1;
569}
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:40
int mc33_process_cube(int c_ndx, float *v)
ADD.
Definition gvl_calc2.c:307
unsigned char m_config
Definition gvl_calc2.c:30
int mc33_test_interior(char s, float *v)
ADD.
Definition gvl_calc2.c:109
unsigned char m_case
Definition gvl_calc2.c:30
unsigned char m_subconfig
Definition gvl_calc2.c:30
#define D
Definition intersect.c:72
OGSF library -.
#define OFFSET_T6_1_1
Definition mc33_table.h:82
#define OFFSET_T12_1_1_S1
Definition mc33_table.h:102
#define OFFSET_T5
Definition mc33_table.h:81
#define OFFSET_T10_1_2
Definition mc33_table.h:98
#define OFFSET_T13_5_2
Definition mc33_table.h:115
#define OFFSET_T13_2_S1
Definition mc33_table.h:109
#define OFFSET_T3_2
Definition mc33_table.h:78
#define OFFSET_T7_3_S1
Definition mc33_table.h:89
#define OFFSET_T10_2_S2
Definition mc33_table.h:100
#define OFFSET_TEST3
#define OFFSET_T10_2_S1
Definition mc33_table.h:99
#define OFFSET_T4_1
Definition mc33_table.h:79
#define OFFSET_T4_2
Definition mc33_table.h:80
#define OFFSET_T13_5_1
Definition mc33_table.h:114
#define OFFSET_T8
Definition mc33_table.h:94
#define OFFSET_T9
Definition mc33_table.h:95
#define OFFSET_T13_4
Definition mc33_table.h:113
#define OFFSET_TEST13
#define OFFSET_T13_1_S1
Definition mc33_table.h:107
#define OFFSET_T13_2_S2
Definition mc33_table.h:110
#define OFFSET_T7_4_1
Definition mc33_table.h:92
#define OFFSET_T10_1_1_S1
Definition mc33_table.h:96
#define OFFSET_TEST4
#define OFFSET_T10_1_1_S2
Definition mc33_table.h:97
#define OFFSET_T1
Definition mc33_table.h:75
#define OFFSET_T12_1_2
Definition mc33_table.h:104
#define OFFSET_T12_1_1_S2
Definition mc33_table.h:103
#define OFFSET_T2
Definition mc33_table.h:76
#define OFFSET_TEST7
#define OFFSET_T13_1_S2
Definition mc33_table.h:108
#define OFFSET_T3_1
Definition mc33_table.h:77
#define OFFSET_T7_3_S3
Definition mc33_table.h:91
#define OFFSET_T7_4_2
Definition mc33_table.h:93
#define OFFSET_T12_2_S2
Definition mc33_table.h:106
#define OFFSET_T7_2_S2
Definition mc33_table.h:87
#define OFFSET_T11
Definition mc33_table.h:101
#define OFFSET_T13_3_S1
Definition mc33_table.h:111
#define OFFSET_TEST6
#define OFFSET_TEST10
#define OFFSET_T12_2_S1
Definition mc33_table.h:105
#define OFFSET_T7_3_S2
Definition mc33_table.h:90
#define OFFSET_TEST12
#define OFFSET_T13_3_S2
Definition mc33_table.h:112
#define OFFSET_T7_2_S1
Definition mc33_table.h:86
#define OFFSET_T7_2_S3
Definition mc33_table.h:88
#define OFFSET_T6_2
Definition mc33_table.h:84
#define OFFSET_T7_1
Definition mc33_table.h:85
#define OFFSET_T6_1_2
Definition mc33_table.h:83
#define OFFSET_T14
Definition mc33_table.h:116