GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71746
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
29 unsigned char m_case, m_config, m_subconfig;
30
31 /*!
33
34  \param face
35  \param v
36
37  \return
38  */
39 int 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 /*!
102
103  \param s
104  \param v
105
106  \return
107  */
108 int 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:
137  edge = cell_table[OFFSET_T6_1_1 + m_config].polys[0];
138  break;
139  case 7:
140  edge = cell_table[OFFSET_T7_4_1 + m_config].polys[13];
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 /*!
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
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 }
#define OFFSET_T12_1_1_S1
Definition: mc33_table.h:323
#define OFFSET_T11
Definition: mc33_table.h:322
#define OFFSET_T13_3_S1
Definition: mc33_table.h:332
OGSF library -.
#define OFFSET_T3_1
Definition: mc33_table.h:298
#define OFFSET_T10_1_1_S2
Definition: mc33_table.h:318
#define OFFSET_TEST12
Definition: mc33_table.h:3264
#define OFFSET_T13_5_2
Definition: mc33_table.h:336
#define OFFSET_TEST13
Definition: mc33_table.h:3265
unsigned char m_case
Definition: gvl_calc2.c:29
#define OFFSET_T12_1_1_S2
Definition: mc33_table.h:324
#define OFFSET_T6_1_1
Definition: mc33_table.h:303
#define OFFSET_TEST10
Definition: mc33_table.h:3263
#define OFFSET_T6_2
Definition: mc33_table.h:305
#define OFFSET_T7_2_S3
Definition: mc33_table.h:309
int mc33_test_interior(char s, float *v)
Definition: gvl_calc2.c:108
int polys[30]
Definition: viz.h:78
#define OFFSET_T7_3_S3
Definition: mc33_table.h:312
#define OFFSET_T12_2_S2
Definition: mc33_table.h:327
#define OFFSET_TEST6
Definition: mc33_table.h:3261
#define OFFSET_T8
Definition: mc33_table.h:315
#define OFFSET_T4_2
Definition: mc33_table.h:301
int mc33_process_cube(int c_ndx, float *v)
Definition: gvl_calc2.c:307
int mc33_test_face(char face, float *v)
Definition: gvl_calc2.c:39
#define OFFSET_T13_2_S2
Definition: mc33_table.h:331
#define OFFSET_T10_1_1_S1
Definition: mc33_table.h:317
#define OFFSET_T7_2_S1
Definition: mc33_table.h:307
#define D
Definition: gis/intersect.c:74
double t
Definition: r_raster.c:39
#define OFFSET_T12_2_S1
Definition: mc33_table.h:326
unsigned char m_config
Definition: gvl_calc2.c:29
#define OFFSET_T9
Definition: mc33_table.h:316
double b
Definition: r_raster.c:39
#define OFFSET_TEST3
Definition: mc33_table.h:3259
unsigned char m_subconfig
Definition: gvl_calc2.c:29
#define OFFSET_T7_4_2
Definition: mc33_table.h:314
#define OFFSET_T7_2_S2
Definition: mc33_table.h:308
#define OFFSET_T6_1_2
Definition: mc33_table.h:304
#define OFFSET_T7_3_S2
Definition: mc33_table.h:311
#define OFFSET_T13_1_S2
Definition: mc33_table.h:329
#define OFFSET_TEST7
Definition: mc33_table.h:3262
#define OFFSET_T13_3_S2
Definition: mc33_table.h:333
#define OFFSET_T4_1
Definition: mc33_table.h:300
#define OFFSET_T13_5_1
Definition: mc33_table.h:335
#define OFFSET_T13_2_S1
Definition: mc33_table.h:330
#define OFFSET_T1
Definition: mc33_table.h:296
#define OFFSET_T7_1
Definition: mc33_table.h:306
#define OFFSET_T7_4_1
Definition: mc33_table.h:313
#define OFFSET_T14
Definition: mc33_table.h:337
#define OFFSET_T2
Definition: mc33_table.h:297
#define OFFSET_T13_1_S1
Definition: mc33_table.h:328
#define OFFSET_T10_2_S2
Definition: mc33_table.h:321
#define OFFSET_T10_2_S1
Definition: mc33_table.h:320
#define OFFSET_T7_3_S1
Definition: mc33_table.h:310
#define OFFSET_T13_4
Definition: mc33_table.h:334
#define OFFSET_T10_1_2
Definition: mc33_table.h:319
#define OFFSET_T3_2
Definition: mc33_table.h:299
CELL_ENTRY cell_table[256]
Definition: cell_table.c:3
#define OFFSET_T5
Definition: mc33_table.h:302
#define OFFSET_TEST4
Definition: mc33_table.h:3260
#define OFFSET_T12_1_2
Definition: mc33_table.h:325