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