PeriDyno 1.0.0
Loading...
Searching...
No Matches
MarchingCubesHelper.cu
Go to the documentation of this file.
1#include "MarchingCubesHelper.h"
2
3#include "Topology/EdgeSet.h"
4
5#include <thrust/sort.h>
6
7namespace dyno
8{
9 __device__
10 uint edgeTable[256] =
11 {
12 0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
13 0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
14 0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
15 0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
16 0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
17 0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
18 0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
19 0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
20 0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
21 0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
22 0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
23 0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
24 0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
25 0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
26 0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
27 0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
28 0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
29 0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
30 0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
31 0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
32 0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
33 0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
34 0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
35 0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
36 0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
37 0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
38 0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
39 0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
40 0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
41 0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
42 0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
43 0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0
44 };
45
46 __device__
47 int triTable[256][16] =
48 {
49 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
50 {0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
51 {0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
52 {1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
53 {1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
54 {0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
55 {9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
56 {2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
57 {3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
58 {0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
59 {1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
60 {1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
61 {3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
62 {0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
63 {3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
64 {9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
65 {4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
66 {4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
67 {0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
68 {4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
69 {1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
70 {3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
71 {9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
72 {2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
73 {8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
74 {11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
75 {9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
76 {4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
77 {3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
78 {1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
79 {4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
80 {4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
81 {9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
82 {9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
83 {0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
84 {8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
85 {1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
86 {3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
87 {5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
88 {2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
89 {9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
90 {0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
91 {0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
92 {2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
93 {10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
94 {4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
95 {5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
96 {5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
97 {9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
98 {9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
99 {0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
100 {1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
101 {9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
102 {10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
103 {8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
104 {2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
105 {7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
106 {9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
107 {2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
108 {11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
109 {9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
110 {5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
111 {11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
112 {11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
113 {10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
114 {0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
115 {9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
116 {1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
117 {1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
118 {1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
119 {9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
120 {5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
121 {2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
122 {11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
123 {0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
124 {5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
125 {6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
126 {0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
127 {3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
128 {6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
129 {5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
130 {4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
131 {1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
132 {10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
133 {6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
134 {1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
135 {8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
136 {7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
137 {3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
138 {5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
139 {0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
140 {9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
141 {8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
142 {5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
143 {0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
144 {6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
145 {10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
146 {4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
147 {10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
148 {8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
149 {1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
150 {3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
151 {0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
152 {8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
153 {10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
154 {0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
155 {3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
156 {6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
157 {9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
158 {8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
159 {3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
160 {6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
161 {7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
162 {0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
163 {10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
164 {10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
165 {1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
166 {2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
167 {7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
168 {7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
169 {2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
170 {2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
171 {1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
172 {11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
173 {8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
174 {0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
175 {7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
176 {7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
177 {7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
178 {3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
179 {0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
180 {8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
181 {10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
182 {1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
183 {2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
184 {6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
185 {7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
186 {7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
187 {2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
188 {1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
189 {10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
190 {10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
191 {0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
192 {7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
193 {6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
194 {3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
195 {8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
196 {9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
197 {6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
198 {1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
199 {4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
200 {10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
201 {8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
202 {0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
203 {1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
204 {1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
205 {8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
206 {10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
207 {4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
208 {10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
209 {4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
210 {0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
211 {5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
212 {11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
213 {9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
214 {6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
215 {7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
216 {3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
217 {7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
218 {9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
219 {3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
220 {6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
221 {9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
222 {1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
223 {4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
224 {7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
225 {6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
226 {3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
227 {0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
228 {6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
229 {1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
230 {0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
231 {11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
232 {6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
233 {5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
234 {9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
235 {1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
236 {1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
237 {1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
238 {10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
239 {0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
240 {10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
241 {11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
242 {11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
243 {5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
244 {10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
245 {11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
246 {0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
247 {9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
248 {7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
249 {2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
250 {8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
251 {9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
252 {9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
253 {1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
254 {0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
255 {9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
256 {9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
257 {5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
258 {5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
259 {0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
260 {10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
261 {2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
262 {0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
263 {0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
264 {9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
265 {2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
266 {5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
267 {3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
268 {5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
269 {8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
270 {0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
271 {8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
272 {9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
273 {4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
274 {0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
275 {1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
276 {3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
277 {4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
278 {9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
279 {11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
280 {11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
281 {2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
282 {9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
283 {3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
284 {1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
285 {4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
286 {4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
287 {4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
288 {4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
289 {9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
290 {3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
291 {0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
292 {3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
293 {1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
294 {3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
295 {0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
296 {3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
297 {2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
298 {9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
299 {2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
300 {1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
301 {1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
302 {0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
303 {0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
304 {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
305 };
306
307 __device__
308 uint numVertsTable[256] =
309 {
310 0,
311 3,
312 3,
313 6,
314 3,
315 6,
316 6,
317 9,
318 3,
319 6,
320 6,
321 9,
322 6,
323 9,
324 9,
325 6,
326 3,
327 6,
328 6,
329 9,
330 6,
331 9,
332 9,
333 12,
334 6,
335 9,
336 9,
337 12,
338 9,
339 12,
340 12,
341 9,
342 3,
343 6,
344 6,
345 9,
346 6,
347 9,
348 9,
349 12,
350 6,
351 9,
352 9,
353 12,
354 9,
355 12,
356 12,
357 9,
358 6,
359 9,
360 9,
361 6,
362 9,
363 12,
364 12,
365 9,
366 9,
367 12,
368 12,
369 9,
370 12,
371 15,
372 15,
373 6,
374 3,
375 6,
376 6,
377 9,
378 6,
379 9,
380 9,
381 12,
382 6,
383 9,
384 9,
385 12,
386 9,
387 12,
388 12,
389 9,
390 6,
391 9,
392 9,
393 12,
394 9,
395 12,
396 12,
397 15,
398 9,
399 12,
400 12,
401 15,
402 12,
403 15,
404 15,
405 12,
406 6,
407 9,
408 9,
409 12,
410 9,
411 12,
412 6,
413 9,
414 9,
415 12,
416 12,
417 15,
418 12,
419 15,
420 9,
421 6,
422 9,
423 12,
424 12,
425 9,
426 12,
427 15,
428 9,
429 6,
430 12,
431 15,
432 15,
433 12,
434 15,
435 6,
436 12,
437 3,
438 3,
439 6,
440 6,
441 9,
442 6,
443 9,
444 9,
445 12,
446 6,
447 9,
448 9,
449 12,
450 9,
451 12,
452 12,
453 9,
454 6,
455 9,
456 9,
457 12,
458 9,
459 12,
460 12,
461 15,
462 9,
463 6,
464 12,
465 9,
466 12,
467 9,
468 15,
469 6,
470 6,
471 9,
472 9,
473 12,
474 9,
475 12,
476 12,
477 15,
478 9,
479 12,
480 12,
481 15,
482 12,
483 15,
484 15,
485 12,
486 9,
487 12,
488 12,
489 9,
490 12,
491 15,
492 15,
493 12,
494 12,
495 9,
496 15,
497 6,
498 15,
499 12,
500 6,
501 3,
502 6,
503 9,
504 9,
505 12,
506 9,
507 12,
508 12,
509 15,
510 9,
511 12,
512 12,
513 15,
514 6,
515 9,
516 9,
517 6,
518 9,
519 12,
520 12,
521 15,
522 12,
523 15,
524 15,
525 6,
526 12,
527 9,
528 15,
529 12,
530 9,
531 6,
532 12,
533 3,
534 9,
535 12,
536 12,
537 15,
538 12,
539 15,
540 9,
541 12,
542 12,
543 15,
544 15,
545 6,
546 9,
547 12,
548 6,
549 3,
550 6,
551 9,
552 9,
553 6,
554 9,
555 12,
556 6,
557 3,
558 9,
559 6,
560 12,
561 3,
562 6,
563 3,
564 3,
565 0,
566 };
567
568 template<typename Real, typename Coord, typename TDataType>
569 __global__ void RconstructSDF(
570 DArray3D<Real> distances,
571 Coord origin,
572 Real h,
573 DistanceField3D<TDataType> sdf)
574 {
575 uint i = threadIdx.x + blockDim.x * blockIdx.x;
576 uint j = threadIdx.y + blockDim.y * blockIdx.y;
577 uint k = threadIdx.z + blockDim.z * blockIdx.z;
578
579 uint nx = distances.nx();
580 uint ny = distances.ny();
581 uint nz = distances.nz();
582
583 if (i >= nx || j >= ny || k >= nz) return;
584
585 Coord p = origin + Coord(i * h, j * h, k * h);
586
587 Real d;
588 Coord normal;
589 sdf.getDistance(p, d, normal);
590 distances(i, j, k) = d;
591 }
592
593 template<typename TDataType>
594 void MarchingCubesHelper<TDataType>::reconstructSDF(
595 DArray3D<Real>& distances,
596 Coord origin,
597 Real h,
598 DistanceField3D<TDataType>& sdf)
599 {
600 cuExecute3D(make_uint3(distances.nx(), distances.ny(), distances.nz()),
601 RconstructSDF,
602 distances,
603 origin,
604 h,
605 sdf);
606 }
607
608 template<typename Real>
609 __global__ void CountVerticeNumber(
610 DArray<int> num,
611 DArray3D<Real> distances,
612 Real isoValue)
613 {
614 uint i = threadIdx.x + blockDim.x * blockIdx.x;
615 uint j = threadIdx.y + blockDim.y * blockIdx.y;
616 uint k = threadIdx.z + blockDim.z * blockIdx.z;
617
618 uint nx = distances.nx();
619 uint ny = distances.ny();
620 uint nz = distances.nz();
621
622 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
623
624 uint cubeindex;
625 cubeindex = uint(distances(i, j, k) < isoValue);
626 cubeindex += uint(distances(i + 1, j, k) < isoValue) * 2;
627 cubeindex += uint(distances(i + 1, j + 1, k) < isoValue) * 4;
628 cubeindex += uint(distances(i, j + 1, k) < isoValue) * 8;
629 cubeindex += uint(distances(i, j, k + 1) < isoValue) * 16;
630 cubeindex += uint(distances(i + 1, j , k + 1) < isoValue) * 32;
631 cubeindex += uint(distances(i + 1, j + 1, k + 1) < isoValue) * 64;
632 cubeindex += uint(distances(i, j + 1, k + 1) < isoValue) * 128;
633
634 num[i + j * (nx - 1) + k * (nx - 1) * (ny - 1)] = numVertsTable[cubeindex];
635 }
636
637 template<typename TDataType>
638 void MarchingCubesHelper<TDataType>::countVerticeNumber(
639 DArray<int>& num,
640 DArray3D<Real>& distances,
641 Real isoValue)
642 {
643 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
644 CountVerticeNumber,
645 num,
646 distances,
647 isoValue);
648 }
649
650 // compute interpolated vertex along an edge
651 template<typename Real, typename Coord>
652 __device__ Coord vertexInterp(Real isolevel, Coord p0, Coord p1, Real f0, Real f1)
653 {
654 if (abs(f1 - f0) < EPSILON)
655 return 0.5 * (p0 + p1);
656
657 Real t = (isolevel - f0) / (f1 - f0);
658
659 t = clamp(t, Real(0), Real(1));
660
661 return (1 - t) * p0 + t * p1;
662 }
663
664 // compute interpolated vertex along an edge
665 template<typename Real, typename Coord>
666 __device__ Coord vertexInterp(Real& value, Coord p0, Coord p1, Real f0, Real f1, Real d0, Real d1)
667 {
668 if (abs(d0 - d1) < EPSILON)
669 {
670 value = 0.5 * (f0 + f1);
671 return 0.5 * (p0 + p1);
672 }
673
674
675 Real t = d0 / (d0 - d1);
676
677 t = clamp(t, Real(0), Real(1));
678
679 value = (1 - t) * f0 + t * f1;
680
681 return (1 - t) * p0 + t * p1;
682 }
683
684 template<typename Coord, typename Triangle>
685 __global__ void ConstructTriangles(
686 DArray<Coord> vertices,
687 DArray<EKey> edgeIds,
688 DArray<Triangle> triangles,
689 DArray<int> vertNum,
690 DArray3D<Real> distances,
691 Coord origin,
692 Real isoValue,
693 Real h)
694 {
695 uint i = threadIdx.x + blockDim.x * blockIdx.x;
696 uint j = threadIdx.y + blockDim.y * blockIdx.y;
697 uint k = threadIdx.z + blockDim.z * blockIdx.z;
698
699 uint nx = distances.nx();
700 uint ny = distances.ny();
701 uint nz = distances.nz();
702
703 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
704
705 Coord v[8];
706 Coord p = origin + h * Coord(i, j, k);
707 v[0] = p;
708 v[1] = p + Coord(h, 0, 0);
709 v[2] = p + Coord(h, h, 0);
710 v[3] = p + Coord(0, h, 0);
711 v[4] = p + Coord(0, 0, h);
712 v[5] = p + Coord(h, 0, h);
713 v[6] = p + Coord(h, h, h);
714 v[7] = p + Coord(0, h, h);
715
716 Real field[8];
717 field[0] = distances(i, j, k);
718 field[1] = distances(i + 1, j, k);
719 field[2] = distances(i + 1, j + 1, k);
720 field[3] = distances(i, j + 1, k);
721 field[4] = distances(i, j, k + 1);
722 field[5] = distances(i + 1, j, k + 1);
723 field[6] = distances(i + 1, j + 1, k + 1);
724 field[7] = distances(i, j + 1, k + 1);
725
726 uint vIds[8];
727 vIds[0] = distances.index(i, j, k);
728 vIds[1] = distances.index(i + 1, j, k);
729 vIds[2] = distances.index(i + 1, j + 1, k);
730 vIds[3] = distances.index(i, j + 1, k);
731 vIds[4] = distances.index(i, j, k + 1);
732 vIds[5] = distances.index(i + 1, j, k + 1);
733 vIds[6] = distances.index(i + 1, j + 1, k + 1);
734 vIds[7] = distances.index(i, j + 1, k + 1);
735
736 Coord vertlist[12];
737 vertlist[0] = vertexInterp(isoValue, v[0], v[1], field[0], field[1]);
738 vertlist[1] = vertexInterp(isoValue, v[1], v[2], field[1], field[2]);
739 vertlist[2] = vertexInterp(isoValue, v[2], v[3], field[2], field[3]);
740 vertlist[3] = vertexInterp(isoValue, v[3], v[0], field[3], field[0]);
741
742 vertlist[4] = vertexInterp(isoValue, v[4], v[5], field[4], field[5]);
743 vertlist[5] = vertexInterp(isoValue, v[5], v[6], field[5], field[6]);
744 vertlist[6] = vertexInterp(isoValue, v[6], v[7], field[6], field[7]);
745 vertlist[7] = vertexInterp(isoValue, v[7], v[4], field[7], field[4]);
746
747 vertlist[8] = vertexInterp(isoValue, v[0], v[4], field[0], field[4]);
748 vertlist[9] = vertexInterp(isoValue, v[1], v[5], field[1], field[5]);
749 vertlist[10] = vertexInterp(isoValue, v[2], v[6], field[2], field[6]);
750 vertlist[11] = vertexInterp(isoValue, v[3], v[7], field[3], field[7]);
751
752 EKey edgelist[12];
753 edgelist[0] = EKey(vIds[0], vIds[1]);
754 edgelist[1] = EKey(vIds[1], vIds[2]);
755 edgelist[2] = EKey(vIds[2], vIds[3]);
756 edgelist[3] = EKey(vIds[3], vIds[0]);
757
758 edgelist[4] = EKey(vIds[4], vIds[5]);
759 edgelist[5] = EKey(vIds[5], vIds[6]);
760 edgelist[6] = EKey(vIds[6], vIds[7]);
761 edgelist[7] = EKey(vIds[7], vIds[4]);
762
763 edgelist[8] = EKey(vIds[0], vIds[4]);
764 edgelist[9] = EKey(vIds[1], vIds[5]);
765 edgelist[10] = EKey(vIds[2], vIds[6]);
766 edgelist[11] = EKey(vIds[3], vIds[7]);
767
768 uint cubeindex;
769 cubeindex = uint(distances(i, j, k) < isoValue);
770 cubeindex += uint(distances(i + 1, j, k) < isoValue) * 2;
771 cubeindex += uint(distances(i + 1, j + 1, k) < isoValue) * 4;
772 cubeindex += uint(distances(i, j + 1, k) < isoValue) * 8;
773 cubeindex += uint(distances(i, j, k + 1) < isoValue) * 16;
774 cubeindex += uint(distances(i + 1, j, k + 1) < isoValue) * 32;
775 cubeindex += uint(distances(i + 1, j + 1, k + 1) < isoValue) * 64;
776 cubeindex += uint(distances(i, j + 1, k + 1) < isoValue) * 128;
777
778 uint numVerts = numVertsTable[cubeindex];
779
780 int index1D = i + j * (nx - 1) + k * (nx - 1) * (ny - 1);
781
782 int radix = vertNum[index1D];
783
784 EKey e[3];
785 for (int n = 0; n < numVerts; n += 3)
786 {
787 uint edge;
788 edge = triTable[cubeindex][n];
789 v[0] = vertlist[edge];
790 e[0] = edgelist[edge];
791
792 edge = triTable[cubeindex][n + 1];
793 v[1] = vertlist[edge];
794 e[1] = edgelist[edge];
795
796 edge = triTable[cubeindex][n + 2];
797 v[2] = vertlist[edge];
798 e[2] = edgelist[edge];
799
800 triangles[radix / 3] = Triangle(radix, radix + 1, radix + 2);
801
802 vertices[radix] = v[0]; edgeIds[radix++] = e[0];
803 vertices[radix] = v[1]; edgeIds[radix++] = e[1];
804 vertices[radix] = v[2]; edgeIds[radix++] = e[2];
805 }
806 }
807
808 __global__ void MCH_InitIndex(
809 DArray<uint> index)
810 {
811 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
812 if (tId >= index.size()) return;
813
814 index[tId] = tId;
815 }
816
817 template<typename EKey>
818 __global__ void MCH_CountEdgeKeys(
819 DArray<int> counter,
820 DArray<EKey> keys)
821 {
822 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
823 if (tId >= keys.size()) return;
824
825 if (tId == keys.size() - 1 || keys[tId] != keys[tId + 1])
826 counter[tId] = 1;
827 else
828 counter[tId] = 0;
829 }
830
831 template<typename Triangle>
832 __global__ void MCH_SetupNewTriangleIndex(
833 DArray<Triangle> triangles,
834 DArray<EKey> eKeys,
835 DArray<int> radix,
836 DArray<uint> ids,
837 DArray<uint> rIds)
838 {
839 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
840 if (tId >= triangles.size()) return;
841
842 Triangle t = triangles[tId];
843
844 uint v0 = rIds[t[0]];
845 uint v1 = rIds[t[1]];
846 uint v2 = rIds[t[2]];
847
848 triangles[tId] = Triangle(v0, v1, v2);
849 }
850
851 template<typename Coord>
852 __global__ void MCH_SetupNewVertices(
853 DArray<Coord> newVertices,
854 DArray<Coord> vertices,
855 DArray<uint> rIds,
856 DArray<EKey> keys,
857 DArray<int> radix,
858 DArray<uint> ids)
859 {
860 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
861 if (tId >= vertices.size()) return;
862
863 if (tId == vertices.size() - 1 || keys[tId] != keys[tId + 1]) {
864 newVertices[radix[tId]] = vertices[ids[tId]];
865 }
866
867 rIds[ids[tId]] = radix[tId];
868
869// printf("%d %d \n", tId, ids[tId]);
870 }
871
872 template<typename TDataType>
873 void MarchingCubesHelper<TDataType>::constructTriangles(
874 DArray<Coord>& vertices,
875 DArray<TopologyModule::Triangle>& triangles,
876 DArray<int>& vertNum,
877 DArray3D<Real>& distances,
878 Coord origin,
879 Real isoValue,
880 Real h)
881 {
882 DArray<EKey> eKeys(vertices.size());
883 DArray<uint> ids(vertices.size());
884 DArray<int> counter(vertices.size());
885 DArray<uint> rIds(vertices.size());
886
887 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
888 ConstructTriangles,
889 vertices,
890 eKeys,
891 triangles,
892 vertNum,
893 distances,
894 origin,
895 isoValue,
896 h);
897
898 cuExecute(ids.size(),
899 MCH_InitIndex,
900 ids);
901
902 thrust::sort_by_key(thrust::device, eKeys.begin(), eKeys.begin() + eKeys.size(), ids.begin());
903
904 cuExecute(eKeys.size(),
905 MCH_CountEdgeKeys,
906 counter,
907 eKeys);
908
909 int num = thrust::reduce(thrust::device, counter.begin(), counter.begin() + counter.size());
910 thrust::exclusive_scan(thrust::device, counter.begin(), counter.begin() + counter.size(), counter.begin());
911
912 DArray<Coord> newVertices(num);
913
914 cuExecute(vertices.size(),
915 MCH_SetupNewVertices,
916 newVertices,
917 vertices,
918 rIds,
919 eKeys,
920 counter,
921 ids);
922
923 cuExecute(triangles.size(),
924 MCH_SetupNewTriangleIndex,
925 triangles,
926 eKeys,
927 counter,
928 ids,
929 rIds);
930
931 vertices.assign(newVertices);
932
933 ids.clear();
934 rIds.clear();
935 eKeys.clear();
936 counter.clear();
937 newVertices.clear();
938 }
939
940 template<typename Real, typename Coord>
941 __global__ void MCH_CountVertexNumberForClipper(
942 DArray<int> num,
943 DArray3D<Real> distances,
944 Coord origin,
945 Real h,
946 TPlane3D<Real> plane)
947 {
948 uint i = threadIdx.x + blockDim.x * blockIdx.x;
949 uint j = threadIdx.y + blockDim.y * blockIdx.y;
950 uint k = threadIdx.z + blockDim.z * blockIdx.z;
951
952 uint nx = distances.nx();
953 uint ny = distances.ny();
954 uint nz = distances.nz();
955
956 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
957
958 Real d;
959 TPoint3D<Real> p;
960
961 Real isoValue = 0.0;
962
963 uint cubeindex;
964
965 p = origin + Coord(i * h, j * h, k * h);
966 d = p.distance(plane);
967 cubeindex = uint(d < isoValue);
968
969 p = origin + Coord((i + 1) * h, j * h, k * h);
970 d = p.distance(plane);
971 cubeindex += uint(d < isoValue) * 2;
972
973 p = origin + Coord((i + 1) * h, (j + 1) * h, k * h);
974 d = p.distance(plane);
975 cubeindex += uint(d < isoValue) * 4;
976
977 p = origin + Coord(i * h, (j + 1) * h, k * h);
978 d = p.distance(plane);
979 cubeindex += uint(d < isoValue) * 8;
980
981 p = origin + Coord(i * h, j * h, (k + 1) * h);
982 d = p.distance(plane);
983 cubeindex += uint(d < isoValue) * 16;
984
985 p = origin + Coord((i + 1) * h, j * h, (k + 1) * h);
986 d = p.distance(plane);
987 cubeindex += uint(d < isoValue) * 32;
988
989 p = origin + Coord((i + 1) * h, (j + 1) * h, (k + 1) * h);
990 d = p.distance(plane);
991 cubeindex += uint(d < isoValue) * 64;
992
993 p = origin + Coord(i * h, (j + 1) * h, (k + 1) * h);
994 d = p.distance(plane);
995 cubeindex += uint(d < isoValue) * 128;
996
997
998 num[i + j * (nx - 1) + k * (nx - 1) * (ny - 1)] = numVertsTable[cubeindex];
999 }
1000
1001 template<typename TDataType>
1002 void MarchingCubesHelper<TDataType>::countVerticeNumberForClipper(DArray<int>& num, DistanceField3D<TDataType>& sdf, TPlane3D<Real> plane)
1003 {
1004 auto& distances = sdf.distances();
1005
1006 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
1007 MCH_CountVertexNumberForClipper,
1008 num,
1009 distances,
1010 sdf.lowerBound(),
1011 sdf.getGridSpacing(),
1012 plane);
1013 }
1014
1015 template<typename Real, typename Coord, typename Triangle>
1016 __global__ void MCH_ConstructTrianglesForClipper(
1017 DArray<Real> sdfValues,
1018 DArray<Coord> vertices,
1019 DArray<Triangle> triangles,
1020 DArray<int> vertRadix,
1021 DArray3D<Real> distances,
1022 Coord origin,
1023 Real h,
1024 TPlane3D<Real> plane)
1025 {
1026 uint i = threadIdx.x + blockDim.x * blockIdx.x;
1027 uint j = threadIdx.y + blockDim.y * blockIdx.y;
1028 uint k = threadIdx.z + blockDim.z * blockIdx.z;
1029
1030 uint nx = distances.nx();
1031 uint ny = distances.ny();
1032 uint nz = distances.nz();
1033
1034 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
1035
1036 Coord v[8];
1037 Coord p = origin + Coord(i * h, j * h, k * h);
1038 v[0] = p;
1039 v[1] = p + Coord(h, 0, 0);
1040 v[2] = p + Coord(h, h, 0);
1041 v[3] = p + Coord(0, h, 0);
1042 v[4] = p + Coord(0, 0, h);
1043 v[5] = p + Coord(h, 0, h);
1044 v[6] = p + Coord(h, h, h);
1045 v[7] = p + Coord(0, h, h);
1046
1047 Real field[8];
1048 field[0] = distances(i, j, k);
1049 field[1] = distances(i + 1, j, k);
1050 field[2] = distances(i + 1, j + 1, k);
1051 field[3] = distances(i, j + 1, k);
1052 field[4] = distances(i, j, k + 1);
1053 field[5] = distances(i + 1, j, k + 1);
1054 field[6] = distances(i + 1, j + 1, k + 1);
1055 field[7] = distances(i, j + 1, k + 1);
1056
1057 Real isoValue = 0.0;
1058
1059 uint cubeindex;
1060
1061 TPoint3D<Real> pos;
1062 pos = v[0];
1063 Real d0 = pos.distance(plane);
1064 cubeindex = uint(d0 < isoValue);
1065
1066 pos = v[1];
1067 Real d1 = pos.distance(plane);
1068 cubeindex += uint(d1 < isoValue) * 2;
1069
1070 pos = v[2];
1071 Real d2 = pos.distance(plane);
1072 cubeindex += uint(d2 < isoValue) * 4;
1073
1074 pos = v[3];
1075 Real d3 = pos.distance(plane);
1076 cubeindex += uint(d3 < isoValue) * 8;
1077
1078 pos = v[4];
1079 Real d4 = pos.distance(plane);
1080 cubeindex += uint(d4 < isoValue) * 16;
1081
1082 pos = v[5];
1083 Real d5 = pos.distance(plane);
1084 cubeindex += uint(d5 < isoValue) * 32;
1085
1086 pos = v[6];
1087 Real d6 = pos.distance(plane);
1088 cubeindex += uint(d6 < isoValue) * 64;
1089
1090 pos = v[7];
1091 Real d7 = pos.distance(plane);
1092 cubeindex += uint(d7 < isoValue) * 128;
1093
1094 Real scalar[12];
1095 Coord vertlist[12];
1096 vertlist[0] = vertexInterp(scalar[0], v[0], v[1], field[0], field[1], d0, d1);
1097 vertlist[1] = vertexInterp(scalar[1], v[1], v[2], field[1], field[2], d1, d2);
1098 vertlist[2] = vertexInterp(scalar[2], v[2], v[3], field[2], field[3], d2, d3);
1099 vertlist[3] = vertexInterp(scalar[3], v[3], v[0], field[3], field[0], d3, d0);
1100
1101 vertlist[4] = vertexInterp(scalar[4], v[4], v[5], field[4], field[5], d4, d5);
1102 vertlist[5] = vertexInterp(scalar[5], v[5], v[6], field[5], field[6], d5, d6);
1103 vertlist[6] = vertexInterp(scalar[6], v[6], v[7], field[6], field[7], d6, d7);
1104 vertlist[7] = vertexInterp(scalar[7], v[7], v[4], field[7], field[4], d7, d4);
1105
1106 vertlist[8] = vertexInterp(scalar[8], v[0], v[4], field[0], field[4], d0, d4);
1107 vertlist[9] = vertexInterp(scalar[9], v[1], v[5], field[1], field[5], d1, d5);
1108 vertlist[10] = vertexInterp(scalar[10], v[2], v[6], field[2], field[6], d2, d6);
1109 vertlist[11] = vertexInterp(scalar[10], v[3], v[7], field[3], field[7], d3, d7);
1110
1111 int index1D = i + j * (nx - 1) + k * (nx - 1) * (ny - 1);
1112
1113 int radix = vertRadix[index1D];
1114
1115 uint numVerts = numVertsTable[cubeindex];
1116
1117 Real c[3];
1118 for (int n = 0; n < numVerts; n += 3)
1119 {
1120 uint edge;
1121 edge = triTable[cubeindex][n];
1122 v[0] = vertlist[edge];
1123 c[0] = scalar[edge];
1124
1125 edge = triTable[cubeindex][n + 1];
1126 v[1] = vertlist[edge];
1127 c[1] = scalar[edge];
1128
1129 edge = triTable[cubeindex][n + 2];
1130 v[2] = vertlist[edge];
1131 c[2] = scalar[edge];
1132
1133 triangles[radix / 3] = Triangle(radix, radix + 1, radix + 2);
1134
1135 vertices[radix] = v[0]; sdfValues[radix] = c[0]; radix++;
1136 vertices[radix] = v[1]; sdfValues[radix] = c[1]; radix++;
1137 vertices[radix] = v[2]; sdfValues[radix] = c[2]; radix++;
1138 }
1139 }
1140
1141 template<typename TDataType>
1142 void MarchingCubesHelper<TDataType>::constructTrianglesForClipper(
1143 DArray<Real>& field,
1144 DArray<Coord>& vertices,
1145 DArray<TopologyModule::Triangle>& triangles,
1146 DArray<int>& vertNum,
1147 DistanceField3D<TDataType>& sdf,
1148 TPlane3D<Real> plane)
1149 {
1150 auto& distances = sdf.distances();
1151
1152 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
1153 MCH_ConstructTrianglesForClipper,
1154 field,
1155 vertices,
1156 triangles,
1157 vertNum,
1158 distances,
1159 sdf.lowerBound(),
1160 sdf.getGridSpacing(),
1161 plane);
1162 }
1163
1164 template<typename Real, typename Coord>
1165 __global__ void CountVerticeNumberForOctree(
1166 DArray<uint> num,
1167 DArray<Coord> vertices,
1168 DArray<Real> sdfs,
1169 Real isoValue)
1170 {
1171 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1172 if (pId >= num.size()) return;
1173
1174 int vIdx = 8 * pId;
1175
1176 uint cubeindex;
1177 cubeindex = uint(sdfs[vIdx] < isoValue);
1178 cubeindex += uint(sdfs[vIdx + 1] < isoValue) * 2;
1179 cubeindex += uint(sdfs[vIdx + 2] < isoValue) * 4;
1180 cubeindex += uint(sdfs[vIdx + 3] < isoValue) * 8;
1181 cubeindex += uint(sdfs[vIdx + 4] < isoValue) * 16;
1182 cubeindex += uint(sdfs[vIdx + 5] < isoValue) * 32;
1183 cubeindex += uint(sdfs[vIdx + 6] < isoValue) * 64;
1184 cubeindex += uint(sdfs[vIdx + 7] < isoValue) * 128;
1185
1186 num[pId] = numVertsTable[cubeindex];
1187 }
1188
1189 template<typename TDataType>
1190 void MarchingCubesHelper<TDataType>::countVerticeNumberForOctree(
1191 DArray<uint>& num,
1192 DArray<Coord>& vertices,
1193 DArray<Real>& sdfs,
1194 Real isoValue)
1195 {
1196 cuExecute(num.size(),
1197 CountVerticeNumberForOctree,
1198 num,
1199 vertices,
1200 sdfs,
1201 isoValue);
1202 }
1203
1204 template<typename Real, typename Coord>
1205 __global__ void ConstructTrianglesForOctree(
1206 DArray<Coord> triangleVertices,
1207 DArray<TopologyModule::Triangle> triangles,
1208 DArray<uint> vertNum,
1209 DArray<Coord> cellVertices,
1210 DArray<Real> sdfs,
1211 Real isoValue)
1212 {
1213 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1214 if (pId >= vertNum.size()) return;
1215
1216 int vIdx = 8 * pId;
1217
1218 Coord v[8];
1219 v[0] = cellVertices[vIdx];
1220 v[1] = cellVertices[vIdx + 1];
1221 v[2] = cellVertices[vIdx + 2];
1222 v[3] = cellVertices[vIdx + 3];
1223 v[4] = cellVertices[vIdx + 4];
1224 v[5] = cellVertices[vIdx + 5];
1225 v[6] = cellVertices[vIdx + 6];
1226 v[7] = cellVertices[vIdx + 7];
1227
1228 Real field[8];
1229 field[0] = sdfs[vIdx];
1230 field[1] = sdfs[vIdx + 1];
1231 field[2] = sdfs[vIdx + 2];
1232 field[3] = sdfs[vIdx + 3];
1233 field[4] = sdfs[vIdx + 4];
1234 field[5] = sdfs[vIdx + 5];
1235 field[6] = sdfs[vIdx + 6];
1236 field[7] = sdfs[vIdx + 7];
1237
1238 Coord vertlist[12];
1239 vertlist[0] = vertexInterp(isoValue, v[0], v[1], field[0], field[1]);
1240 vertlist[1] = vertexInterp(isoValue, v[1], v[2], field[1], field[2]);
1241 vertlist[2] = vertexInterp(isoValue, v[2], v[3], field[2], field[3]);
1242 vertlist[3] = vertexInterp(isoValue, v[3], v[0], field[3], field[0]);
1243
1244 vertlist[4] = vertexInterp(isoValue, v[4], v[5], field[4], field[5]);
1245 vertlist[5] = vertexInterp(isoValue, v[5], v[6], field[5], field[6]);
1246 vertlist[6] = vertexInterp(isoValue, v[6], v[7], field[6], field[7]);
1247 vertlist[7] = vertexInterp(isoValue, v[7], v[4], field[7], field[4]);
1248
1249 vertlist[8] = vertexInterp(isoValue, v[0], v[4], field[0], field[4]);
1250 vertlist[9] = vertexInterp(isoValue, v[1], v[5], field[1], field[5]);
1251 vertlist[10] = vertexInterp(isoValue, v[2], v[6], field[2], field[6]);
1252 vertlist[11] = vertexInterp(isoValue, v[3], v[7], field[3], field[7]);
1253
1254 uint cubeindex;
1255 cubeindex = uint(sdfs[vIdx] < isoValue);
1256 cubeindex += uint(sdfs[vIdx + 1] < isoValue) * 2;
1257 cubeindex += uint(sdfs[vIdx + 2] < isoValue) * 4;
1258 cubeindex += uint(sdfs[vIdx + 3] < isoValue) * 8;
1259 cubeindex += uint(sdfs[vIdx + 4] < isoValue) * 16;
1260 cubeindex += uint(sdfs[vIdx + 5] < isoValue) * 32;
1261 cubeindex += uint(sdfs[vIdx + 6] < isoValue) * 64;
1262 cubeindex += uint(sdfs[vIdx + 7] < isoValue) * 128;
1263
1264 uint numVerts = numVertsTable[cubeindex];
1265
1266 int radix = vertNum[pId];
1267
1268 for (int n = 0; n < numVerts; n += 3)
1269 {
1270 uint edge;
1271 edge = triTable[cubeindex][n];
1272 v[0] = vertlist[edge];
1273
1274 edge = triTable[cubeindex][n + 1];
1275 v[1] = vertlist[edge];
1276
1277 edge = triTable[cubeindex][n + 2];
1278 v[2] = vertlist[edge];
1279
1280 triangles[radix / 3] = TopologyModule::Triangle(radix, radix + 1, radix + 2);
1281
1282 triangleVertices[radix++] = v[0];
1283 triangleVertices[radix++] = v[1];
1284 triangleVertices[radix++] = v[2];
1285 }
1286 }
1287
1288 template<typename TDataType>
1289 void MarchingCubesHelper<TDataType>::constructTrianglesForOctree(
1290 DArray<Coord>& triangleVertices,
1291 DArray<TopologyModule::Triangle>& triangles,
1292 DArray<uint>& num,
1293 DArray<Coord>& cellVertices,
1294 DArray<Real>& sdfs,
1295 Real isoValue)
1296 {
1297 cuExecute(num.size(),
1298 ConstructTrianglesForOctree,
1299 triangleVertices,
1300 triangles,
1301 num,
1302 cellVertices,
1303 sdfs,
1304 isoValue);
1305 }
1306
1307 template<typename Real, typename Coord>
1308 __global__ void MCH_CountVertexNumberForOctreeClipper(
1309 DArray<uint> num,
1310 DArray<Coord> vertices,
1311 TPlane3D<Real> plane)
1312 {
1313 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1314 if (pId >= num.size()) return;
1315
1316 int vIdx = 8 * pId;
1317
1318 Real d;
1319 TPoint3D<Real> p;
1320
1321 Real isoValue = 0.0;
1322
1323 uint cubeindex;
1324
1325 p = vertices[vIdx];
1326 d = p.distance(plane);
1327 cubeindex = uint(d < isoValue);
1328
1329 p = vertices[vIdx + 1];
1330 d = p.distance(plane);
1331 cubeindex += uint(d < isoValue) * 2;
1332
1333 p = vertices[vIdx + 2];
1334 d = p.distance(plane);
1335 cubeindex += uint(d < isoValue) * 4;
1336
1337 p = vertices[vIdx + 3];
1338 d = p.distance(plane);
1339 cubeindex += uint(d < isoValue) * 8;
1340
1341 p = vertices[vIdx + 4];
1342 d = p.distance(plane);
1343 cubeindex += uint(d < isoValue) * 16;
1344
1345 p = vertices[vIdx + 5];
1346 d = p.distance(plane);
1347 cubeindex += uint(d < isoValue) * 32;
1348
1349 p = vertices[vIdx + 6];
1350 d = p.distance(plane);
1351 cubeindex += uint(d < isoValue) * 64;
1352
1353 p = vertices[vIdx + 7];
1354 d = p.distance(plane);
1355 cubeindex += uint(d < isoValue) * 128;
1356
1357 num[pId] = numVertsTable[cubeindex];
1358 }
1359
1360 template<typename TDataType>
1361 void MarchingCubesHelper<TDataType>::countVerticeNumberForOctreeClipper(
1362 DArray<uint>& num,
1363 DArray<Coord>& vertices,
1364 TPlane3D<Real> plane)
1365 {
1366 cuExecute(num.size(),
1367 MCH_CountVertexNumberForOctreeClipper,
1368 num,
1369 vertices,
1370 plane);
1371 }
1372
1373 template<typename Real, typename Coord>
1374 __global__ void MCH_ConstructTrianglesForOctreeClipper(
1375 DArray<Real> vertSDFs,
1376 DArray<Coord> triangleVertices,
1377 DArray<TopologyModule::Triangle> triangles,
1378 DArray<uint> vertNum,
1379 DArray<Coord> cellVertices,
1380 DArray<Real> sdfs,
1381 TPlane3D<Real> plane)
1382 {
1383 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1384 if (pId >= vertNum.size()) return;
1385
1386 int vIdx = 8 * pId;
1387
1388 Coord v[8];
1389 v[0] = cellVertices[vIdx];
1390 v[1] = cellVertices[vIdx + 1];
1391 v[2] = cellVertices[vIdx + 2];
1392 v[3] = cellVertices[vIdx + 3];
1393 v[4] = cellVertices[vIdx + 4];
1394 v[5] = cellVertices[vIdx + 5];
1395 v[6] = cellVertices[vIdx + 6];
1396 v[7] = cellVertices[vIdx + 7];
1397
1398 Real isoValue = Real(0);
1399
1400 Real field[8];
1401 field[0] = sdfs[vIdx];
1402 field[1] = sdfs[vIdx + 1];
1403 field[2] = sdfs[vIdx + 2];
1404 field[3] = sdfs[vIdx + 3];
1405 field[4] = sdfs[vIdx + 4];
1406 field[5] = sdfs[vIdx + 5];
1407 field[6] = sdfs[vIdx + 6];
1408 field[7] = sdfs[vIdx + 7];
1409
1410 uint cubeindex;
1411 TPoint3D<Real> pos;
1412 pos = cellVertices[vIdx];
1413 Real d0 = pos.distance(plane);
1414 cubeindex = uint(d0 < isoValue);
1415
1416 pos = cellVertices[vIdx + 1];
1417 Real d1 = pos.distance(plane);
1418 cubeindex += uint(d1 < isoValue) * 2;
1419
1420 pos = cellVertices[vIdx + 2];
1421 Real d2 = pos.distance(plane);
1422 cubeindex += uint(d2 < isoValue) * 4;
1423
1424 pos = cellVertices[vIdx + 3];
1425 Real d3 = pos.distance(plane);
1426 cubeindex += uint(d3 < isoValue) * 8;
1427
1428 pos = cellVertices[vIdx + 4];
1429 Real d4 = pos.distance(plane);
1430 cubeindex += uint(d4 < isoValue) * 16;
1431
1432 pos = cellVertices[vIdx + 5];
1433 Real d5 = pos.distance(plane);
1434 cubeindex += uint(d5 < isoValue) * 32;
1435
1436 pos = cellVertices[vIdx + 6];
1437 Real d6 = pos.distance(plane);
1438 cubeindex += uint(d6 < isoValue) * 64;
1439
1440 pos = cellVertices[vIdx + 7];
1441 Real d7 = pos.distance(plane);
1442 cubeindex += uint(d7 < isoValue) * 128;
1443
1444 Real scalar[12];
1445 Coord vertlist[12];
1446 vertlist[0] = vertexInterp(scalar[0], v[0], v[1], field[0], field[1], d0, d1);
1447 vertlist[1] = vertexInterp(scalar[1], v[1], v[2], field[1], field[2], d1, d2);
1448 vertlist[2] = vertexInterp(scalar[2], v[2], v[3], field[2], field[3], d2, d3);
1449 vertlist[3] = vertexInterp(scalar[3], v[3], v[0], field[3], field[0], d3, d0);
1450
1451 vertlist[4] = vertexInterp(scalar[4], v[4], v[5], field[4], field[5], d4, d5);
1452 vertlist[5] = vertexInterp(scalar[5], v[5], v[6], field[5], field[6], d5, d6);
1453 vertlist[6] = vertexInterp(scalar[6], v[6], v[7], field[6], field[7], d6, d7);
1454 vertlist[7] = vertexInterp(scalar[7], v[7], v[4], field[7], field[4], d7, d4);
1455
1456 vertlist[8] = vertexInterp(scalar[8], v[0], v[4], field[0], field[4], d0, d4);
1457 vertlist[9] = vertexInterp(scalar[9], v[1], v[5], field[1], field[5], d1, d5);
1458 vertlist[10] = vertexInterp(scalar[10], v[2], v[6], field[2], field[6], d2, d6);
1459 vertlist[11] = vertexInterp(scalar[10], v[3], v[7], field[3], field[7], d3, d7);
1460
1461 uint numVerts = numVertsTable[cubeindex];
1462
1463 int radix = vertNum[pId];
1464
1465 Real c[3];
1466 for (int n = 0; n < numVerts; n += 3)
1467 {
1468 uint edge;
1469 edge = triTable[cubeindex][n];
1470 v[0] = vertlist[edge];
1471 c[0] = scalar[edge];
1472
1473 edge = triTable[cubeindex][n + 1];
1474 v[1] = vertlist[edge];
1475 c[1] = scalar[edge];
1476
1477 edge = triTable[cubeindex][n + 2];
1478 v[2] = vertlist[edge];
1479 c[2] = scalar[edge];
1480
1481 triangles[radix / 3] = TopologyModule::Triangle(radix, radix + 1, radix + 2);
1482
1483 triangleVertices[radix] = v[0]; vertSDFs[radix] = c[0]; radix++;
1484 triangleVertices[radix] = v[1]; vertSDFs[radix] = c[1]; radix++;
1485 triangleVertices[radix] = v[2]; vertSDFs[radix] = c[2]; radix++;
1486 }
1487 }
1488
1489 template<typename TDataType>
1490 void MarchingCubesHelper<TDataType>::constructTrianglesForOctreeClipper(
1491 DArray<Real>& vertSDFs,
1492 DArray<Coord>& triangleVertices,
1493 DArray<TopologyModule::Triangle>& triangles,
1494 DArray<uint>& num,
1495 DArray<Coord>& cellVertices,
1496 DArray<Real>& sdfs,
1497 TPlane3D<Real> plane)
1498 {
1499 cuExecute(num.size(),
1500 MCH_ConstructTrianglesForOctreeClipper,
1501 vertSDFs,
1502 triangleVertices,
1503 triangles,
1504 num,
1505 cellVertices,
1506 sdfs,
1507 plane);
1508 }
1509
1510 DEFINE_CLASS(MarchingCubesHelper);
1511}