MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
iso_surface.cc
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 *
9 * The code assumes an ordering of cube vertices and cube edges as follows:
10 *
11 * 2------3 +---3--+ +------+
12 * /| /| 5| 11| /| /| y
13 * 6-+----7 | +-+ 9--+ | +-+----+ | |
14 * | | | | | | | | | 4 | 1 |
15 * | 0----+-1 | +--0-+-+ | +----+-+ +------ x
16 * |/ |/ |8 |2 7/ 10/ /
17 * 4------5 +--6---+ +------+ z
18 * Vertex Order Edge Order 1 Edge Order 2
19 */
20
21#include <iostream>
22#include <bitset>
23
24#include "util/timer.h"
25#include "fssr/octree.h"
26#include "fssr/iso_surface.h"
27#include "fssr/triangulation.h"
28
29#define ISO_VALUE 0.0f
30#define CUBE_CORNERS 8
31#define CUBE_EDGES 12
32#define CUBE_FACES 6
33
35
37namespace
38{
40 enum CubeFace
41 {
42 POSITIVE_X = 0,
43 NEGATIVE_X = 1,
44 POSITIVE_Y = 2,
45 NEGATIVE_Y = 3,
46 POSITIVE_Z = 4,
47 NEGATIVE_Z = 5
48 };
49
50 CubeFace const face_opposite[6] =
51 {
52 NEGATIVE_X, POSITIVE_X,
53 NEGATIVE_Y, POSITIVE_Y,
54 NEGATIVE_Z, POSITIVE_Z
55 };
56
58 int const face_corners[6][4] =
59 {
60 { 1, 3, 5, 7 }, /* Positive X. */
61 { 0, 2, 4, 6 }, /* Negative X. */
62 { 2, 3, 6, 7 }, /* Positive Y. */
63 { 0, 1, 4, 5 }, /* Negative Y. */
64 { 4, 5, 6, 7 }, /* Positive Z. */
65 { 0, 1, 2, 3 }, /* Negative Z. */
66 };
67
68 bool const is_edge_on_face[12][6] =
69 {
70 /* PosX NegX PosY NegY PosZ NegZ */
71 { false, false, false, true , false, true }, /* Edge 0. */
72 { true , false, false, false, false, true }, /* Edge 1. */
73 { true , false, false, true , false, false }, /* Edge 2. */
74 { false, false, true , false, false, true }, /* Edge 3. */
75 { false, true , false, false, false, true }, /* Edge 4. */
76 { false, true , true , false, false, false }, /* Edge 5. */
77 { false, false, false, true , true , false }, /* Edge 6. */
78 { false, true , false, false, true , false }, /* Edge 7. */
79 { false, true , false, true , false, false }, /* Edge 8. */
80 { false, false, true , false, true , false }, /* Edge 9. */
81 { true , false, false, false, true , false }, /* Edge 10. */
82 { true , false, true , false, false, false }, /* Edge 11. */
83 };
84
86 CubeFace const edge_neighbors[12][2] =
87 {
88 { NEGATIVE_Y, NEGATIVE_Z }, /* Edge 0. */
89 { POSITIVE_X, NEGATIVE_Z }, /* Edge 1. */
90 { POSITIVE_X, NEGATIVE_Y }, /* Edge 2. */
91 { POSITIVE_Y, NEGATIVE_Z }, /* Edge 3. */
92 { NEGATIVE_X, NEGATIVE_Z }, /* Edge 4. */
93 { NEGATIVE_X, POSITIVE_Y }, /* Edge 5. */
94 { NEGATIVE_Y, POSITIVE_Z }, /* Edge 6. */
95 { NEGATIVE_X, POSITIVE_Z }, /* Edge 7. */
96 { NEGATIVE_X, NEGATIVE_Y }, /* Edge 8. */
97 { POSITIVE_Y, POSITIVE_Z }, /* Edge 9. */
98 { POSITIVE_X, POSITIVE_Z }, /* Edge 10. */
99 { POSITIVE_X, POSITIVE_Y } /* Edge 11. */
100 };
101
103 int const edge_corners[12][2] =
104 {
105 { 0, 1 }, /* Edge 0. */
106 { 1, 3 }, /* Edge 1. */
107 { 1, 5 }, /* Edge 2. */
108 { 2, 3 }, /* Edge 3. */
109 { 0, 2 }, /* Edge 4. */
110 { 2, 6 }, /* Edge 5. */
111 { 4, 5 }, /* Edge 6. */
112 { 4, 6 }, /* Edge 7. */
113 { 0, 4 }, /* Edge 8. */
114 { 6, 7 }, /* Edge 9. */
115 { 5, 7 }, /* Edge 10. */
116 { 3, 7 }, /* Edge 11. */
117 };
118
120 int const edge_reflections[12][3] =
121 {
122 { 3, 6, 9 }, /* Edge 0. */
123 { 4, 10, 7 }, /* Edge 1. */
124 { 8, 11, 5 }, /* Edge 2. */
125 { 0, 9, 6 }, /* Edge 3. */
126 { 1, 7, 10 }, /* Edge 4. */
127 { 11, 8, 2 }, /* Edge 5. */
128 { 9, 0, 3 }, /* Edge 6. */
129 { 10, 4, 1 }, /* Edge 7. */
130 { 2, 5, 11 }, /* Edge 8. */
131 { 6, 3, 0 }, /* Edge 9. */
132 { 7, 1, 4 }, /* Edge 10. */
133 { 5, 2, 8 }, /* Edge 11. */
134 };
135
137 int const mc_polygons[256][17] =
138 {
139 { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
140 { 0, 4, 8, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
141 { 0, 2, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
142 { 4, 8, 2, 1, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
143 { 3, 5, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
144 { 8, 0, 3, 5, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
145 { 0, 2, 1, 0, 4, 3, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
146 { 2, 1, 3, 5, 8, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
147 { 1, 11, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
148 { 1, 11, 3, 1, 0, 4, 8, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
149 { 11, 3, 0, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
150 { 11, 3, 4, 8, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
151 { 5, 4, 1, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
152 { 8, 0, 1, 11, 5, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
153 { 5, 4, 0, 2, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
154 { 2, 11, 5, 8, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
155 { 8, 7, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
156 { 7, 6, 0, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
157 { 2, 1, 0, 2, 6, 8, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
158 { 7, 6, 2, 1, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
159 { 4, 3, 5, 4, 8, 7, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
160 { 3, 5, 7, 6, 0, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
161 { 1, 0, 2, 1, 4, 3, 5, 4, 6, 8, 7, 6, -1, -1, -1, -1, -1 },
162 { 3, 5, 7, 6, 2, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
163 { 1, 11, 3, 1, 6, 8, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
164 { 7, 6, 0, 4, 7, 3, 1, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1 },
165 { 11, 3, 0, 2, 11, 6, 8, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1 },
166 { 6, 2, 11, 3, 4, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
167 { 1, 11, 5, 4, 1, 8, 7, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1 },
168 { 1, 11, 5, 7, 6, 0, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
169 { 6, 8, 7, 6, 0, 2, 11, 5, 4, 0, -1, -1, -1, -1, -1, -1, -1 },
170 { 7, 6, 2, 11, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
171 { 2, 6, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
172 { 0, 4, 8, 0, 2, 6, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
173 { 6, 10, 1, 0, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
174 { 4, 8, 6, 10, 1, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
175 { 2, 6, 10, 2, 4, 3, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
176 { 3, 5, 8, 0, 3, 2, 6, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1 },
177 { 6, 10, 1, 0, 6, 4, 3, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1 },
178 { 6, 10, 1, 3, 5, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
179 { 2, 6, 10, 2, 1, 11, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
180 { 8, 0, 4, 8, 2, 6, 10, 2, 3, 1, 11, 3, -1, -1, -1, -1, -1 },
181 { 6, 10, 11, 3, 0, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
182 { 6, 10, 11, 3, 4, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
183 { 5, 4, 1, 11, 5, 10, 2, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1 },
184 { 10, 2, 6, 10, 1, 11, 5, 8, 0, 1, -1, -1, -1, -1, -1, -1, -1 },
185 { 4, 0, 6, 10, 11, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
186 { 6, 10, 11, 5, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
187 { 10, 2, 8, 7, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
188 { 10, 2, 0, 4, 7, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
189 { 1, 0, 8, 7, 10, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
190 { 1, 4, 7, 10, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
191 { 10, 2, 8, 7, 10, 5, 4, 3, 5, -1, -1, -1, -1, -1, -1, -1, -1 },
192 { 7, 10, 2, 0, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
193 { 3, 5, 4, 3, 1, 0, 8, 7, 10, 1, -1, -1, -1, -1, -1, -1, -1 },
194 { 3, 5, 7, 10, 1, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
195 { 8, 7, 10, 2, 8, 1, 11, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1 },
196 { 3, 1, 11, 3, 0, 4, 7, 10, 2, 0, -1, -1, -1, -1, -1, -1, -1 },
197 { 11, 3, 0, 8, 7, 10, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
198 { 11, 3, 4, 7, 10, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
199 { 1, 11, 5, 4, 1, 10, 2, 8, 7, 10, -1, -1, -1, -1, -1, -1, -1 },
200 { 11, 5, 7, 10, 2, 0, 1, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
201 { 7, 10, 11, 5, 4, 0, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
202 { 7, 10, 11, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
203 { 5, 9, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
204 { 8, 0, 4, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
205 { 0, 2, 1, 0, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
206 { 2, 1, 4, 8, 2, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1 },
207 { 4, 3, 9, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
208 { 9, 7, 8, 0, 3, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
209 { 9, 7, 4, 3, 9, 1, 0, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1 },
210 { 1, 3, 9, 7, 8, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
211 { 3, 1, 11, 3, 5, 9, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
212 { 0, 4, 8, 0, 3, 1, 11, 3, 7, 5, 9, 7, -1, -1, -1, -1, -1 },
213 { 0, 2, 11, 3, 0, 5, 9, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1 },
214 { 7, 5, 9, 7, 4, 8, 2, 11, 3, 4, -1, -1, -1, -1, -1, -1, -1 },
215 { 1, 11, 9, 7, 4, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
216 { 9, 7, 8, 0, 1, 11, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
217 { 7, 4, 0, 2, 11, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
218 { 9, 7, 8, 2, 11, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
219 { 6, 8, 5, 9, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
220 { 0, 4, 5, 9, 6, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
221 { 5, 9, 6, 8, 5, 0, 2, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1 },
222 { 9, 6, 2, 1, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
223 { 6, 8, 4, 3, 9, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
224 { 0, 3, 9, 6, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
225 { 2, 1, 0, 2, 6, 8, 4, 3, 9, 6, -1, -1, -1, -1, -1, -1, -1 },
226 { 2, 1, 3, 9, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
227 { 6, 8, 5, 9, 6, 11, 3, 1, 11, -1, -1, -1, -1, -1, -1, -1, -1 },
228 { 1, 11, 3, 1, 0, 4, 5, 9, 6, 0, -1, -1, -1, -1, -1, -1, -1 },
229 { 9, 6, 8, 5, 9, 2, 11, 3, 0, 2, -1, -1, -1, -1, -1, -1, -1 },
230 { 9, 6, 2, 11, 3, 4, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
231 { 4, 1, 11, 9, 6, 8, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
232 { 1, 11, 9, 6, 0, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
233 { 2, 11, 9, 6, 8, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
234 { 2, 11, 9, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
235 { 10, 2, 6, 10, 9, 7, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
236 { 2, 6, 10, 2, 8, 0, 4, 8, 9, 7, 5, 9, -1, -1, -1, -1, -1 },
237 { 1, 0, 6, 10, 1, 9, 7, 5, 9, -1, -1, -1, -1, -1, -1, -1, -1 },
238 { 9, 7, 5, 9, 6, 10, 1, 4, 8, 6, -1, -1, -1, -1, -1, -1, -1 },
239 { 4, 3, 9, 7, 4, 6, 10, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1 },
240 { 2, 6, 10, 2, 8, 0, 3, 9, 7, 8, -1, -1, -1, -1, -1, -1, -1 },
241 { 7, 4, 3, 9, 7, 0, 6, 10, 1, 0, -1, -1, -1, -1, -1, -1, -1 },
242 { 10, 1, 3, 9, 7, 8, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
243 { 1, 11, 3, 1, 10, 2, 6, 10, 5, 9, 7, 5, -1, -1, -1, -1, -1 },
244 { 9, 7, 5, 9, 2, 6, 10, 2, 0, 4, 8, 0, 3, 1, 11, 3, -1 },
245 { 5, 9, 7, 5, 11, 3, 0, 6, 10, 11, -1, -1, -1, -1, -1, -1, -1 },
246 { 3, 4, 8, 6, 10, 11, 3, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1 },
247 { 2, 6, 10, 2, 1, 11, 9, 7, 4, 1, -1, -1, -1, -1, -1, -1, -1 },
248 { 7, 8, 0, 1, 11, 9, 7, 2, 6, 10, 2, -1, -1, -1, -1, -1, -1 },
249 { 7, 4, 0, 6, 10, 11, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
250 { 9, 7, 8, 6, 10, 11, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
251 { 5, 9, 10, 2, 8, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
252 { 5, 9, 10, 2, 0, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
253 { 0, 8, 5, 9, 10, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
254 { 5, 9, 10, 1, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
255 { 2, 8, 4, 3, 9, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
256 { 10, 2, 0, 3, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
257 { 3, 9, 10, 1, 0, 8, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
258 { 10, 1, 3, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
259 { 1, 11, 3, 1, 10, 2, 8, 5, 9, 10, -1, -1, -1, -1, -1, -1, -1 },
260 { 9, 10, 2, 0, 4, 5, 9, 1, 11, 3, 1, -1, -1, -1, -1, -1, -1 },
261 { 3, 0, 8, 5, 9, 10, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
262 { 5, 9, 10, 11, 3, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
263 { 2, 8, 4, 1, 11, 9, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
264 { 10, 2, 0, 1, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
265 { 0, 8, 4, 0, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
266 { 10, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
267 { 10, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
268 { 0, 4, 8, 0, 11, 10, 9, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
269 { 1, 0, 2, 1, 11, 10, 9, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
270 { 4, 8, 2, 1, 4, 11, 10, 9, 11, -1, -1, -1, -1, -1, -1, -1, -1 },
271 { 11, 10, 9, 11, 3, 5, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
272 { 8, 0, 3, 5, 8, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1, -1 },
273 { 2, 1, 0, 2, 11, 10, 9, 11, 4, 3, 5, 4, -1, -1, -1, -1, -1 },
274 { 10, 9, 11, 10, 2, 1, 3, 5, 8, 2, -1, -1, -1, -1, -1, -1, -1 },
275 { 3, 1, 10, 9, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
276 { 10, 9, 3, 1, 10, 0, 4, 8, 0, -1, -1, -1, -1, -1, -1, -1, -1 },
277 { 0, 2, 10, 9, 3, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
278 { 4, 8, 2, 10, 9, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
279 { 10, 9, 5, 4, 1, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
280 { 5, 8, 0, 1, 10, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
281 { 10, 9, 5, 4, 0, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
282 { 10, 9, 5, 8, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
283 { 6, 8, 7, 6, 10, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
284 { 0, 4, 7, 6, 0, 10, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1, -1 },
285 { 0, 2, 1, 0, 6, 8, 7, 6, 11, 10, 9, 11, -1, -1, -1, -1, -1 },
286 { 11, 10, 9, 11, 2, 1, 4, 7, 6, 2, -1, -1, -1, -1, -1, -1, -1 },
287 { 8, 7, 6, 8, 5, 4, 3, 5, 10, 9, 11, 10, -1, -1, -1, -1, -1 },
288 { 11, 10, 9, 11, 3, 5, 7, 6, 0, 3, -1, -1, -1, -1, -1, -1, -1 },
289 { 0, 2, 1, 0, 3, 5, 4, 3, 6, 8, 7, 6, 11, 10, 9, 11, -1 },
290 { 6, 2, 1, 3, 5, 7, 6, 11, 10, 9, 11, -1, -1, -1, -1, -1, -1 },
291 { 3, 1, 10, 9, 3, 7, 6, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1 },
292 { 0, 4, 7, 6, 0, 3, 1, 10, 9, 3, -1, -1, -1, -1, -1, -1, -1 },
293 { 8, 7, 6, 8, 0, 2, 10, 9, 3, 0, -1, -1, -1, -1, -1, -1, -1 },
294 { 9, 3, 4, 7, 6, 2, 10, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
295 { 6, 8, 7, 6, 10, 9, 5, 4, 1, 10, -1, -1, -1, -1, -1, -1, -1 },
296 { 6, 0, 1, 10, 9, 5, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
297 { 4, 0, 2, 10, 9, 5, 4, 6, 8, 7, 6, -1, -1, -1, -1, -1, -1 },
298 { 7, 6, 2, 10, 9, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
299 { 9, 11, 2, 6, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
300 { 9, 11, 2, 6, 9, 8, 0, 4, 8, -1, -1, -1, -1, -1, -1, -1, -1 },
301 { 9, 11, 1, 0, 6, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
302 { 11, 1, 4, 8, 6, 9, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
303 { 2, 6, 9, 11, 2, 3, 5, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1 },
304 { 2, 6, 9, 11, 2, 8, 0, 3, 5, 8, -1, -1, -1, -1, -1, -1, -1 },
305 { 4, 3, 5, 4, 1, 0, 6, 9, 11, 1, -1, -1, -1, -1, -1, -1, -1 },
306 { 5, 8, 6, 9, 11, 1, 3, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
307 { 3, 1, 2, 6, 9, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
308 { 8, 0, 4, 8, 2, 6, 9, 3, 1, 2, -1, -1, -1, -1, -1, -1, -1 },
309 { 0, 6, 9, 3, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
310 { 4, 8, 6, 9, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
311 { 2, 6, 9, 5, 4, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
312 { 6, 9, 5, 8, 0, 1, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
313 { 5, 4, 0, 6, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
314 { 5, 8, 6, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
315 { 8, 7, 9, 11, 2, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
316 { 0, 4, 7, 9, 11, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
317 { 8, 7, 9, 11, 1, 0, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
318 { 9, 11, 1, 4, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
319 { 4, 3, 5, 4, 8, 7, 9, 11, 2, 8, -1, -1, -1, -1, -1, -1, -1 },
320 { 11, 2, 0, 3, 5, 7, 9, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
321 { 11, 1, 0, 8, 7, 9, 11, 4, 3, 5, 4, -1, -1, -1, -1, -1, -1 },
322 { 9, 11, 1, 3, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
323 { 9, 3, 1, 2, 8, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
324 { 4, 7, 9, 3, 1, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
325 { 8, 7, 9, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
326 { 9, 3, 4, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
327 { 4, 1, 2, 8, 7, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
328 { 0, 1, 2, 0, 7, 9, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
329 { 5, 4, 0, 8, 7, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
330 { 5, 7, 9, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
331 { 11, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
332 { 11, 10, 7, 5, 11, 4, 8, 0, 4, -1, -1, -1, -1, -1, -1, -1, -1 },
333 { 7, 5, 11, 10, 7, 2, 1, 0, 2, -1, -1, -1, -1, -1, -1, -1, -1 },
334 { 5, 11, 10, 7, 5, 1, 4, 8, 2, 1, -1, -1, -1, -1, -1, -1, -1 },
335 { 4, 3, 11, 10, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
336 { 0, 3, 11, 10, 7, 8, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
337 { 0, 2, 1, 0, 4, 3, 11, 10, 7, 4, -1, -1, -1, -1, -1, -1, -1 },
338 { 10, 7, 8, 2, 1, 3, 11, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
339 { 7, 5, 3, 1, 10, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
340 { 0, 4, 8, 0, 3, 1, 10, 7, 5, 3, -1, -1, -1, -1, -1, -1, -1 },
341 { 2, 10, 7, 5, 3, 0, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
342 { 8, 2, 10, 7, 5, 3, 4, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
343 { 1, 10, 7, 4, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
344 { 8, 0, 1, 10, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
345 { 0, 2, 10, 7, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
346 { 8, 2, 10, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
347 { 11, 10, 6, 8, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
348 { 6, 0, 4, 5, 11, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
349 { 1, 0, 2, 1, 11, 10, 6, 8, 5, 11, -1, -1, -1, -1, -1, -1, -1 },
350 { 1, 4, 5, 11, 10, 6, 2, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
351 { 11, 10, 6, 8, 4, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
352 { 11, 10, 6, 0, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
353 { 10, 6, 8, 4, 3, 11, 10, 0, 2, 1, 0, -1, -1, -1, -1, -1, -1 },
354 { 11, 10, 6, 2, 1, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
355 { 1, 10, 6, 8, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
356 { 1, 10, 6, 0, 4, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
357 { 8, 5, 3, 0, 2, 10, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
358 { 2, 10, 6, 2, 4, 5, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
359 { 6, 8, 4, 1, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
360 { 1, 10, 6, 0, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
361 { 6, 8, 4, 0, 2, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
362 { 2, 10, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
363 { 2, 6, 7, 5, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
364 { 0, 4, 8, 0, 2, 6, 7, 5, 11, 2, -1, -1, -1, -1, -1, -1, -1 },
365 { 5, 11, 1, 0, 6, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
366 { 5, 11, 1, 4, 8, 6, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
367 { 11, 2, 6, 7, 4, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
368 { 0, 3, 11, 2, 6, 7, 8, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
369 { 0, 6, 7, 4, 3, 11, 1, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
370 { 1, 3, 11, 1, 6, 7, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
371 { 7, 5, 3, 1, 2, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
372 { 5, 3, 1, 2, 6, 7, 5, 0, 4, 8, 0, -1, -1, -1, -1, -1, -1 },
373 { 7, 5, 3, 0, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
374 { 7, 5, 3, 4, 8, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
375 { 2, 6, 7, 4, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
376 { 8, 0, 1, 2, 6, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
377 { 0, 6, 7, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
378 { 8, 6, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
379 { 2, 8, 5, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
380 { 0, 4, 5, 11, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
381 { 1, 0, 8, 5, 11, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
382 { 1, 4, 5, 11, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
383 { 4, 3, 11, 2, 8, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
384 { 0, 3, 11, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
385 { 4, 3, 11, 1, 0, 8, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
386 { 1, 3, 11, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
387 { 3, 1, 2, 8, 5, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
388 { 3, 1, 2, 0, 4, 5, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
389 { 3, 0, 8, 5, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
390 { 3, 4, 5, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
391 { 2, 8, 4, 1, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
392 { 0, 1, 2, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
393 { 0, 8, 4, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 },
394 { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 }
395 };
396
398 bool
399 modify_path (CubeFace const& dir, uint8_t level, uint64_t* path)
400 {
401 bool const subtract = static_cast<int>(dir) % 2;
402 int const bit_offset = static_cast<int>(dir) / 2;
403 bool overflow = true;
404 for (int i = 0; overflow && i < level; ++i)
405 {
406 uint64_t const mask = uint64_t(1) << (i * 3 + bit_offset);
407 if (*path & mask)
408 {
409 *path ^= mask;
410 overflow = !subtract;
411 }
412 else
413 {
414 *path |= mask;
415 overflow = subtract;
416 }
417 }
418
419 return !overflow;
420 }
421
422 Octree::Iterator
423 modify_iterator (CubeFace const& dir, Octree::Iterator const& iter)
424 {
425 uint64_t path = iter.path;
426 if (!modify_path(dir, iter.level, &path))
427 return Octree::Iterator();
428 return iter.descend(iter.level, path);
429 }
430
431 Octree::Iterator
432 modify_iterator (CubeFace const& dir1, CubeFace const& dir2,
433 Octree::Iterator const& iter)
434 {
435 uint64_t path = iter.path;
436 if (!modify_path(dir1, iter.level, &path))
437 return Octree::Iterator();
438 if (!modify_path(dir2, iter.level, &path))
439 return Octree::Iterator();
440 return iter.descend(iter.level, path);
441 }
442}
443
445IsoSurface::extract_mesh (void)
446{
447 std::cout << " Sanity-checking input data..." << std::flush;
448 util::WallTimer timer;
449 this->sanity_checks();
450 std::cout << " took " << timer.get_elapsed() << " ms." << std::endl;
451
452 /*
453 * Assign MC index to every octree node. This can be done in two ways:
454 * (1) Iterate all nodes, query corner values, and determine MC index.
455 * (2) Iterate leafs only, query corner values, propagate to parents.
456 * Strategy (1) is implemented, it is simpler but slightly more expensive.
457 */
458 std::cout << " Computing Marching Cubes indices..." << std::flush;
459 Octree::Iterator iter = this->octree->get_iterator_for_root();
460 for (iter.first_node(); iter.current != nullptr; iter.next_node())
461 this->compute_mc_index(iter);
462 std::cout << " took " << timer.get_elapsed() << " ms." << std::endl;
463
464 /*
465 * Compute isovertices on the octree edges for every leaf node.
466 * This locates for every leaf edge the finest unique edge which
467 * contains an isovertex. The vertex is stored in the vertex vector,
468 * while the edge is stored in the map, mapping edge to vertex ID.
469 */
470 std::cout << " Computing isovertices..." << std::flush;
471 timer.reset();
472 EdgeVertexMap edgemap;
473 IsoVertexVector isovertices;
474 for (iter.first_leaf(); iter.current != nullptr; iter.next_leaf())
475 this->compute_isovertices(iter, &edgemap, &isovertices);
476 std::cout << " took " << timer.get_elapsed() << " ms." << std::endl;
477
478 /*
479 * Compute polygons for every leaf node. For every leaf face the
480 * list of isoedges is retrieved. The isoedges are linked to form
481 * one or more closed polygons per node. In some cases open polygons
482 * are created, which need to be linked by retrieving twin vertices.
483 */
484 std::cout << " Computing isopolygons..." << std::flush;
485 timer.reset();
486 PolygonList polygons;
487 for (iter.first_leaf(); iter.current != nullptr; iter.next_leaf())
488 this->compute_isopolygons(iter, edgemap, &polygons);
489 std::cout << " took " << timer.get_elapsed() << " ms." << std::endl;
490
491 /*
492 * The vertices are transferred to a mesh and the polygons are
493 * triangulated using the minimum area triangulation.
494 */
495 std::cout << " Computing triangulation..." << std::flush;
496 timer.reset();
498 this->compute_triangulation(isovertices, polygons, mesh);
499 std::cout << " took " << timer.get_elapsed() << " ms." << std::endl;
500
501 return mesh;
502}
503
504void
505IsoSurface::sanity_checks (void)
506{
507 if (this->voxels == nullptr || this->octree == nullptr)
508 throw std::runtime_error("sanity_checks(): Null octree/voxels");
509 for (std::size_t i = 1; i < this->voxels->size(); ++i)
510 if (this->voxels->at(i).first < this->voxels->at(i - 1).first)
511 throw std::runtime_error("sanity_checks(): Voxels unsorted");
512
513 Octree::Iterator iter = this->octree->get_iterator_for_root();
514 for (iter.first_node(); iter.current != nullptr; iter.next_node())
515 {
516 /* Check if parent pointer of children is correct. */
517 if (iter.current->children == nullptr)
518 continue;
519 for (int i = 0; i < CUBE_CORNERS; ++i)
520 {
521 if (iter.current->children[i].parent != iter.current)
522 throw std::runtime_error("sanity_checks(): Wrong parent");
523 }
524 }
525}
526
527void
528IsoSurface::compute_mc_index (Octree::Iterator const& iter)
529{
530 iter.current->mc_index = 0;
531 for (int i = 0; i < CUBE_CORNERS; ++i)
532 {
533 VoxelIndex vi;
534 vi.from_path_and_corner(iter.level, iter.path, i);
535 VoxelData const* vd = this->get_voxel_data(vi);
536 if (vd->value < ISO_VALUE)
537 iter.current->mc_index |= (1 << i);
538 }
539}
540
541void
542IsoSurface::compute_isovertices (Octree::Iterator const& iter,
543 EdgeVertexMap* edgemap, IsoVertexVector* isovertices)
544{
545 /* This should always be a leaf node. */
546 if (iter.current == nullptr || iter.current->children != nullptr)
547 throw std::runtime_error("compute_isovertices(): Invalid node");
548
549 /* Check if cube contains an isosurface. */
550 if (iter.current->mc_index == 0x00 || iter.current->mc_index == 0xff)
551 return;
552
553 for (int i = 0; i < CUBE_EDGES; ++i)
554 {
555 /* Check if edge contains an isovertex. */
556 if (!this->is_isovertex_on_edge(iter.current->mc_index, i))
557 continue;
558
559 /* Get the finest edge that contains an isovertex. */
560 EdgeIndex edge_index;
561 this->get_finest_cube_edge(iter, i, &edge_index, nullptr);
562 if (edgemap->find(edge_index) != edgemap->end())
563 continue;
564
565 /* Interpolate isovertex and add index to map. */
566 IsoVertex isovertex;
567 this->get_isovertex(edge_index, i, &isovertex);
568 edgemap->insert(std::make_pair(edge_index, isovertices->size()));
569 isovertices->push_back(isovertex);
570 }
571}
572
573void
574IsoSurface::get_finest_cube_edge (Octree::Iterator const& iter,
575 int edge_id, EdgeIndex* edge_index, EdgeInfo* edge_info)
576{
577 if (edge_id < 0 || edge_id > 11)
578 throw std::runtime_error("get_finest_cube_edge(): Invalid edge ID");
579
580 if (!this->is_isovertex_on_edge(iter.current->mc_index, edge_id))
581 throw std::runtime_error("get_finest_cube_edge(): Invalid isoedge");
582
583 /* Remember the finest node found so far. */
584 Octree::Iterator finest_iter = iter;
585 int finest_edge_id = edge_id;
586 bool found_node = false;
587
588 /* Check if the current node has children. */
589 if (iter.current->children != nullptr)
590 found_node = true;
591
592 /* Check if the two face-adjacent nodes have children. */
593 for (int i = 0; !found_node && i < 2; ++i)
594 {
595 Octree::Iterator temp_iter;
596 temp_iter = modify_iterator(edge_neighbors[edge_id][i], iter);
597 if (temp_iter.current != nullptr && temp_iter.current->children != nullptr)
598 {
599 found_node = true;
600 finest_iter = temp_iter;
601 finest_edge_id = edge_reflections[edge_id][i];
602 }
603 }
604
605 /* Check if edge-adjacent node has children. */
606 if (found_node == false)
607 {
608 Octree::Iterator temp_iter = modify_iterator(
609 edge_neighbors[edge_id][0], edge_neighbors[edge_id][1], iter);
610 if (temp_iter.current != nullptr && temp_iter.current->children != nullptr)
611 {
612 found_node = true;
613 finest_iter = temp_iter;
614 finest_edge_id = edge_reflections[edge_id][2];
615 }
616 }
617
618 if (finest_iter.current == nullptr)
619 throw std::runtime_error("get_finest_cube_edge(): Error finding edge");
620
621 /* If the node has no children, we found the finest node. */
622 if (finest_iter.current->children == nullptr)
623 {
624 if (edge_index != nullptr)
625 {
626 VoxelIndex vi1, vi2;
627 vi1.from_path_and_corner(finest_iter.level, finest_iter.path,
628 edge_corners[finest_edge_id][0]);
629 vi2.from_path_and_corner(finest_iter.level, finest_iter.path,
630 edge_corners[finest_edge_id][1]);
631 edge_index->first = std::min(vi1.index, vi2.index);
632 edge_index->second = std::max(vi1.index, vi2.index);
633 }
634 if (edge_info != nullptr)
635 {
636 edge_info->iter = finest_iter;
637 edge_info->edge_id = finest_edge_id;
638 }
639
640 return;
641 }
642
643 /* Find unique child with isoedge and recurse. */
644 int const child_1_id = edge_corners[finest_edge_id][0];
645 int const child_2_id = edge_corners[finest_edge_id][1];
646 Octree::Node const& child1 = finest_iter.current->children[child_1_id];
647 Octree::Node const& child2 = finest_iter.current->children[child_2_id];
648 bool c1_iso = this->is_isovertex_on_edge(child1.mc_index, finest_edge_id);
649 bool c2_iso = this->is_isovertex_on_edge(child2.mc_index, finest_edge_id);
650 if (c1_iso && !c2_iso)
651 this->get_finest_cube_edge(finest_iter.descend(child_1_id),
652 finest_edge_id, edge_index, edge_info);
653 else if (c2_iso && !c1_iso)
654 this->get_finest_cube_edge(finest_iter.descend(child_2_id),
655 finest_edge_id, edge_index, edge_info);
656 else
657 throw std::runtime_error("get_finest_cube_edge(): Invalid parent edge");
658}
659
660void
661IsoSurface::get_isovertex (EdgeIndex const& edge_index,
662 int edge_id, IsoVertex* iso_vertex)
663{
664 /* Get voxel data. */
665 VoxelIndex vi1, vi2;
666 vi1.index = edge_index.first;
667 vi2.index = edge_index.second;
668
669#if FSSR_USE_DERIVATIVES
670
671 int const edge_axis = edge_id % 3;
672 if ((edge_axis == 0 && vi1.get_offset_x() > vi2.get_offset_x())
673 || (edge_axis == 1 && vi1.get_offset_y() > vi2.get_offset_y())
674 || (edge_axis == 2 && vi1.get_offset_z() > vi2.get_offset_z()))
675 std::swap(vi1, vi2);
676
677#endif // FSSR_USE_DERIVATIVES
678
679 VoxelData const* vd1 = this->get_voxel_data(vi1);
680 VoxelData const* vd2 = this->get_voxel_data(vi2);
681 /* Get voxel positions. */
682 math::Vec3d pos1 = vi1.compute_position(
683 this->octree->get_root_node_center(),
684 this->octree->get_root_node_size());
685 math::Vec3d pos2 = vi2.compute_position(
686 this->octree->get_root_node_center(),
687 this->octree->get_root_node_size());
688
689#if FSSR_USE_DERIVATIVES
690
691 /* Interpolate voxel data and position. */
692 double const norm = pos2[edge_axis] - pos1[edge_axis];
693 double const weight = interpolate_root(
694 vd1->value - ISO_VALUE, vd2->value - ISO_VALUE,
695 vd1->deriv[edge_axis] * norm, vd2->deriv[edge_axis] * norm,
696 this->interpolation_type);
697
698#else // FSSR_USE_DERIVATIVES
699
700 /* Interpolate voxel data and position. */
701 double const weight = (vd1->value - ISO_VALUE) / (vd1->value - vd2->value);
702
703#endif // FSSR_USE_DERIVATIVES
704
705 iso_vertex->data = interpolate_voxel(*vd1, (1.0 - weight), *vd2, weight);
706 iso_vertex->pos = pos1 * (1.0 - weight) + pos2 * weight;
707}
708
709bool
710IsoSurface::is_isovertex_on_edge (int mc_index, int edge_id)
711{
712 int const* bit = edge_corners[edge_id];
713 return ((mc_index >> bit[0]) & 1) ^ ((mc_index >> bit[1]) & 1);
714}
715
716void
717IsoSurface::compute_isopolygons (Octree::Iterator const& iter,
718 EdgeVertexMap const& edgemap, PolygonList* polygons)
719{
720 /*
721 * Step 1: Collect iso edges for all faces of this node.
722 */
723 IsoEdgeList isoedges;
724 for (int i = 0; i < CUBE_FACES; ++i)
725 this->get_finest_isoedges(iter, i, &isoedges, false);
726
727 /* Even cubes with MC index 0x0 or 0xff can have isoedges on the faces. */
728 if (isoedges.empty())
729 return;
730
731 /*
732 * Step 2: Find open vertices by computing vertex valences.
733 */
734 std::map<EdgeIndex, int> vertex_valence;
735 for (std::size_t i = 0; i < isoedges.size(); ++i)
736 {
737 std::map<EdgeIndex, int>::iterator valence_iter;
738 valence_iter = vertex_valence.find(isoedges[i].first);
739 if (valence_iter == vertex_valence.end())
740 vertex_valence[isoedges[i].first] = 1;
741 else
742 valence_iter->second += 1;
743 valence_iter = vertex_valence.find(isoedges[i].second);
744 if (valence_iter == vertex_valence.end())
745 vertex_valence[isoedges[i].second] = -1;
746 else
747 valence_iter->second -= 1;
748 }
749
750 /*
751 * Step 3: Close open polygons by connecting open twin vertices.
752 */
753 for (std::size_t i = 0; i < isoedges.size(); ++i)
754 {
755 IsoEdge const& isoedge = isoedges[i];
756 if (vertex_valence[isoedge.first] != 0)
757 {
758 EdgeIndex twin;
759 EdgeInfo twin_info;
760 this->find_twin_vertex(isoedge.first_info, &twin, &twin_info);
761
762 IsoEdge new_edge;
763 new_edge.first = twin;
764 new_edge.first_info = twin_info;
765 new_edge.second = isoedge.first;
766 new_edge.second_info = isoedge.first_info;
767 isoedges.push_back(new_edge);
768
769 vertex_valence[new_edge.first] += 1;
770 vertex_valence[new_edge.second] -= 1;
771 }
772
773 if (vertex_valence[isoedge.second] != 0)
774 {
775 EdgeIndex twin;
776 EdgeInfo twin_info;
777 this->find_twin_vertex(isoedge.second_info, &twin, &twin_info);
778
779 IsoEdge new_edge;
780 new_edge.first = isoedges[i].second;
781 new_edge.first_info = isoedges[i].second_info;
782 new_edge.second = twin;
783 new_edge.second_info = twin_info;
784 isoedges.push_back(new_edge);
785
786 vertex_valence[new_edge.first] += 1;
787 vertex_valence[new_edge.second] -= 1;
788 }
789 }
790
791 /*
792 * Step 4: Join edges to form closed polygons.
793 */
794 std::size_t poly_start = 0;
795 for (std::size_t i = 0; i < isoedges.size(); ++i)
796 {
797 /* Once joined edges close, issue a new polygon. */
798 if (isoedges[i].second == isoedges[poly_start].first)
799 {
800 polygons->push_back(std::vector<std::size_t>());
801 std::vector<std::size_t>& poly = polygons->back();
802 for (std::size_t j = poly_start; j <= i; ++j)
803 {
804 poly.push_back(this->lookup_edge_vertex(edgemap,
805 isoedges[j].first));
806 }
807 poly_start = i + 1;
808 continue;
809 }
810
811 /* Fix successive edge. */
812 bool found_edge = false;
813 for (std::size_t j = i + 1; !found_edge && j < isoedges.size(); ++j)
814 {
815 if (isoedges[i].second == isoedges[j].first)
816 {
817 std::swap(isoedges[i + 1], isoedges[j]);
818 found_edge = true;
819 }
820 }
821
822 if (!found_edge)
823 throw std::runtime_error("Cannot find next edge");
824 }
825}
826
827std::size_t
828IsoSurface::lookup_edge_vertex (EdgeVertexMap const& edgemap,
829 EdgeIndex const& edge)
830{
831 EdgeVertexMap::const_iterator iter = edgemap.find(edge);
832 if (iter == edgemap.end())
833 throw std::runtime_error("lookup_edge_vertex(): No such edge vertex");
834 return iter->second;
835}
836
837void
838IsoSurface::find_twin_vertex (EdgeInfo const& edge_info,
839 EdgeIndex* twin, EdgeInfo* twin_info)
840{
841 /*
842 * Find the twin isovertex for the given edge. Go upwards in the tree
843 * through the parent edges until an edge with no iso crossing is found.
844 * Get the twin isovertex on the second child adjacent to the edge.
845 */
846 Octree::Iterator iter = edge_info.iter;
847 int const edge_id = edge_info.edge_id;
848 while (iter.current->parent != nullptr)
849 {
850 int const node_octant
851 = static_cast<int>(iter.current - iter.current->parent->children);
852
853 /* The octant of this node in the parent must be on the same edge. */
854 int descend_octant;
855 if (edge_corners[edge_id][0] == node_octant)
856 descend_octant = edge_corners[edge_id][1];
857 else if (edge_corners[edge_id][1] == node_octant)
858 descend_octant = edge_corners[edge_id][0];
859 else
860 throw std::runtime_error("find_twin_vertex(): Invalid parent edge");
861
862 /* If the parent edge has no isocrossing, descend to find twin. */
863 iter = iter.ascend();
864 if (!this->is_isovertex_on_edge(iter.current->mc_index, edge_id))
865 {
866 iter = iter.descend(descend_octant);
867 this->get_finest_cube_edge(iter, edge_id, twin, twin_info);
868 return;
869 }
870 }
871
872 throw std::runtime_error("find_twin_vertex(): Reached octree root");
873
874}
875
876void
877IsoSurface::get_finest_isoedges (Octree::Iterator const& iter,
878 int face_id, IsoEdgeList* isoedges, bool descend_only)
879{
880 /* If descend only is set, face-neighboring nodes are not considered. */
881 if (descend_only)
882 {
883 if (iter.current->children != nullptr)
884 {
885 /* Recursively descend to obtain iso edges for this face. */
886 this->get_finest_isoedges(iter.descend(face_corners[face_id][0]),
887 face_id, isoedges, true);
888 this->get_finest_isoedges(iter.descend(face_corners[face_id][1]),
889 face_id, isoedges, true);
890 this->get_finest_isoedges(iter.descend(face_corners[face_id][2]),
891 face_id, isoedges, true);
892 this->get_finest_isoedges(iter.descend(face_corners[face_id][3]),
893 face_id, isoedges, true);
894 }
895 else
896 {
897 /* Create list of isoedges for this face. */
898 int const mc_index = iter.current->mc_index;
899 int const* edge_table = mc_polygons[mc_index];
900
901 int first_id = -1;
902 for (int i = 0; edge_table[i] != -1; ++i)
903 {
904 if (first_id == -1)
905 {
906 first_id = edge_table[i];
907 continue;
908 }
909 if (edge_table[i] == first_id)
910 first_id = -1;
911
912 if (!is_edge_on_face[edge_table[i - 1]][face_id]
913 || !is_edge_on_face[edge_table[i]][face_id])
914 continue;
915
916 IsoEdge isoedge;
917 this->get_finest_cube_edge(iter, edge_table[i - 1],
918 &isoedge.first, &isoedge.first_info);
919 this->get_finest_cube_edge(iter, edge_table[i],
920 &isoedge.second, &isoedge.second_info);
921 isoedges->push_back(isoedge);
922 }
923 }
924
925 return;
926 }
927
928 /* Check if face-neighboring node has finer subdivision. */
929 Octree::Iterator niter = modify_iterator((CubeFace)face_id, iter);
930 if (niter.current != nullptr && niter.current->children != nullptr)
931 {
932 /* Face-neighboring node has finer subdivision. */
933 CubeFace const opposite_face_id = face_opposite[face_id];
934 std::size_t const last_isoedge_index = isoedges->size();
935 this->get_finest_isoedges(niter, opposite_face_id, isoedges, true);
936
937 /* Flip orientation for face-neighboring iso-edges. */
938 for (std::size_t i = last_isoedge_index; i < isoedges->size(); ++i)
939 {
940 IsoEdge& isoedge = isoedges->at(i);
941 std::swap(isoedge.first, isoedge.second);
942 std::swap(isoedge.first_info, isoedge.second_info);
943 }
944 }
945 else
946 {
947 /* Find the isoedges for this node face. */
948 this->get_finest_isoedges(iter, face_id, isoedges, true);
949 }
950}
951
952void
953IsoSurface::compute_triangulation(IsoVertexVector const& isovertices,
954 PolygonList const& polygons, mve::TriangleMesh::Ptr mesh)
955{
956 /* Populate per-vertex attributes. */
957 mve::TriangleMesh::VertexList& verts = mesh->get_vertices();
958 mve::TriangleMesh::ColorList& colors = mesh->get_vertex_colors();
959 mve::TriangleMesh::ValueList& values = mesh->get_vertex_values();
960 mve::TriangleMesh::ConfidenceList& cfs = mesh->get_vertex_confidences();
961 verts.reserve(isovertices.size());
962 colors.reserve(isovertices.size());
963 values.reserve(isovertices.size());
964 cfs.reserve(isovertices.size());
965
966 for (std::size_t i = 0; i < isovertices.size(); ++i)
967 {
968 IsoVertex const& vertex = isovertices[i];
969 verts.push_back(vertex.pos);
970 colors.push_back(math::Vec4f(vertex.data.color, 1.0f));
971 values.push_back(vertex.data.scale);
972 cfs.push_back(vertex.data.conf);
973 }
974
975 /* Triangulate isopolygons. */
976 mve::TriangleMesh::FaceList& triangles = mesh->get_faces();
978 for (std::size_t i = 0; i < polygons.size(); i++)
979 {
980 std::vector<math::Vector<float, 3> > loop;
981 loop.resize(polygons[i].size());
982 for (std::size_t j = 0; j < polygons[i].size(); ++j)
983 loop[j] = verts[polygons[i][j]];
984 std::vector<unsigned int> result;
985 tri.triangulate(loop, &result);
986 for (std::size_t j = 0; j < result.size(); ++j)
987 triangles.push_back(polygons[i][result[j]]);
988 }
989}
990
Computes the minimum area triangulation of a polygon.
void triangulate(std::vector< math::Vec3f > const &verts, std::vector< unsigned int > *indices)
std::vector< math::Vec3f > VertexList
Definition mesh.h:33
std::vector< float > ConfidenceList
Definition mesh.h:35
std::vector< math::Vec4f > ColorList
Definition mesh.h:34
std::vector< float > ValueList
Definition mesh.h:36
std::shared_ptr< TriangleMesh > Ptr
Definition mesh.h:92
std::vector< VertexID > FaceList
Definition mesh.h:97
static Ptr create(void)
Definition mesh.h:295
Cross-platform high-resolution real-time timer.
Definition timer.h:30
std::size_t get_elapsed(void) const
Returns the milli seconds since last reset.
Definition timer.h:94
void reset(void)
Definition timer.h:88
#define FSSR_NAMESPACE_END
Definition defines.h:14
#define FSSR_NAMESPACE_BEGIN
Definition defines.h:13
#define CUBE_EDGES
#define CUBE_FACES
#define CUBE_CORNERS
#define ISO_VALUE
std::size_t second
Definition mesh_info.cc:47
std::size_t first
Definition mesh_info.cc:46
std::size_t face_id
Definition mesh_info.cc:45
VoxelData interpolate_voxel(VoxelData const &d1, float w1, VoxelData const &d2, float w2)
Interpolates between two VoxelData objects for Marching Cubes.
Definition voxel.cc:74
double interpolate_root(double v0, double v1, double d0, double d1, InterpolationType type)
Interpolates the root of an unknown function f(x) given value and derivative constraints: f(0) = v0,...
Definition hermite.cc:145
void swap(mve::Image< T > &a, mve::Image< T > &b)
Specialization of std::swap for efficient image swapping.
Definition image.h:478
Octree iterator that keeps track of level and path through the octree.
Definition octree.h:58
Node * first_leaf(void)
Definition octree.cc:31
Node * next_leaf(void)
Definition octree.cc:78
Node * first_node(void)
Definition octree.cc:22
Node * next_node(void)
Definition octree.cc:44
Node * children
Definition octree.h:45