1#include "MarchingCubesHelper.h"
3#include "Topology/EdgeSet.h"
5#include <thrust/sort.h>
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
47 int triTable[256][16] =
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}
308 uint numVertsTable[256] =
568 template<typename Real, typename Coord, typename TDataType>
569 __global__ void RconstructSDF(
570 DArray3D<Real> distances,
573 DistanceField3D<TDataType> sdf)
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;
579 uint nx = distances.nx();
580 uint ny = distances.ny();
581 uint nz = distances.nz();
583 if (i >= nx || j >= ny || k >= nz) return;
585 Coord p = origin + Coord(i * h, j * h, k * h);
589 sdf.getDistance(p, d, normal);
590 distances(i, j, k) = d;
593 template<typename TDataType>
594 void MarchingCubesHelper<TDataType>::reconstructSDF(
595 DArray3D<Real>& distances,
598 DistanceField3D<TDataType>& sdf)
600 cuExecute3D(make_uint3(distances.nx(), distances.ny(), distances.nz()),
608 template<typename Real>
609 __global__ void CountVerticeNumber(
611 DArray3D<Real> distances,
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;
618 uint nx = distances.nx();
619 uint ny = distances.ny();
620 uint nz = distances.nz();
622 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
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;
634 num[i + j * (nx - 1) + k * (nx - 1) * (ny - 1)] = numVertsTable[cubeindex];
637 template<typename TDataType>
638 void MarchingCubesHelper<TDataType>::countVerticeNumber(
640 DArray3D<Real>& distances,
643 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
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)
654 if (abs(f1 - f0) < EPSILON)
655 return 0.5 * (p0 + p1);
657 Real t = (isolevel - f0) / (f1 - f0);
659 t = clamp(t, Real(0), Real(1));
661 return (1 - t) * p0 + t * p1;
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)
668 if (abs(d0 - d1) < EPSILON)
670 value = 0.5 * (f0 + f1);
671 return 0.5 * (p0 + p1);
675 Real t = d0 / (d0 - d1);
677 t = clamp(t, Real(0), Real(1));
679 value = (1 - t) * f0 + t * f1;
681 return (1 - t) * p0 + t * p1;
684 template<typename Coord, typename Triangle>
685 __global__ void ConstructTriangles(
686 DArray<Coord> vertices,
687 DArray<EKey> edgeIds,
688 DArray<Triangle> triangles,
690 DArray3D<Real> distances,
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;
699 uint nx = distances.nx();
700 uint ny = distances.ny();
701 uint nz = distances.nz();
703 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
706 Coord p = origin + h * Coord(i, j, k);
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);
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);
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);
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]);
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]);
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]);
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]);
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]);
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]);
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;
778 uint numVerts = numVertsTable[cubeindex];
780 int index1D = i + j * (nx - 1) + k * (nx - 1) * (ny - 1);
782 int radix = vertNum[index1D];
785 for (int n = 0; n < numVerts; n += 3)
788 edge = triTable[cubeindex][n];
789 v[0] = vertlist[edge];
790 e[0] = edgelist[edge];
792 edge = triTable[cubeindex][n + 1];
793 v[1] = vertlist[edge];
794 e[1] = edgelist[edge];
796 edge = triTable[cubeindex][n + 2];
797 v[2] = vertlist[edge];
798 e[2] = edgelist[edge];
800 triangles[radix / 3] = Triangle(radix, radix + 1, radix + 2);
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];
808 __global__ void MCH_InitIndex(
811 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
812 if (tId >= index.size()) return;
817 template<typename EKey>
818 __global__ void MCH_CountEdgeKeys(
822 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
823 if (tId >= keys.size()) return;
825 if (tId == keys.size() - 1 || keys[tId] != keys[tId + 1])
831 template<typename Triangle>
832 __global__ void MCH_SetupNewTriangleIndex(
833 DArray<Triangle> triangles,
839 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
840 if (tId >= triangles.size()) return;
842 Triangle t = triangles[tId];
844 uint v0 = rIds[t[0]];
845 uint v1 = rIds[t[1]];
846 uint v2 = rIds[t[2]];
848 triangles[tId] = Triangle(v0, v1, v2);
851 template<typename Coord>
852 __global__ void MCH_SetupNewVertices(
853 DArray<Coord> newVertices,
854 DArray<Coord> vertices,
860 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
861 if (tId >= vertices.size()) return;
863 if (tId == vertices.size() - 1 || keys[tId] != keys[tId + 1]) {
864 newVertices[radix[tId]] = vertices[ids[tId]];
867 rIds[ids[tId]] = radix[tId];
869// printf("%d %d \n", tId, ids[tId]);
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,
882 DArray<EKey> eKeys(vertices.size());
883 DArray<uint> ids(vertices.size());
884 DArray<int> counter(vertices.size());
885 DArray<uint> rIds(vertices.size());
887 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
898 cuExecute(ids.size(),
902 thrust::sort_by_key(thrust::device, eKeys.begin(), eKeys.begin() + eKeys.size(), ids.begin());
904 cuExecute(eKeys.size(),
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());
912 DArray<Coord> newVertices(num);
914 cuExecute(vertices.size(),
915 MCH_SetupNewVertices,
923 cuExecute(triangles.size(),
924 MCH_SetupNewTriangleIndex,
931 vertices.assign(newVertices);
940 template<typename Real, typename Coord>
941 __global__ void MCH_CountVertexNumberForClipper(
943 DArray3D<Real> distances,
946 TPlane3D<Real> plane)
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;
952 uint nx = distances.nx();
953 uint ny = distances.ny();
954 uint nz = distances.nz();
956 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
965 p = origin + Coord(i * h, j * h, k * h);
966 d = p.distance(plane);
967 cubeindex = uint(d < isoValue);
969 p = origin + Coord((i + 1) * h, j * h, k * h);
970 d = p.distance(plane);
971 cubeindex += uint(d < isoValue) * 2;
973 p = origin + Coord((i + 1) * h, (j + 1) * h, k * h);
974 d = p.distance(plane);
975 cubeindex += uint(d < isoValue) * 4;
977 p = origin + Coord(i * h, (j + 1) * h, k * h);
978 d = p.distance(plane);
979 cubeindex += uint(d < isoValue) * 8;
981 p = origin + Coord(i * h, j * h, (k + 1) * h);
982 d = p.distance(plane);
983 cubeindex += uint(d < isoValue) * 16;
985 p = origin + Coord((i + 1) * h, j * h, (k + 1) * h);
986 d = p.distance(plane);
987 cubeindex += uint(d < isoValue) * 32;
989 p = origin + Coord((i + 1) * h, (j + 1) * h, (k + 1) * h);
990 d = p.distance(plane);
991 cubeindex += uint(d < isoValue) * 64;
993 p = origin + Coord(i * h, (j + 1) * h, (k + 1) * h);
994 d = p.distance(plane);
995 cubeindex += uint(d < isoValue) * 128;
998 num[i + j * (nx - 1) + k * (nx - 1) * (ny - 1)] = numVertsTable[cubeindex];
1001 template<typename TDataType>
1002 void MarchingCubesHelper<TDataType>::countVerticeNumberForClipper(DArray<int>& num, DistanceField3D<TDataType>& sdf, TPlane3D<Real> plane)
1004 auto& distances = sdf.distances();
1006 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
1007 MCH_CountVertexNumberForClipper,
1011 sdf.getGridSpacing(),
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,
1024 TPlane3D<Real> plane)
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;
1030 uint nx = distances.nx();
1031 uint ny = distances.ny();
1032 uint nz = distances.nz();
1034 if (i >= nx - 1 || j >= ny - 1 || k >= nz - 1) return;
1037 Coord p = origin + Coord(i * h, j * h, k * h);
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);
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);
1057 Real isoValue = 0.0;
1063 Real d0 = pos.distance(plane);
1064 cubeindex = uint(d0 < isoValue);
1067 Real d1 = pos.distance(plane);
1068 cubeindex += uint(d1 < isoValue) * 2;
1071 Real d2 = pos.distance(plane);
1072 cubeindex += uint(d2 < isoValue) * 4;
1075 Real d3 = pos.distance(plane);
1076 cubeindex += uint(d3 < isoValue) * 8;
1079 Real d4 = pos.distance(plane);
1080 cubeindex += uint(d4 < isoValue) * 16;
1083 Real d5 = pos.distance(plane);
1084 cubeindex += uint(d5 < isoValue) * 32;
1087 Real d6 = pos.distance(plane);
1088 cubeindex += uint(d6 < isoValue) * 64;
1091 Real d7 = pos.distance(plane);
1092 cubeindex += uint(d7 < isoValue) * 128;
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);
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);
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);
1111 int index1D = i + j * (nx - 1) + k * (nx - 1) * (ny - 1);
1113 int radix = vertRadix[index1D];
1115 uint numVerts = numVertsTable[cubeindex];
1118 for (int n = 0; n < numVerts; n += 3)
1121 edge = triTable[cubeindex][n];
1122 v[0] = vertlist[edge];
1123 c[0] = scalar[edge];
1125 edge = triTable[cubeindex][n + 1];
1126 v[1] = vertlist[edge];
1127 c[1] = scalar[edge];
1129 edge = triTable[cubeindex][n + 2];
1130 v[2] = vertlist[edge];
1131 c[2] = scalar[edge];
1133 triangles[radix / 3] = Triangle(radix, radix + 1, radix + 2);
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++;
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)
1150 auto& distances = sdf.distances();
1152 cuExecute3D(make_uint3(distances.nx() - 1, distances.ny() - 1, distances.nz() - 1),
1153 MCH_ConstructTrianglesForClipper,
1160 sdf.getGridSpacing(),
1164 template<typename Real, typename Coord>
1165 __global__ void CountVerticeNumberForOctree(
1167 DArray<Coord> vertices,
1171 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1172 if (pId >= num.size()) return;
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;
1186 num[pId] = numVertsTable[cubeindex];
1189 template<typename TDataType>
1190 void MarchingCubesHelper<TDataType>::countVerticeNumberForOctree(
1192 DArray<Coord>& vertices,
1196 cuExecute(num.size(),
1197 CountVerticeNumberForOctree,
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,
1213 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1214 if (pId >= vertNum.size()) return;
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];
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];
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]);
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]);
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]);
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;
1264 uint numVerts = numVertsTable[cubeindex];
1266 int radix = vertNum[pId];
1268 for (int n = 0; n < numVerts; n += 3)
1271 edge = triTable[cubeindex][n];
1272 v[0] = vertlist[edge];
1274 edge = triTable[cubeindex][n + 1];
1275 v[1] = vertlist[edge];
1277 edge = triTable[cubeindex][n + 2];
1278 v[2] = vertlist[edge];
1280 triangles[radix / 3] = TopologyModule::Triangle(radix, radix + 1, radix + 2);
1282 triangleVertices[radix++] = v[0];
1283 triangleVertices[radix++] = v[1];
1284 triangleVertices[radix++] = v[2];
1288 template<typename TDataType>
1289 void MarchingCubesHelper<TDataType>::constructTrianglesForOctree(
1290 DArray<Coord>& triangleVertices,
1291 DArray<TopologyModule::Triangle>& triangles,
1293 DArray<Coord>& cellVertices,
1297 cuExecute(num.size(),
1298 ConstructTrianglesForOctree,
1307 template<typename Real, typename Coord>
1308 __global__ void MCH_CountVertexNumberForOctreeClipper(
1310 DArray<Coord> vertices,
1311 TPlane3D<Real> plane)
1313 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1314 if (pId >= num.size()) return;
1321 Real isoValue = 0.0;
1326 d = p.distance(plane);
1327 cubeindex = uint(d < isoValue);
1329 p = vertices[vIdx + 1];
1330 d = p.distance(plane);
1331 cubeindex += uint(d < isoValue) * 2;
1333 p = vertices[vIdx + 2];
1334 d = p.distance(plane);
1335 cubeindex += uint(d < isoValue) * 4;
1337 p = vertices[vIdx + 3];
1338 d = p.distance(plane);
1339 cubeindex += uint(d < isoValue) * 8;
1341 p = vertices[vIdx + 4];
1342 d = p.distance(plane);
1343 cubeindex += uint(d < isoValue) * 16;
1345 p = vertices[vIdx + 5];
1346 d = p.distance(plane);
1347 cubeindex += uint(d < isoValue) * 32;
1349 p = vertices[vIdx + 6];
1350 d = p.distance(plane);
1351 cubeindex += uint(d < isoValue) * 64;
1353 p = vertices[vIdx + 7];
1354 d = p.distance(plane);
1355 cubeindex += uint(d < isoValue) * 128;
1357 num[pId] = numVertsTable[cubeindex];
1360 template<typename TDataType>
1361 void MarchingCubesHelper<TDataType>::countVerticeNumberForOctreeClipper(
1363 DArray<Coord>& vertices,
1364 TPlane3D<Real> plane)
1366 cuExecute(num.size(),
1367 MCH_CountVertexNumberForOctreeClipper,
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,
1381 TPlane3D<Real> plane)
1383 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
1384 if (pId >= vertNum.size()) return;
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];
1398 Real isoValue = Real(0);
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];
1412 pos = cellVertices[vIdx];
1413 Real d0 = pos.distance(plane);
1414 cubeindex = uint(d0 < isoValue);
1416 pos = cellVertices[vIdx + 1];
1417 Real d1 = pos.distance(plane);
1418 cubeindex += uint(d1 < isoValue) * 2;
1420 pos = cellVertices[vIdx + 2];
1421 Real d2 = pos.distance(plane);
1422 cubeindex += uint(d2 < isoValue) * 4;
1424 pos = cellVertices[vIdx + 3];
1425 Real d3 = pos.distance(plane);
1426 cubeindex += uint(d3 < isoValue) * 8;
1428 pos = cellVertices[vIdx + 4];
1429 Real d4 = pos.distance(plane);
1430 cubeindex += uint(d4 < isoValue) * 16;
1432 pos = cellVertices[vIdx + 5];
1433 Real d5 = pos.distance(plane);
1434 cubeindex += uint(d5 < isoValue) * 32;
1436 pos = cellVertices[vIdx + 6];
1437 Real d6 = pos.distance(plane);
1438 cubeindex += uint(d6 < isoValue) * 64;
1440 pos = cellVertices[vIdx + 7];
1441 Real d7 = pos.distance(plane);
1442 cubeindex += uint(d7 < isoValue) * 128;
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);
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);
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);
1461 uint numVerts = numVertsTable[cubeindex];
1463 int radix = vertNum[pId];
1466 for (int n = 0; n < numVerts; n += 3)
1469 edge = triTable[cubeindex][n];
1470 v[0] = vertlist[edge];
1471 c[0] = scalar[edge];
1473 edge = triTable[cubeindex][n + 1];
1474 v[1] = vertlist[edge];
1475 c[1] = scalar[edge];
1477 edge = triTable[cubeindex][n + 2];
1478 v[2] = vertlist[edge];
1479 c[2] = scalar[edge];
1481 triangles[radix / 3] = TopologyModule::Triangle(radix, radix + 1, radix + 2);
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++;
1489 template<typename TDataType>
1490 void MarchingCubesHelper<TDataType>::constructTrianglesForOctreeClipper(
1491 DArray<Real>& vertSDFs,
1492 DArray<Coord>& triangleVertices,
1493 DArray<TopologyModule::Triangle>& triangles,
1495 DArray<Coord>& cellVertices,
1497 TPlane3D<Real> plane)
1499 cuExecute(num.size(),
1500 MCH_ConstructTrianglesForOctreeClipper,
1510 DEFINE_CLASS(MarchingCubesHelper);