GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-87d2d42d3b
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 
30 unsigned char m_case, m_config, m_subconfig;
31 
32 /*!
33  \brief ADD
34 
35  \param face
36  \param v
37 
38  \return
39  */
40 int 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  */
109 int 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:
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  */
307 int 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
#define D
Definition: gis/intersect.c:72
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
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
Definition: mc33_table.h:16870
#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
Definition: mc33_table.h:16876
#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
Definition: mc33_table.h:16871
#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
Definition: mc33_table.h:16873
#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
Definition: mc33_table.h:16872
#define OFFSET_TEST10
Definition: mc33_table.h:16874
#define OFFSET_T12_2_S1
Definition: mc33_table.h:105
#define OFFSET_T7_3_S2
Definition: mc33_table.h:90
#define OFFSET_TEST12
Definition: mc33_table.h:16875
#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
double b
Definition: r_raster.c:39
double t
Definition: r_raster.c:39
int polys[30]
Definition: viz.h:71