NDEVR
API Documentation
CubeTriangulation.h
1#pragma once
2#include "Base/Headers/BaseValues.hpp"
3#include "Base/Headers/Vector.hpp"
4#include "BlockModel/Headers/CubeLookupTable.h"
5#define MC33C_VERSION_MAJOR 5
6#define MC33C_VERSION_MINOR 2
7#define GRD_orthogonal
8#if defined(integer_GRD)
9#if size_type_GRD == 4
10/*
11GRD_data_type is the variable type of the grid data, by default it is float.
12*/
13typedef uint04 GRD_data_type;
14#elif size_type_GRD == 2
15typedef unsigned short int GRD_data_type;
16#elif size_type_GRD == 1
17typedef unsigned char GRD_data_type;
18#else
19#error "Incorrect size of the data type. size_type_GRD permitted values: 1, 2 or 4."
20typedef float GRD_data_type;
21#endif
22#elif size_type_GRD == 8
23typedef double GRD_data_type;
24#else
25#undef size_type_GRD
26#define size_type_GRD 4
27#endif
28#if defined (__SSE__) || ((defined (_M_IX86) || defined (_M_X64)) && !defined(_CHPE_ONLY_))
29// https://stackoverflow.com/questions/59644197/inverse-square-root-intrinsics
30// faster than 1.0f/std::sqrt, but with little accuracy.
31#include <immintrin.h>
32inline float invSqrt(float f) {
33 __m128 temp = _mm_set_ss(f);
34 temp = _mm_rsqrt_ss(temp);
35 return _mm_cvtss_f32(temp);
36}
37#else
38#include <math.h>
39inline float invSqrt(float f) {
40 return 1.0 / sqrt(f);
41}
42#endif
43namespace NDEVR
44{
45 typedef float GRD_data_type;
49 template<class t_point_type>
50 struct _GRD
51 {
52 _GRD(unsigned x, unsigned y, unsigned z)
53 {
54 N[0] = x;
55 N[1] = y;
56 N[2] = z;
57 uint04 j, k;
58 values = (GRD_data_type***)malloc((N[2] + 1) * sizeof(void*));
59 smoothed_values = (GRD_data_type***)malloc((N[2] + 1) * sizeof(void*));
60 nodes = (t_point_type***)malloc((N[2] + 1) * sizeof(void*));
61 for (k = 0; k <= N[2]; ++k)
62 {
63 values[k] = (GRD_data_type**)malloc((N[1] + 1) * sizeof(void*));
64 smoothed_values[k] = (GRD_data_type**)malloc((N[1] + 1) * sizeof(void*));
65 nodes[k] = (t_point_type**)malloc((N[1] + 1) * sizeof(void*));
66 for (j = 0; j <= N[1]; ++j)
67 {
68 values[k][j] = (GRD_data_type*)malloc((N[0] + 1) * sizeof(GRD_data_type));
69 nodes[k][j] = (t_point_type*)malloc((N[0] + 1) * sizeof(t_point_type));
70 smoothed_values[k][j] = (GRD_data_type*)malloc((N[0] + 1) * sizeof(GRD_data_type));
71 }
72 }
73 }
74 ~_GRD()
75 {
76 for (uint04 k = 0; k <= N[2]; ++k)
77 {
78 for (uint04 j = 0; j <= N[1]; ++j)
79 {
80 free(values[k][j]);
81 free(nodes[k][j]);
82 free(smoothed_values[k][j]);
83 }
84 free(values[k]);
85 free(nodes[k]);
86 free(smoothed_values[k]);
87 }
88 }
89 void smooth()
90 {
91 for (uint04 z = 0; z <= N[2]; ++z)
92 {
93 for (uint04 y = 0; y <= N[1]; ++y)
94 {
95 for (uint04 x = 0; x <= N[0]; ++x)
96 {
97 smoothed_values[z][y][x] = values[z][y][x] * (1.0f / 14.0f);
98 }
99 }
100 }
101 for (uint04 z = 1; z < N[2]; ++z)
102 {
103 for (uint04 y = 1; y < N[1]; ++y)
104 {
105 for (uint04 x = 1; x < N[0]; ++x)
106 {
107 values[z][y][x] *= (1.0f / 14.0f);
108 values[z][y][x] += smoothed_values[z + 0][y + 0][x + 1];
109 values[z][y][x] += smoothed_values[z + 0][y + 1][x + 0];
110 values[z][y][x] += smoothed_values[z + 0][y + 1][x + 1];
111 values[z][y][x] += smoothed_values[z + 1][y + 0][x + 0];
112 values[z][y][x] += smoothed_values[z + 1][y + 0][x + 1];
113 values[z][y][x] += smoothed_values[z + 1][y + 1][x + 0];
114 values[z][y][x] += smoothed_values[z + 1][y + 1][x + 1];
115 values[z][y][x] += smoothed_values[z - 0][y - 0][x - 1];
116 values[z][y][x] += smoothed_values[z - 0][y - 1][x - 0];
117 values[z][y][x] += smoothed_values[z - 0][y - 1][x - 1];
118 values[z][y][x] += smoothed_values[z - 1][y - 0][x - 0];
119 values[z][y][x] += smoothed_values[z - 1][y - 0][x - 1];
120 values[z][y][x] += smoothed_values[z - 1][y - 1][x - 0];
121 values[z][y][x] += smoothed_values[z - 1][y - 1][x - 1];
122 }
123 }
124 }
125 }
126 GRD_data_type*** values;
127 GRD_data_type*** smoothed_values;
128 t_point_type*** nodes;
129 uint04 N[3];
130 };
131 template<class t_point_type>
132 struct MC33
133 {
134 MC33(const _GRD<t_point_type>& G)
135 {
136 uint04 x, y;
137 nx = G.N[0];
138 ny = G.N[1];
139 nz = G.N[2];
140 Lz = (uint04**)malloc((ny + 1) * sizeof(int*));//edges 1, 3, 5 (only write) and 7
141 Dy = (uint04**)malloc(ny * sizeof(int*));//edges 0 (only read) and 4
142 Uy = (uint04**)malloc(ny * sizeof(int*));//edges 2 and 6 (only write)
143 Dx = (uint04**)malloc((ny + 1) * sizeof(int*));//edges 8 and 9
144 Ux = (uint04**)malloc((ny + 1) * sizeof(int*));//edges 10 (only write) and 11
145 values = (const GRD_data_type***)G.values;
146 nodes = (const t_point_type***)G.nodes;
147 x = nx * sizeof(int);
148 for (y = 0; y != ny; ++y)
149 {
150 Dx[y] = (uint04*)malloc(x);
151 Ux[y] = (uint04*)malloc(x);
152 Lz[y] = (uint04*)malloc(x + sizeof(int));
153 Dy[y] = (uint04*)malloc(x + sizeof(int));
154 Uy[y] = (uint04*)malloc(x + sizeof(int));
155 }
156 if (Uy[y - 1])
157 {
158 Dx[y] = (uint04*)malloc(x);
159 Ux[y] = (uint04*)malloc(x);
160 Lz[y] = (uint04*)malloc(x + sizeof(int));
161 }
162 }
163 uint04 store(float* r, const t_point_type& color)
164 {
165 /*fltp04 t = invSqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
166 if (t_point_type::HasNormal())
167 {
168 Vector<3, fltp04> norm;
169 for (uint01 i = 0; i < 3; ++i)
170 norm[2 - i] = t * r[3 + i];
171 if (dot(norm, color.normal()) < 0)
172 return Constant<uint04>::Invalid;
173 }*/
174 uint04 nv = nV++;
175 if (nv == capv)
176 {
177 N = (float(*)[3])realloc(N, nv * 6 * sizeof(float)); // the memory space is duplicated
178 V = (float(*)[3])realloc(V, nv * 6 * sizeof(float)); // the memory space is duplicated
179 C = (t_point_type*)realloc(C, nv * 2 * sizeof(t_point_type)); // the memory space is duplicated
180 }
181 for (uint01 i = 0; i < 3; ++i)
182 V[nv][i] = r[i];
183 C[nv] = color;
184 return nv;
185 }
186 ~MC33()
187 {
188 for (uint04 y = 0; y != ny; ++y)
189 {
190 free(Dx[y]);
191 free(Ux[y]);
192 free(Dy[y]);
193 free(Uy[y]);
194 free(Lz[y]);
195 }
196 free(Dx[ny]);
197 free(Ux[ny]);
198 free(Lz[ny]);
199 free(Dx);
200 free(Ux);
201 free(Dy);
202 free(Uy);
203 free(Lz);
204 free(V);
205 }
206 uint04 surfint(uint04 x, uint04 y, uint04 z, fltp04* r)
207 {
208 t_point_type pt = nodes[z][y][x];
209 r[0] = cast<fltp04>(x);
210 r[1] = cast<fltp04>(y);
211 r[2] = cast<fltp04>(z);
212 if (x == 0)
213 r[3] = values[z][y][0] - values[z][y][1];
214 else if (x == nx)
215 r[3] = values[z][y][x - 1] - values[z][y][x];
216 else
217 r[3] = 0.5f * (values[z][y][x - 1] - values[z][y][x + 1]);
218 if (y == 0)
219 r[4] = values[z][0][x] - values[z][1][x];
220 else if (y == ny)
221 r[4] = values[z][y - 1][x] - values[z][y][x];
222 else
223 r[4] = 0.5f * (values[z][y - 1][x] - values[z][y + 1][x]);
224 if (z == 0)
225 r[5] = values[0][y][x] - values[1][y][x];
226 else if (z == nz)
227 r[5] = values[z - 1][y][x] - values[z][y][x];
228 else
229 r[5] = 0.5f * (values[z - 1][y][x] - values[z + 1][y][x]);
230 return store(r, pt);
231 }
232 // copy of some variables of the surface structure:
233 uint04(*T)[3] = nullptr;
234 float(*V)[3] = nullptr;
235 float(*N)[3] = nullptr;
236 t_point_type(*C) = nullptr;
237 uint04 nV = 0U;
238 uint04 nT = 0U;
239 float iso;
240 uint04 capt = 0U;
241 uint04 capv = 0U;
242 // copy of some variables of the _GRD structure:
243 const GRD_data_type*** values;
244 const t_point_type*** nodes;
245 float ca, cb;
246 //int cubicflag;
247 uint04 nx, ny, nz;
248 // temporary structures that store the indices of triangle vertices:
249 uint04** Dx, ** Dy, ** Ux, ** Uy, ** Lz;
250 };
251 /***************** Marching cubes 33 functions ****************/
252 /******************************************************************
253 Vertices: Faces:
254 3 __________2 ___________
255 /| /| /| /|
256 / | / | / | 2 / |
257 7/__________6/ | / | 4 / |
258 | | | | |¯¯¯¯¯¯¯¯¯¯¯| 1 | z
259 | 0_______|___1 | 3 |_______|___| |
260 | / | / | / 5 | / |____y
261 | / | / | / 0 | / /
262 4/__________5/ |/__________|/ x
263 This function return a vector with all six test face results (face[6]). Each
264 result value is 1 if the positive face vertices are joined, -1 if the negative
265 vertices are joined, and 0 (unchanged) if the test must no be applied. The
266 return value of this function is the the sum of all six results.*/
267 static inline int MC33_faceTests(int* face, int ind, const float* v)
268 {
269 if (ind & 0x80)//vertex 0
270 {
271 face[0] = ((ind & 0xCC) == 0x84 ? (v[0] * v[5] < v[1] * v[4] ? -1 : 1) : 0);//0x84 = 10000100, vertices 0 and 5
272 face[3] = ((ind & 0x99) == 0x81 ? (v[0] * v[7] < v[3] * v[4] ? -1 : 1) : 0);//0x81 = 10000001, vertices 0 and 7
273 face[4] = ((ind & 0xF0) == 0xA0 ? (v[0] * v[2] < v[1] * v[3] ? -1 : 1) : 0);//0xA0 = 10100000, vertices 0 and 2
274 }
275 else
276 {
277 face[0] = ((ind & 0xCC) == 0x48 ? (v[0] * v[5] < v[1] * v[4] ? 1 : -1) : 0);//0x48 = 01001000, vertices 1 and 4
278 face[3] = ((ind & 0x99) == 0x18 ? (v[0] * v[7] < v[3] * v[4] ? 1 : -1) : 0);//0x18 = 00011000, vertices 3 and 4
279 face[4] = ((ind & 0xF0) == 0x50 ? (v[0] * v[2] < v[1] * v[3] ? 1 : -1) : 0);//0x50 = 01010000, vertices 1 and 3
280 }
281 if (ind & 0x02)//vertex 6
282 {
283 face[1] = ((ind & 0x66) == 0x42 ? (v[1] * v[6] < v[2] * v[5] ? -1 : 1) : 0);//0x42 = 01000010, vertices 1 and 6
284 face[2] = ((ind & 0x33) == 0x12 ? (v[3] * v[6] < v[2] * v[7] ? -1 : 1) : 0);//0x12 = 00010010, vertices 3 and 6
285 face[5] = ((ind & 0x0F) == 0x0A ? (v[4] * v[6] < v[5] * v[7] ? -1 : 1) : 0);//0x0A = 00001010, vertices 4 and 6
286 }
287 else
288 {
289 face[1] = ((ind & 0x66) == 0x24 ? (v[1] * v[6] < v[2] * v[5] ? 1 : -1) : 0);//0x24 = 00100100, vertices 2 and 5
290 face[2] = ((ind & 0x33) == 0x21 ? (v[3] * v[6] < v[2] * v[7] ? 1 : -1) : 0);//0x21 = 00100001, vertices 2 and 7
291 face[5] = ((ind & 0x0F) == 0x05 ? (v[4] * v[6] < v[5] * v[7] ? 1 : -1) : 0);//0x05 = 00000101, vertices 5 and 7
292 }
293 return face[0] + face[1] + face[2] + face[3] + face[4] + face[5];
294 }
295 /* Faster function for the face test, the test is applied to only one face
296 (int face). This function is only used for the cases 3 and 6 of MC33*/
297 static inline int MC33_faceTest1(int face, const float* v)
298 {
299 switch (face) {
300 case 0: return (v[0] * v[5] < v[1] * v[4] ? 0x48 : 0x84);
301 case 1: return (v[1] * v[6] < v[2] * v[5] ? 0x24 : 0x42);
302 case 2: return (v[3] * v[6] < v[2] * v[7] ? 0x21 : 0x12);
303 case 3: return (v[0] * v[7] < v[3] * v[4] ? 0x18 : 0x81);
304 case 4: return (v[0] * v[2] < v[1] * v[3] ? 0x50 : 0xA0);
305 default: return (v[4] * v[6] < v[5] * v[7] ? 0x05 : 0x0A);
306 }
307 }
308 // an ugly signbit for float type with
309 // warning: dereferencing type-punned pointer will break strict-aliasing rules
310 // Silence dereferencing type-punned pointer warning in GCC
311#ifdef __GNUC__
312#pragma GCC diagnostic push
313#pragma GCC diagnostic ignored "-Wstrict-aliasing"
314#endif
315 inline uint04 signbf(float x) {
316 return ((*(uint04*)(void*)(&x)) & 0x80000000);
317 }
318#ifdef __GNUC__
319#pragma GCC diagnostic pop
320#endif
321 /******************************************************************
322 Interior test function. If the test is positive, the function returns a value
323 different from 0. The integer i must be 0 to test if the vertices 0 and 6 are
324 joined. 1 for vertices 1 and 7, 2 for vertices 2 and 4, and 3 for 3 and 5.
325 For case 13, the integer flag13 must be 1, and the function returns 2 if one
326 of the vertices 0, 1, 2 or 3 is joined to the center point of the cube (case
327 13.5.2), returns 1 if one of the vertices 4, 5, 6 or 7 is joined to the
328 center point of the cube (case 13.5.2 too), and it returns 0 if the vertices
329 are no joined (case 13.5.1)*/
330 static inline int MC33_interiorTest(int i, int flag13, const float* v)
331 {
332 //Signs of cube vertices were changed to use signbit function in calc_isosurface
333 //A0 = -v[0], B0 = -v[1], C0 = -v[2], D0 = -v[3]
334 //A1 = -v[4], B1 = -v[5], C1 = -v[6], D1 = -v[7]
335 //But the function still works
336 float At = v[4] - v[0], Bt = v[5] - v[1], Ct = v[6] - v[2], Dt = v[7] - v[3];
337 float t = At * Ct - Bt * Dt; // the "a" value.
338 if (signbf(t)) {
339 if (i & 0x01)
340 return 0;
341 }
342 else {
343 if (!(i & 0x01) || t == 0)
344 return 0;
345 }
346 t = 0.5f * (v[3] * Bt - v[2] * At + v[1] * Dt - v[0] * Ct) / t; // t = -b/2a
347 if (t > 0 && t < 1) {
348 At = v[0] + At * t;
349 Bt = v[1] + Bt * t;
350 Ct = v[2] + Ct * t;
351 Dt = v[3] + Dt * t;
352 Ct *= At;
353 Dt *= Bt;
354 if (i & 0x01) {
355 if (Ct < Dt && signbf(Dt) == 0)
356 return cast<int>(signbf(Bt) == signbf(v[i])) + flag13;
357 }
358 else {
359 if (Ct > Dt && signbf(Ct) == 0)
360 return cast<int>(signbf(At) == signbf(v[i])) + flag13;
361 }
362 }
363 return 0;
364 }
365 /******************************************************************
366 This function find the MC33 case (using the index i, and the face and interior
367 tests). The correct triangle pattern is obtained from the arrays contained in
368 the MC33_LookUpTable.h file. The necessary vertices (intersection points) are
369 also calculated here (using trilinear interpolation).
370 _____2_____
371 /| /|
372 11 |<-3 10 |
373 /____6_____ / 1 z
374 | | | | |
375 | |_____0_|___| |____y
376 7 / 5 / /
377 | 8 | 9 x
378 |/____4_____|/
379 The temporary matrices: M->Lz, M->Dx, M->Dy, M->Ux and M->Uy are filled
380 and used here.*/
381 template<class t_point_type>
382 void MC33_findCase(MC33<t_point_type>* M, uint04 x, uint04 y, uint04 z, uint04 i, const float* v)
383 {
384 Vector<13, uint04> p(Constant<uint04>::Invalid);
385 uint04 ti[3];//for vertex indices of a triangle
386 union { // memory saving
387 int f[6];//for the face tests
388 fltp04 r[6];//for intercept and normal coordinates
389 } u;
390 const unsigned short int* pcase = MC33_all_tables;
391 float t;
392 uint04 c, m, k;
393 if (i & 0x80) {
394 c = pcase[i ^ 0xFF];
395 m = (c & 0x800) == 0;
396 }
397 else {
398 c = pcase[i];
399 m = (c & 0x800) != 0;
400 }
401 k = c & 0x7FF;
402 switch (c >> 12) { //find the MC33 case
403 case 0: // cases 1, 2, 5, 8, 9, 11 and 14
404 pcase += k;
405 break;
406 case 1: // case 3
407 pcase += ((m ? i : i ^ 0xFF) & MC33_faceTest1(k >> 2, v) ? 183 + (k << 1) : 159 + k);
408 break;
409 case 2: // case 4
410 pcase += (MC33_interiorTest(k, 0, v) ? 239 + 6 * k : 231 + (k << 1));
411 break;
412 case 3: // case 6
413 if ((m ? i : i ^ 0xFF) & MC33_faceTest1(k % 6, v))
414 pcase += 575 + 5 * k; //6.2
415 else
416 pcase += (MC33_interiorTest(k / 6, 0, v) ? 407 + 7 * k : 335 + 3 * k); //6.1
417 break;
418 case 4: // case 7
419 switch (MC33_faceTests(u.f, (m ? i : i ^ 0xFF), v))
420 {
421 case -3:
422 pcase += 695 + 3 * k; //7.1
423 break;
424 case -1: //7.2
425 pcase += (u.f[4] + u.f[5] < 0 ? (u.f[0] + u.f[2] < 0 ? 759 : 799) : 719) + 5 * k;
426 break;
427 case 1: //7.3
428 pcase += (u.f[4] + u.f[5] < 0 ? 983 : (u.f[0] + u.f[2] < 0 ? 839 : 911)) + 9 * k;
429 break;
430 default: //7.4
431 pcase += (MC33_interiorTest(k >> 1, 0, v) ? 1095 + 9 * k : 1055 + 5 * k);
432 }
433 break;
434 case 5: // case 10
435 switch (MC33_faceTests(u.f, (m ? i : i ^ 0xFF), v)) {
436 case -2:
437 if (k == 2 ? MC33_interiorTest(0, 0, v) : MC33_interiorTest(0, 0, v) != 0 || MC33_interiorTest(k ? 1 : 3, 0, v) != 0)
438 pcase += 1213 + (k << 3); //10.1.2
439 else
440 pcase += 1189 + (k << 2); //10.1.1
441 break;
442 case 0: //10.2
443 pcase += (u.f[2 + k] < 0 ? 1261 : 1285) + (k << 3);
444 break;
445 default:
446 if (k == 2 ? MC33_interiorTest(1, 0, v) : MC33_interiorTest(2, 0, v) != 0 || MC33_interiorTest(k ? 3 : 1, 0, v) != 0)
447 pcase += 1237 + (k << 3); //10.1.2
448 else
449 pcase += 1201 + (k << 2); //10.1.1
450 }
451 break;
452 case 6: // case 12
453 switch (MC33_faceTests(u.f, (m ? i : i ^ 0xFF), v)) {
454 case -2: //12.1
455 pcase += (MC33_interiorTest((0xDA010C >> (k << 1)) & 3, 0, v) ? 1453 + (k << 3) : 1357 + (k << 2));
456 break;
457 case 0: //12.2
458 pcase += (u.f[k >> 1] < 0 ? 1645 : 1741) + (k << 3);
459 break;
460 default: //12.1
461 pcase += (MC33_interiorTest((0xA7B7E5 >> (k << 1)) & 3, 0, v) ? 1549 + (k << 3) : 1405 + (k << 2));
462 }
463 break;
464 default: // case 13
465 switch (abs(MC33_faceTests(u.f, 165, v))) {
466 case 0:
467 k = ((u.f[1] < 0) << 1) | (u.f[5] < 0);
468 if (u.f[0] * u.f[1] == u.f[5]) //13.4
469 pcase += 2157 + 12 * k;
470 else {
471 c = MC33_interiorTest(k, 1, v); // 13.5.1 if c == 0 else 13.5.2
472 pcase += 2285 + (c ? 10 * k - 40 * c : 6 * k);
473 }
474 break;
475 case 2: //13.3
476 pcase += 1917 + 10 * ((u.f[0] < 0 ? u.f[2] > 0: 12 + (u.f[2] < 0)) + (u.f[1] < 0 ? u.f[3] < 0 : 6 + (u.f[3] > 0)));
477 if (u.f[4] > 0)
478 pcase += 30;
479 break;
480 case 4: //13.2
481 k = 21 + 11 * u.f[0] + 4 * u.f[1] + 3 * u.f[2] + 2 * u.f[3] + u.f[4];
482 if (k >> 4)
483 k -= (k & 32 ? 20 : 10);
484 pcase += 1845 + 3 * k;
485 break;
486 default: //13.1
487 pcase += 1839 + 2 * u.f[0];
488 }
489 }
490 while (i) {
491 i = *(++pcase);
492 for (k = 3; k;) {
493 c = i & 0x0F;
494 i >>= 4;
495 if (IsInvalid(p[cast<uint01>(c)])) {
496 switch (c) { // the vertices r[3] and normals (r + 3)[3] are calculated here
497 case 0:
498 if (z || x)
499 p[0] = M->Dy[y][x];
500 else {
501 if (v[0] == 0) {
502 if (IsValid(p[3]))
503 p[0] = p[3];
504 else if (IsValid(p[8]))
505 p[0] = p[8];
506 else if (y && signbf(v[3]))
507 p[0] = M->Lz[y][0];
508 else if (y && signbf(v[4]))
509 p[0] = M->Dx[y][0];
510 else if (y ? signbf(M->iso - M->values[0][y - 1][0]) : 0)
511 p[0] = M->Dy[y - 1][0];
512 else
513 p[0] = M->surfint(0, y, 0, u.r);
514 }
515 else if (v[1] == 0) {
516 if (IsValid(p[1]))
517 p[0] = p[1];
518 else if (IsValid(p[9]))
519 p[0] = p[9];
520 else
521 p[0] = M->surfint(0, y + 1, 0, u.r);
522 }
523 else {
524 t = v[0] / (v[0] - v[1]);
525 u.r[0] = u.r[2] = 0;
526 u.r[1] = y + t;
527 u.r[3] = (v[4] - v[0]) * (1 - t) + (v[5] - v[1]) * t;
528 u.r[4] = v[1] - v[0];
529 u.r[5] = (v[3] - v[0]) * (1 - t) + (v[2] - v[1]) * t;
530 p[0] = M->store(u.r, combine(M->nodes[0][y][0], 1.0f - t, M->nodes[0][y + 1][0], t));
531 }
532 M->Dy[y][0] = p[0];
533 }
534 break;
535 case 1:
536 if (x)
537 p[1] = M->Lz[y + 1][x];
538 else {
539 if (v[1] == 0) {
540 if (IsValid(p[0]))
541 p[1] = p[0];
542 else if (IsValid(p[9]))
543 p[1] = p[9];
544 else if (z && signbf(v[0]))
545 p[1] = M->Dy[y][0];
546 else if (z && signbf(v[5]))
547 p[1] = M->Dx[y + 1][0];
548 else if (z && y + 1 < M->ny ? signbf(M->iso - M->values[z][y + 2][0]) : 0)
549 p[1] = M->Dy[y + 1][0];
550 else if (z ? signbf(M->iso - M->values[z - 1][y + 1][0]) : 0)
551 p[1] = M->Lz[y + 1][0]; // value of previous slice
552 else
553 p[1] = M->surfint(0, y + 1, z, u.r);
554 }
555 else if (v[2] == 0) {
556 if (IsValid(p[10]))
557 p[1] = p[10];
558 else if (IsValid(p[2]))
559 p[1] = p[2];
560 else
561 p[1] = M->surfint(0, y + 1, z + 1, u.r);
562 }
563 else {
564 t = v[1] / (v[1] - v[2]);
565 u.r[0] = 0;
566 u.r[1] = cast<fltp04>(y + 1);
567 u.r[2] = z + t;
568 u.r[3] = (v[5] - v[1]) * (1 - t) + (v[6] - v[2]) * t;
569 u.r[4] = (y + 1 < M->ny ? 0.5f * ((M->values[z][y][0] - M->values[z][y + 2][0]) * (1 - t)
570 + (M->values[z + 1][y][0] - M->values[z + 1][y + 2][0]) * t) :
571 (v[1] - v[0]) * (1 - t) + (v[2] - v[3]) * t);
572 u.r[5] = v[2] - v[1];
573 p[1] = M->store(u.r, combine(M->nodes[z][y + 1][0], 1.0f - t, M->nodes[z + 1][y + 1][0], 1.0f));
574 }
575 M->Lz[y + 1][0] = p[1];
576 }
577 break;
578 case 2:
579 if (x)
580 p[2] = M->Uy[y][x];
581 else {
582 if (v[3] == 0) {
583 if (IsValid(p[3]))
584 p[2] = p[3];
585 else if (IsValid(p[11]))
586 p[2] = p[11];
587 else if (y && signbf(v[0]))
588 p[2] = M->Lz[y][0];
589 else if (y && signbf(v[7]))
590 p[2] = M->Ux[y][0];
591 else if (y ? signbf(M->iso - M->values[z + 1][y - 1][0]) : 0)
592 p[2] = M->Uy[y - 1][0];
593 else
594 p[2] = M->surfint(0, y, z + 1, u.r);
595 }
596 else if (v[2] == 0) {
597 if (IsValid(p[10]))
598 p[2] = p[10];
599 else if (IsValid(p[1]))
600 p[2] = p[1];
601 else
602 p[2] = M->surfint(0, y + 1, z + 1, u.r);
603 }
604 else {
605 t = v[3] / (v[3] - v[2]);
606 u.r[0] = 0;
607 u.r[2] = cast<fltp04>(z + 1);
608 u.r[1] = y + t;
609 u.r[3] = (v[7] - v[3]) * (1 - t) + (v[6] - v[2]) * t;
610 u.r[4] = v[2] - v[3];
611 u.r[5] = (z + 1 < M->nz ? 0.5f * ((M->values[z][y][0] - M->values[z + 2][y][0]) * (1 - t)
612 + (M->values[z][y + 1][0] - M->values[z + 2][y + 1][0]) * t) :
613 (v[3] - v[0]) * (1 - t) + (v[2] - v[1]) * t);
614 p[2] = M->store(u.r, combine(M->nodes[z + 1][y][0], 1.0f - t, M->nodes[z + 1][y + 1][0], t));
615 }
616 M->Uy[y][0] = p[2];
617 }
618 break;
619 case 3:
620 if (y || x)
621 p[3] = M->Lz[y][x];
622 else {
623 if (v[0] == 0) {
624 if (IsValid(p[0]))
625 p[3] = p[0];
626 else if (IsValid(p[8]))
627 p[3] = p[8];
628 else if (z && signbf(v[1]))
629 p[3] = M->Dy[0][0];
630 else if (z && signbf(v[4]))
631 p[3] = M->Dx[0][0];
632 else if (z ? signbf(M->iso - M->values[z - 1][0][0]) : 0)
633 p[3] = M->Lz[0][0]; // value of previous slice
634 else
635 p[3] = M->surfint(0, 0, z, u.r);
636 }
637 else if (v[3] == 0) {
638 if (IsValid(p[2]))
639 p[3] = p[2];
640 else if (IsValid(p[11]))
641 p[3] = p[11];
642 else
643 p[3] = M->surfint(0, 0, z + 1, u.r);
644 }
645 else {
646 t = v[0] / (v[0] - v[3]);
647 u.r[0] = u.r[1] = 0;
648 u.r[2] = z + t;
649 u.r[3] = (v[4] - v[0]) * (1 - t) + (v[7] - v[3]) * t;
650 u.r[4] = (v[1] - v[0]) * (1 - t) + (v[2] - v[3]) * t;
651 u.r[5] = v[3] - v[0];
652 p[3] = M->store(u.r, combine(M->nodes[z][0][0], 1.0f - t, M->nodes[z + 1][0][0], t));
653 }
654 M->Lz[0][0] = p[3];
655 }
656 break;
657 case 4:
658 if (z)
659 p[4] = M->Dy[y][x + 1];
660 else {
661 if (v[4] == 0) {
662 if (IsValid(p[8]))
663 p[4] = p[8];
664 //TODO: TP Check this?
665 //else if (p[7] != FF)
666 // p[4] = p[7];
667 else if (y && signbf(v[0]))
668 p[4] = M->Dx[y][x];
669 else if (y && signbf(v[7]))
670 p[4] = M->Lz[y][x + 1];
671 else if (y ? signbf(M->iso - M->values[0][y - 1][x + 1]) : 0)
672 p[4] = M->Dy[y - 1][x + 1];
673 else if (y && x + 1 < M->nx ? signbf(M->iso - M->values[0][y][x + 2]) : 0)
674 p[4] = M->Dx[y][x + 1];
675 else
676 p[4] = M->surfint(x + 1, y, 0, u.r);
677 }
678 else if (v[5] == 0) {
679 if (IsValid(p[5]))
680 p[4] = p[5];
681 else if (IsValid(p[9]))
682 p[4] = p[9];
683 else
684 p[4] = M->surfint(x + 1, y + 1, 0, u.r);
685 }
686 else {
687 t = v[4] / (v[4] - v[5]);
688 u.r[0] = cast<fltp04>(x + 1);
689 u.r[2] = 0.0f;
690 u.r[1] = cast<fltp04>(y) + t;
691 u.r[3] = (x + 1 < M->nx ? 0.5f * ((M->values[0][y][x] - M->values[0][y][x + 2]) * (1 - t)
692 + (M->values[0][y + 1][x] - M->values[0][y + 1][x + 2]) * t) :
693 (v[4] - v[0]) * (1 - t) + (v[5] - v[1]) * t);
694 u.r[4] = v[5] - v[4];
695 u.r[5] = (v[7] - v[4]) * (1 - t) + (v[6] - v[5]) * t;
696 p[4] = M->store(u.r, combine(M->nodes[0][y][x + 1], 1.0f - t, M->nodes[0][y + 1][x + 1], t));
697 }
698 M->Dy[y][x + 1] = p[4];
699 }
700 break;
701 case 5:
702 if (v[5] == 0) {
703 if (z) {
704 if (signbf(v[4]))
705 p[5] = p[4] = M->Dy[y][x + 1];
706 else if (signbf(v[1]))
707 p[5] = p[9] = M->Dx[y + 1][x];
708 else if (x + 1 < M->nx ? signbf(M->iso - M->values[z][y + 1][x + 2]) : 0)
709 p[5] = M->Dx[y + 1][x + 1];
710 else if (y + 1 < M->ny ? signbf(M->iso - M->values[z][y + 2][x + 1]) : 0)
711 p[5] = M->Dy[y + 1][x + 1];
712 else if (signbf(M->iso - M->values[z - 1][y + 1][x + 1]))
713 p[5] = M->Lz[y + 1][x + 1]; // value of previous slice
714 else
715 p[5] = M->surfint(x + 1, y + 1, z, u.r);
716 }
717 else
718 p[5] = M->surfint(x + 1, y + 1, 0, u.r);
719 }
720 else if (v[6] == 0)
721 p[5] = M->surfint(x + 1, y + 1, z + 1, u.r);
722 else {
723 t = v[5] / (v[5] - v[6]);
724 u.r[0] = cast<fltp04>(x + 1);
725 u.r[1] = cast<fltp04>(y + 1);
726 u.r[2] = cast<fltp04>(z) + t;
727 u.r[3] = (x + 1 < M->nx ? 0.5f * ((M->values[z][y + 1][x] - M->values[z][y + 1][x + 2]) * (1 - t)
728 + (M->values[z + 1][y + 1][x] - M->values[z + 1][y + 1][x + 2]) * t) :
729 (v[5] - v[1]) * (1 - t) + (v[6] - v[2]) * t);
730 u.r[4] = (y + 1 < M->ny ? 0.5f * ((M->values[z][y][x + 1] - M->values[z][y + 2][x + 1]) * (1 - t)
731 + (M->values[z + 1][y][x + 1] - M->values[z + 1][y + 2][x + 1]) * t) :
732 (v[5] - v[4]) * (1 - t) + (v[6] - v[7]) * t);
733 u.r[5] = v[6] - v[5];
734 p[5] = M->store(u.r, combine(M->nodes[z][y + 1][x + 1], 1.0f - t, M->nodes[z + 1][y + 1][x + 1], t));
735 }
736 M->Lz[y + 1][x + 1] = p[5];
737 break;
738 case 6:
739 if (v[7] == 0) {
740 if (y) {
741 if (signbf(v[3]))
742 p[6] = p[11] = M->Ux[y][x];
743 else if (signbf(v[4]))
744 p[6] = p[7] = M->Lz[y][x + 1];
745 else if (signbf(M->iso - M->values[z + 1][y - 1][x + 1]))
746 p[6] = M->Uy[y - 1][x + 1];
747 else if (x + 1 < M->nx ? signbf(M->iso - M->values[z + 1][y][x + 2]) : 0)
748 p[6] = M->Ux[y][x + 1];
749 else
750 p[6] = M->surfint(x + 1, y, z + 1, u.r);
751 }
752 else if (IsValid(p[11]))
753 p[6] = p[11];
754 //TODO: TP Check this?
755 //else if (p[7] != FF)
756 // p[6] = p[7];
757 else
758 p[6] = M->surfint(x + 1, 0, z + 1, u.r);
759 }
760 else if (v[6] == 0) {
761 if (IsInvalid(p[5]))
762 p[6] = (IsInvalid(p[10]) ? M->surfint(x + 1, y + 1, z + 1, u.r) : p[10]);
763 else
764 p[6] = p[5];
765 }
766 else {
767 t = v[7] / (v[7] - v[6]);
768 u.r[0] = cast<fltp04>(x + 1);
769 u.r[1] = y + t;
770 u.r[2] = cast<fltp04>(z + 1);
771 u.r[3] = (x + 1 < M->nx ? 0.5f * ((M->values[z + 1][y][x] - M->values[z + 1][y][x + 2]) * (1 - t)
772 + (M->values[z + 1][y + 1][x] - M->values[z + 1][y + 1][x + 2]) * t) :
773 (v[7] - v[3]) * (1 - t) + (v[6] - v[2]) * t);
774 u.r[4] = v[6] - v[7];
775 u.r[5] = (z + 1 < M->nz ? 0.5f * ((M->values[z][y][x + 1] - M->values[z + 2][y][x + 1]) * (1 - t)
776 + (M->values[z][y + 1][x + 1] - M->values[z + 2][y + 1][x + 1]) * t) :
777 (v[7] - v[4]) * (1 - t) + (v[6] - v[5]) * t);
778 p[6] = M->store(u.r, combine(M->nodes[z + 1][y][x + 1], 1.0f - t, M->nodes[z + 1][y + 1][x + 1], t));
779 }
780 M->Uy[y][x + 1] = p[6];
781 break;
782 case 7:
783 if (y)
784 p[7] = M->Lz[y][x + 1];
785 else {
786 if (v[4] == 0) {
787 if (IsValid(p[8]))
788 p[7] = p[8];
789 else if (IsValid(p[4]))
790 p[7] = p[4];
791 else if (z && signbf(v[0]))
792 p[7] = M->Dx[0][x];
793 else if (z && signbf(v[5]))
794 p[7] = M->Dy[0][x + 1];
795 else if (z && x + 1 < M->nx ? signbf(M->iso - M->values[z][0][x + 2]) : 0)
796 p[7] = M->Dx[0][x + 1];
797 else if (z ? signbf(M->iso - M->values[z - 1][0][x + 1]) : 0)
798 p[7] = M->Lz[0][x + 1]; // value of previous slice
799 else
800 p[7] = M->surfint(x + 1, 0, z, u.r);
801 }
802 else if (v[7] == 0) {
803 if (IsValid(p[6]))
804 p[7] = p[6];
805 else if (IsValid(p[11]))
806 p[7] = p[11];
807 else
808 p[7] = M->surfint(x + 1, 0, z + 1, u.r);
809 }
810 else {
811 t = v[4] / (v[4] - v[7]);
812 u.r[0] = cast<fltp04>(x + 1);
813 u.r[1] = 0;
814 u.r[2] = z + t;
815 u.r[3] = (x + 1 < M->nx ? 0.5f * ((M->values[z][0][x] - M->values[z][0][x + 2]) * (1 - t)
816 + (M->values[z + 1][0][x] - M->values[z + 1][0][x + 2]) * t) :
817 (v[4] - v[0]) * (1 - t) + (v[7] - v[3]) * t);
818 u.r[4] = (v[5] - v[4]) * (1 - t) + (v[6] - v[7]) * t;
819 u.r[5] = v[7] - v[4];
820 p[7] = M->store(u.r, combine(M->nodes[z][0][x + 1], 1.0f - t, M->nodes[z + 1][0][x + 1], t));
821 }
822 M->Lz[0][x + 1] = p[7];
823 }
824 break;
825 case 8:
826 if (z || y)
827 p[8] = M->Dx[y][x];
828 else {
829 if (v[0] == 0) {
830 if (IsValid(p[3]))
831 p[8] = p[3];
832 else if (IsValid(p[0]))
833 p[8] = p[0];
834 else if (x && signbf(v[3]))
835 p[8] = M->Lz[0][x];
836 else if (x && signbf(v[1]))
837 p[8] = M->Dy[0][x];
838 else if (x ? signbf(M->iso - M->values[0][0][x - 1]) : 0)
839 p[8] = M->Dx[0][x - 1];
840 else
841 p[8] = M->surfint(x, 0, 0, u.r);
842 }
843 else if (v[4] == 0) {
844 if (IsValid(p[4]))
845 p[8] = p[4];
846 else if (IsValid(p[7]))
847 p[8] = p[7];
848 else
849 p[8] = M->surfint(x + 1, 0, 0, u.r);
850 }
851 else {
852 t = v[0] / (v[0] - v[4]);
853 u.r[1] = u.r[2] = 0;
854 u.r[0] = x + t;
855 u.r[3] = v[4] - v[0];
856 u.r[4] = (v[1] - v[0]) * (1 - t) + (v[5] - v[4]) * t;
857 u.r[5] = (v[3] - v[0]) * (1 - t) + (v[7] - v[4]) * t;
858 p[8] = M->store(u.r, combine(M->nodes[0][0][x], 1.0f - t, M->nodes[0][0][x + 1], t));
859 }
860 M->Dx[0][x] = p[8];
861 }
862 break;
863 case 9:
864 if (z)
865 p[9] = M->Dx[y + 1][x];
866 else {
867 if (v[1] == 0) {
868 if (IsValid(p[0]))
869 p[9] = p[0];
870 //else if (p[1] != FF)
871 // p[9] = p[1];
872 else if (x && signbf(v[0]))
873 p[9] = M->Dy[y][x];
874 else if (x && signbf(v[2]))
875 p[9] = M->Lz[y + 1][x];
876 else if (x ? signbf(M->iso - M->values[0][y + 1][x - 1]) : 0)
877 p[9] = M->Dx[y + 1][x - 1];
878 else
879 p[9] = M->surfint(x, y + 1, 0, u.r);
880 }
881 else if (v[5] == 0) {
882 if (IsValid(p[5]))
883 p[9] = p[5];
884 else if (IsValid(p[4]))
885 p[9] = p[4];
886 else
887 p[9] = M->surfint(x + 1, y + 1, 0, u.r);
888 }
889 else {
890 t = v[1] / (v[1] - v[5]);
891 u.r[1] = cast<fltp04>(y + 1);
892 u.r[2] = 0.0f;
893 u.r[0] = cast<fltp04>(x) + t;
894 u.r[3] = v[5] - v[1];
895 u.r[4] = (y + 1 < M->ny ? 0.5f * ((M->values[0][y][x] - M->values[0][y + 2][x]) * (1 - t)
896 + (M->values[0][y][x + 1] - M->values[0][y + 2][x + 1]) * t) :
897 (v[1] - v[0]) * (1 - t) + (v[5] - v[4]) * t);
898 u.r[5] = (v[2] - v[1]) * (1 - t) + (v[6] - v[5]) * t;
899 p[9] = M->store(u.r, combine(M->nodes[0][y + 1][x], 1.0f - t, M->nodes[0][y + 1][x + 1], t));
900 }
901 M->Dx[y + 1][x] = p[9];
902 }
903 break;
904 case 10:
905 if (v[2] == 0) {
906 if (x) {
907 if (signbf(v[1]))
908 p[10] = p[1] = M->Lz[y + 1][x];
909 else if (signbf(v[3]))
910 p[10] = p[2] = M->Uy[y][x];
911 else if (signbf(M->iso - M->values[z + 1][y + 1][x - 1]))
912 p[10] = M->Ux[y + 1][x - 1];
913 else
914 p[10] = M->surfint(x, y + 1, z + 1, u.r);
915 }
916 else if (IsValid(p[2]))
917 p[10] = p[2];
918 //else if (p[1] != FF)
919 // p[10] = p[1];
920 else
921 p[10] = M->surfint(0, y + 1, z + 1, u.r);
922 }
923 else if (v[6] == 0) {
924 if (IsInvalid(p[5]))
925 p[10] = (IsInvalid(p[6]) ? M->surfint(x + 1, y + 1, z + 1, u.r) : p[6]);
926 else
927 p[10] = p[5];
928 }
929 else {
930 t = v[2] / (v[2] - v[6]);
931 u.r[0] = x + t;
932 u.r[1] = cast<fltp04>(y + 1);
933 u.r[2] = cast<fltp04>(z + 1);
934 u.r[3] = v[6] - v[2];
935 u.r[4] = (y + 1 < M->ny ? 0.5f * ((M->values[z + 1][y][x] - M->values[z + 1][y + 2][x]) * (1 - t)
936 + (M->values[z + 1][y][x + 1] - M->values[z + 1][y + 2][x + 1]) * t) :
937 (v[2] - v[3]) * (1 - t) + (v[6] - v[7]) * t);
938 u.r[5] = (z + 1 < M->nz ? 0.5f * ((M->values[z][y + 1][x] - M->values[z + 2][y + 1][x]) * (1 - t)
939 + (M->values[z][y + 1][x + 1] - M->values[z + 2][y + 1][x + 1]) * t) :
940 (v[2] - v[1]) * (1 - t) + (v[6] - v[5]) * t);
941 p[10] = M->store(u.r, combine(M->nodes[z + 1][y + 1][x], 1.0f - t, M->nodes[z + 1][y + 1][x + 1], t));
942 }
943 M->Ux[y + 1][x] = p[10];
944 break;
945 case 11:
946 if (y)
947 p[11] = M->Ux[y][x];
948 else {
949 if (v[3] == 0) {
950 if (IsValid(p[3]))
951 p[11] = p[3];
952 else if (IsValid(p[2]))
953 p[11] = p[2];
954 else if (x && signbf(v[0]))
955 p[11] = M->Lz[0][x];
956 else if (x && signbf(v[2]))
957 p[11] = M->Uy[0][x];
958 else if (x ? signbf(M->iso - M->values[z + 1][0][x - 1]) : 0)
959 p[11] = M->Ux[0][x - 1];
960 else
961 p[11] = M->surfint(x, 0, z + 1, u.r);
962 }
963 else if (v[7] == 0) {
964 if (IsValid(p[6]))
965 p[11] = p[6];
966 else if (IsValid(p[7]))
967 p[11] = p[7];
968 else
969 p[11] = M->surfint(x + 1, 0, z + 1, u.r);
970 }
971 else {
972 t = v[3] / (v[3] - v[7]);
973 u.r[1] = 0;
974 u.r[2] = cast<fltp04>(z + 1);
975 u.r[0] = x + t;
976 u.r[3] = v[7] - v[3];
977 u.r[4] = (v[2] - v[3]) * (1 - t) + (v[6] - v[7]) * t;
978 u.r[5] = (z + 1 < M->nz ? 0.5f * ((M->values[z][0][x] - M->values[z + 2][0][x]) * (1 - t)
979 + (M->values[z][0][x + 1] - M->values[z + 2][0][x + 1]) * t) :
980 (v[3] - v[0]) * (1 - t) + (v[7] - v[4]) * t);
981 p[11] = M->store(u.r, combine(M->nodes[z + 1][0][x], 1.0f - t, M->nodes[z + 1][0][x + 1], t));
982 }
983 M->Ux[0][x] = p[11];
984 }
985 break;
986 default:
987 u.r[0] = x + 0.5f;
988 u.r[1] = y + 0.5f;
989 u.r[2] = z + 0.5f;
990 u.r[3] = v[4] + v[5] + v[6] + v[7] - v[0] - v[1] - v[2] - v[3];
991 u.r[4] = v[1] + v[2] + v[5] + v[6] - v[0] - v[3] - v[4] - v[7];
992 u.r[5] = v[2] + v[3] + v[6] + v[7] - v[0] - v[1] - v[4] - v[5];
993 p[12] = M->store(u.r, combine(M->nodes[z][y][x], 0.5f, M->nodes[z + 1][y + 1][x + 1], 0.5f));
994 }
995 }
996 ti[--k] = p[cast<uint01>(c)];//now ti contains the vertex indices of the triangle
997 }
998 if (ti[0] != ti[1] && ti[0] != ti[2] && ti[1] != ti[2])//to avoid zero area triangles
999 {
1000 if (M->nT == M->capt)
1001 {
1002 uint04(*pt)[3] = M->T;
1003 M->capt <<= 1;
1004 M->T = (uint04(*)[3])realloc(pt, M->capt * 3 * sizeof(int));
1005 }
1006 uint04* vp = M->T[M->nT++];
1007#ifndef MC_Normal_neg
1008 * vp = ti[!m]; *(++vp) = ti[m]; *(++vp) = ti[2];
1009#else
1010 * vp = ti[m]; *(++vp) = ti[!m]; *(++vp) = ti[2];
1011#endif
1012 }
1013 }
1014 }
1015 class CubeTriangulation
1016 {
1017 public:
1018 /******************************************************************
1019 Calculates the isosurface (iso is the isovalue) using a MC33 structure
1020 pointed by M. The return value is a pointer to surface. The pointer will
1021 be NULL if there is not enough memory. The isovalue is stored in the
1022 member user.f[3] of surface struct.*/
1023 template<class t_point_type>
1024 static void calculate_isosurface(MC33<t_point_type>& M, float iso)
1025 {
1026 uint04 x, y, z, Nx = M.nx;
1027 float Vt[12];
1028 float* v1 = Vt, * v2 = Vt + 4;
1029 const GRD_data_type*** F = M.values, ** F0, ** F1, * V00, * V01, * V11, * V10;
1030 M.nT = M.nV = 0U;
1031 if (M.capt == 0)
1032 {
1033 M.capt = 4096;
1034 M.T = (uint04(*)[3])malloc(3 * 4096 * sizeof(int));
1035 }
1036 if (M.capv == 0)
1037 {
1038 M.capv = 4096;
1039 M.N = (float(*)[3])malloc(3 * 4096 * sizeof(float));
1040 M.V = (float(*)[3])malloc(3 * 4096 * sizeof(float));
1041 M.C = (t_point_type(*))malloc(4096 * sizeof(t_point_type));
1042 }
1043 M.iso = iso;
1044 for (z = 0; z != M.nz; ++z)
1045 {
1046 F0 = *F;
1047 F1 = *(++F);
1048 for (y = 0; y != M.ny; ++y)
1049 {
1050 V00 = *F0;
1051 V01 = *(++F0);
1052 V10 = *F1;
1053 V11 = *(++F1);
1054 v2[0] = iso - *V00;//the difference was inverted to use signbit function
1055 v2[1] = iso - *V01;
1056 v2[2] = iso - *V11;
1057 v2[3] = iso - *V10;
1058 //the eight least significant bits of i correspond to vertex indices. (x...x01234567)
1059 //If the bit is 1 then the vertex value is greater than zero.
1060 uint04 i = signbf(v2[3]) != 0;
1061 if (signbf(v2[2])) i |= 2;
1062 if (signbf(v2[1])) i |= 4;
1063 if (signbf(v2[0])) i |= 8;
1064 for (x = 0; x != Nx; ++x)
1065 {
1066 {float* P = v1; v1 = v2; v2 = P; }//v1 and v2 are exchanged
1067 v2[0] = iso - *(++V00);
1068 v2[1] = iso - *(++V01);
1069 v2[2] = iso - *(++V11);
1070 v2[3] = iso - *(++V10);
1071 i = ((i & 0x0F) << 4) | (signbf(v2[3]) != 0);
1072 if (signbf(v2[2])) i |= 2;
1073 if (signbf(v2[1])) i |= 4;
1074 if (signbf(v2[0])) i |= 8;
1075 if (i && i ^ 0xFF)
1076 {
1077 if (v1 > v2)
1078 {
1079 float* t = v2;
1080 float* s = t + 8;
1081 *s = *t;
1082 *(++s) = *(++t);
1083 *(++s) = *(++t);
1084 *(++s) = *(++t);
1085 }
1086 MC33_findCase(&M, x, y, z, i, v1);
1087 }
1088 }
1089 }
1090 {uint04** P = M.Dx; M.Dx = M.Ux; M.Ux = P; }//M.Dx and M.Ux are exchanged
1091 {uint04** P = M.Dy; M.Dy = M.Uy; M.Uy = P; }//M.Dy and M.Uy are exchanged
1092 }
1093 }
1094 };
1095}
Vector
A fixed-size array with N dimensions used as the basis for geometric and mathematical types.
Definition Vector.hpp:62
NDEVR
The primary namespace for the NDEVR SDK.
Definition ArialTileFetcherModule.h:35
fltp04
float fltp04
Defines an alias representing a 4 byte floating-point number Bit layout is as follows: -Sign: 1 bit a...
Definition BaseValues.hpp:128
IsValid
static constexpr bool IsValid(const Angle< t_type > &value)
Checks whether the given Angle holds a valid value.
Definition Angle.h:398
abs
constexpr Angle< t_angle_type > abs(const Angle< t_angle_type > &value)
Changes an input with a negative sign, to a positive sign.
Definition AngleFunctions.h:753
uint04
uint32_t uint04
-Defines an alias representing a 4 byte, unsigned integer -Can represent exact integer values 0 throu...
Definition BaseValues.hpp:97
uint01
uint8_t uint01
-Defines an alias representing a 1 byte, unsigned integer -Can represent exact integer values 0 throu...
Definition BaseValues.hpp:81
IsInvalid
static constexpr bool IsInvalid(const Angle< t_type > &value)
Checks whether the given Angle holds an invalid value.
Definition Angle.h:388
sqrt
t_type sqrt(const t_type &value)
Definition VectorFunctions.hpp:1181
cast
constexpr t_to cast(const Angle< t_from > &value)
Casts an Angle from one backing type to another.
Definition Angle.h:408
_GRD
3D grid data structure holding scalar field values and node data for marching cubes isosurface extrac...
Definition CubeTriangulation.h:51