COVISE Core
IsoSurfaceGPMUtil.h
Go to the documentation of this file.
1/* This file is part of COVISE.
2
3 You can use it under the terms of the GNU Lesser General Public License
4 version 2.1 or later, see lgpl-2.1.txt.
5
6 * License: LGPL 2+ */
7
8#ifndef _ISOSURFACEGPMUTIL_H_
9#define _ISOSURFACEGPMUTIL_H_
10
11#include <math.h>
12#include <vector>
13#include <algorithm>
14#include <cassert>
15#include <do/Triangulate.h>
16
17#if defined(__GNUC__) && (__GNUC__ > 4 || __GNUC__ == 4 && __GNUC_MINOR__ > 2)
18#pragma GCC diagnostic warning "-Wuninitialized"
19#endif
20
21/************************************************************/
22/* Structs and methods used in the IsoSurfaceComp module */
23/* for supporting Generalized Polyhedral Meshes */
24/************************************************************/
25namespace covise
26{
27
28typedef struct
29{
30 float v[3];
31 int dimension;
32} S_V_DATA;
33
34typedef struct
35{
36 float x;
37 float y;
38 float z;
39} EDGE_VECTOR;
40
41typedef struct
42{
43 int vertex1;
44 int vertex2;
45 int vertex3;
46
47 EDGE_VECTOR normal;
48} TRIANGLE;
49
50typedef struct
51{
52 bool intersection_at_vertex1;
53 bool intersection_at_vertex2;
54
55 int int_flag;
56 int vertex1;
57 int vertex2;
58
59 float data_vertex1;
60 float data_vertex2;
61
62 S_V_DATA data_vertex_int;
63 EDGE_VECTOR intersection;
64} ISOSURFACE_EDGE_INTERSECTION;
65
66typedef struct
67{
68 std::vector<int> ring;
69 std::vector<int> ring_index;
70 std::vector<int> polyhedron_faces;
71} CONTOUR;
72
73typedef std::vector<ISOSURFACE_EDGE_INTERSECTION> ISOSURFACE_EDGE_INTERSECTION_VECTOR;
74
75typedef std::vector<TRIANGLE> TESSELATION;
76
77typedef std::vector<int> CONNECTIVITY_VECTOR;
78typedef std::vector<int> POLYGON;
79typedef std::vector<int> PROCESSED_ELEMENTS;
80
81typedef std::vector<int>::iterator POLYGON_ITERATOR;
82typedef std::vector<int>::iterator ITERATOR;
83
84float dot_product(EDGE_VECTOR vector1, EDGE_VECTOR vector2)
85{
86 float scalar;
87 scalar = vector1.x * vector2.x + vector1.y * vector2.y + vector1.z * vector2.z;
88
89 return scalar;
90}
91
92EDGE_VECTOR cross_product(EDGE_VECTOR vector1, EDGE_VECTOR vector2)
93{
94 EDGE_VECTOR normal;
95
96 normal.x = (vector1.y * vector2.z) - (vector2.y * vector1.z);
97 normal.y = (vector2.x * vector1.z) - (vector1.x * vector2.z);
98 normal.z = (vector1.x * vector2.y) - (vector2.x * vector1.y);
99
100 return normal;
101}
102
103double vec_length(EDGE_VECTOR &vector)
104{
105 return sqrt(vector.x * vector.x + vector.y * vector.y + vector.z * vector.z);
106}
107
108ISOSURFACE_EDGE_INTERSECTION VertexInterpolate(float x1, float y1, float z1, float x2, float y2, float z2, float isovalue, float data1, float data2, int v1, int v2)
109{
110 ISOSURFACE_EDGE_INTERSECTION edge;
111
112 /************************************************/
113 /* Case 1: Isosurface passes through Vertex 1 */
114 /************************************************/
115
116 if (isovalue == data1)
117 {
118 edge.int_flag = 1;
119 edge.vertex1 = v1;
120 edge.vertex2 = v2;
121 edge.data_vertex1 = data1;
122 edge.data_vertex2 = data2;
123 edge.intersection_at_vertex1 = true;
124 edge.intersection_at_vertex2 = false;
125
126 if (x1 < x2)
127 edge.intersection.x = x1 + 0.01f * fabs(x2 - x1);
128 else if (x1 > x2)
129 edge.intersection.x = x1 - 0.01f * fabs(-x2 + x1);
130 else if (x1 == x2)
131 edge.intersection.x = x1;
132
133 if (y1 < y2)
134 edge.intersection.y = y1 + 0.01f * fabs(y2 - y1);
135 else if (y1 > y2)
136 edge.intersection.y = y1 - 0.01f * fabs(-y2 + y1);
137 else if (y1 == y2)
138 edge.intersection.y = y1;
139
140 if (z1 < z2)
141 edge.intersection.z = z1 + 0.01f * fabs(z2 - z1);
142 else if (z1 > z2)
143 edge.intersection.z = z1 - 0.01f * fabs(-z2 + z1);
144 else if (z1 == z2)
145 edge.intersection.z = z1;
146
147 return (edge);
148 }
149
150 /************************************************/
151 /* Case 2: Isosurface passes through Vertex 2 */
152 /************************************************/
153
154 if (isovalue == data2)
155 {
156 edge.int_flag = 1;
157 edge.vertex1 = v1;
158 edge.vertex2 = v2;
159 edge.data_vertex1 = data1;
160 edge.data_vertex2 = data2;
161 edge.intersection_at_vertex1 = false;
162 edge.intersection_at_vertex2 = true;
163
164 if (x1 < x2)
165 edge.intersection.x = x2 - 0.01f * fabs(x2 - x1);
166 else if (x1 > x2)
167 edge.intersection.x = x2 + 0.01f * fabs(-x2 + x1);
168 else if (x1 == x2)
169 edge.intersection.x = x2;
170
171 if (y1 < y2)
172 edge.intersection.y = y2 - 0.01f * fabs(y2 - y1);
173 else if (y1 > y2)
174 edge.intersection.y = y2 + 0.01f * fabs(-y2 + y1);
175 else if (y1 == y2)
176 edge.intersection.y = y2;
177
178 if (z1 < z2)
179 edge.intersection.z = z2 - 0.01f * fabs(z2 - z1);
180 else if (z1 > z2)
181 edge.intersection.z = z2 + 0.01f * fabs(-z2 + z1);
182 else if (z1 == z2)
183 edge.intersection.z = z2;
184
185 return (edge);
186 }
187
188 /****************************************************/
189 /* Case 3: Isosurface passes between both vertices */
190 /****************************************************/
191
192 else
193 {
194 edge.int_flag = 1;
195 edge.vertex1 = v1;
196 edge.vertex2 = v2;
197 edge.data_vertex1 = data1;
198 edge.data_vertex2 = data2;
199 edge.intersection_at_vertex1 = false;
200 edge.intersection_at_vertex2 = false;
201 edge.intersection.x = x1 + (x2 - x1) * ((isovalue - data1) / (data2 - data1));
202 edge.intersection.y = y1 + (y2 - y1) * ((isovalue - data1) / (data2 - data1));
203 edge.intersection.z = z1 + (z2 - z1) * ((isovalue - data1) / (data2 - data1));
204
205 return (edge);
206 }
207}
208
209bool test_intersection(ISOSURFACE_EDGE_INTERSECTION_VECTOR &intsec_vector, ISOSURFACE_EDGE_INTERSECTION &intsec, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool &improper_topology)
210{
211 /***********************************************************************************************************************/
212 /* Test if intersection has already been calculated */
213 /* */
214 /* According to O'Rourke a valid polyhedral surface must satisfy the following three conditions: */
215 /* */
216 /* a) components intersect "properly" */
217 /* b) local topology is "proper" */
218 /* c) global topology is "proper" */
219 /* */
220 /* It has been noted that some datasets contain cells that violate the second and consequently also third conditions, */
221 /* therefore for an adequate analysis it has to be determined if the cell is a "proper" polyhedron or not. Examples of */
222 /* these cells are transition cells found in multi-resolution grids. For more information refer to "Computational Geometry */
223 /* in C" by Joseph O'Rourke. */
224 /***********************************************************************************************************************/
225
226 for (std::vector<ISOSURFACE_EDGE_INTERSECTION>::iterator existent_intsec = intsec_vector.begin(); existent_intsec < intsec_vector.end(); ++existent_intsec)
227 {
228 /***********************************************/
229 /* Case 1: Local topology of the cell is "proper" */
230 /***********************************************/
231
232 /* Edge intersections share both vertices */
233 if (existent_intsec->vertex1 == intsec.vertex1 && existent_intsec->vertex2 == intsec.vertex2)
234 {
235 return true;
236 }
237
238 /* Edge intersections share both vertices (swapped) */
239 if (existent_intsec->vertex1 == intsec.vertex2 && existent_intsec->vertex2 == intsec.vertex1)
240 {
241 return true;
242 }
243
244 /*************************************************/
245 /* Case 2: Local topology of the cell is "improper" */
246 /*************************************************/
247
248 /* Check for T-vertices */
249 if (existent_intsec->vertex1 == intsec.vertex1 || existent_intsec->vertex2 == intsec.vertex2 || existent_intsec->vertex1 == intsec.vertex2 || existent_intsec->vertex2 == intsec.vertex1)
250 {
251 int vtx1;
252 int vtx2;
253
254 if (intsec.vertex1 == existent_intsec->vertex1 || intsec.vertex2 == existent_intsec->vertex2)
255 {
256 vtx1 = existent_intsec->vertex1;
257 vtx2 = existent_intsec->vertex2;
258 }
259
260 if (intsec.vertex1 == existent_intsec->vertex2 || intsec.vertex2 == existent_intsec->vertex1)
261 {
262 vtx1 = existent_intsec->vertex2;
263 vtx2 = existent_intsec->vertex1;
264 }
265
266 /* Check direction of edges */
267 EDGE_VECTOR d1;
268 d1.x = x_coord_in[vtx1] - x_coord_in[vtx2];
269 d1.y = y_coord_in[vtx1] - y_coord_in[vtx2];
270 d1.z = z_coord_in[vtx1] - z_coord_in[vtx2];
271
272 EDGE_VECTOR d2;
273 d2.x = x_coord_in[intsec.vertex1] - x_coord_in[intsec.vertex2];
274 d2.y = y_coord_in[intsec.vertex1] - y_coord_in[intsec.vertex2];
275 d2.z = z_coord_in[intsec.vertex1] - z_coord_in[intsec.vertex2];
276
277 double length1 = vec_length(d1);
278 double length2 = vec_length(d2);
279
280 if ((length1 > 0) && (length2 > 0))
281 {
282 double cosangle = dot_product(d1, d2) / (length1 * length2);
283
284 /******************************************************************************************************/
285 /* Cell is topologically "improper" --> Two edges have the same direction (and share a common vertex) */
286 /* */
287 /* The chosen tolerance has brought good results however the possibility always remains that a certain */
288 /* dataset contains extremely small triangles (polygons) which are practically degenerate. Under these */
289 /* circumstances a "proper" cell could be treated as "improper". */
290 /******************************************************************************************************/
291
292 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
293 {
294 improper_topology = true;
295 return true;
296 }
297 }
298 }
299 }
300
301 return false;
302}
303
304float map_to_isosurface(float coord_x1, float coord_x2, float coord_y1, float coord_y2, float coord_z1, float coord_z2, float coord_isox, float coord_isoy, float coord_isoz, float data_1, float data_2, bool int_vertex1, bool int_vertex2)
305{
306 float mapped_value;
307 float dist_x1x2;
308 float dist_x1xiso;
309
310 if (int_vertex1 == true)
311 {
312 mapped_value = data_1;
313 }
314
315 else if (int_vertex2 == true)
316 {
317 mapped_value = data_2;
318 }
319
320 else
321 {
322 dist_x1x2 = sqrt(pow(coord_x1 - coord_x2, 2.0f) + pow(coord_y1 - coord_y2, 2.0f) + pow(coord_z1 - coord_z2, 2.0f));
323
324 // Avoid division by zero
325 if (dist_x1x2 == 0)
326 {
327 mapped_value = data_1;
328 }
329
330 else
331 {
332 dist_x1xiso = sqrt(pow(coord_x1 - coord_isox, 2.0f) + pow(coord_y1 - coord_isoy, 2.0f) + pow(coord_z1 - coord_isoz, 2.0f));
333 mapped_value = data_1 + ((data_2 - data_1) / dist_x1x2) * dist_x1xiso;
334 }
335 }
336
337 return mapped_value;
338}
339
340ISOSURFACE_EDGE_INTERSECTION_VECTOR calculate_intersections(int num_elem_in, int *elem_in, int num_conn_in, int *conn_in, float *x_coord_in, float *y_coord_in, float *z_coord_in, float *isodata_in, float isovalue, float *sdata_in, float *udata_in, float *vdata_in, float *wdata_in, bool isomap, bool &improper_topology)
341{
342 int i;
343 int j;
344
345 float data_vertex1;
346 float data_vertex2;
347
348 ISOSURFACE_EDGE_INTERSECTION intersec;
349 ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector;
350
351 /* Avoid Unnecessary Reallocation */
352 intsec_vector.reserve(15);
353
354 for (i = 0; i < num_elem_in; i++)
355 {
356 if (i < num_elem_in - 1)
357 {
358 for (j = elem_in[i]; j <= elem_in[i + 1] - 1; j++)
359 {
360 if (j < elem_in[i + 1] - 1)
361 {
362 if (((isodata_in[conn_in[j]] <= isovalue) && (isovalue < isodata_in[conn_in[j + 1]])) || (isodata_in[conn_in[j]] > isovalue && isovalue >= isodata_in[conn_in[j + 1]]) || ((isodata_in[conn_in[j]] < isovalue) && (isovalue <= isodata_in[conn_in[j + 1]])) || (isodata_in[conn_in[j]] >= isovalue && isovalue > isodata_in[conn_in[j + 1]]))
363 {
364 data_vertex1 = isodata_in[conn_in[j]];
365 data_vertex2 = isodata_in[conn_in[j + 1]];
366
367 intersec = VertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j + 1]], isovalue, data_vertex1, data_vertex2, conn_in[j], conn_in[j + 1]);
368
369 if (intersec.int_flag == 1)
370 {
371 /* Check for repeated intersections */
372 if (!test_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in, improper_topology))
373 {
374 if (sdata_in != NULL)
375 {
376 intersec.data_vertex_int.dimension = 1;
377
378 /* Map scalar data only if required */
379 if (isomap == true)
380 {
381 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, sdata_in[conn_in[j]], sdata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
382 }
383 }
384
385 else if (udata_in != NULL)
386 {
387 intersec.data_vertex_int.dimension = 3;
388 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, udata_in[conn_in[j]], udata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
389 intersec.data_vertex_int.v[1] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, vdata_in[conn_in[j]], vdata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
390 intersec.data_vertex_int.v[2] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, wdata_in[conn_in[j]], wdata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
391 }
392
393 intsec_vector.push_back(intersec);
394 }
395 }
396 }
397 }
398
399 else if (j == elem_in[i + 1] - 1)
400 {
401 if ((isodata_in[conn_in[j]] <= isovalue && isovalue < isodata_in[conn_in[elem_in[i]]]) || (isodata_in[conn_in[j]] > isovalue && isovalue >= isodata_in[conn_in[elem_in[i]]]) || (isodata_in[conn_in[j]] < isovalue && isovalue <= isodata_in[conn_in[elem_in[i]]]) || (isodata_in[conn_in[j]] >= isovalue && isovalue > isodata_in[conn_in[elem_in[i]]]))
402 {
403 data_vertex1 = isodata_in[conn_in[j]];
404 data_vertex2 = isodata_in[conn_in[elem_in[i]]];
405
406 intersec = VertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[elem_in[i]]], isovalue, data_vertex1, data_vertex2, conn_in[j], conn_in[elem_in[i]]);
407
408 if (intersec.int_flag == 1)
409 {
410 /* Check for repeated intersections */
411 if (!test_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in, improper_topology))
412 {
413 if (sdata_in != NULL)
414 {
415 intersec.data_vertex_int.dimension = 1;
416
417 /* Map scalar data only if required */
418 if (isomap == true)
419 {
420 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, sdata_in[conn_in[j]], sdata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
421 }
422 }
423
424 else if (udata_in != NULL)
425 {
426 intersec.data_vertex_int.dimension = 3;
427 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, udata_in[conn_in[j]], udata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
428 intersec.data_vertex_int.v[1] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, vdata_in[conn_in[j]], vdata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
429 intersec.data_vertex_int.v[2] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, wdata_in[conn_in[j]], wdata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
430 }
431
432 intsec_vector.push_back(intersec);
433 }
434 }
435 }
436 }
437 }
438 }
439
440 else
441 {
442 for (j = elem_in[i]; j <= (num_conn_in - 1); j++)
443 {
444 if (j < (num_conn_in - 1))
445 {
446 if ((isodata_in[conn_in[j]] <= isovalue && isovalue < isodata_in[conn_in[j + 1]]) || (isodata_in[conn_in[j]] > isovalue && isovalue >= isodata_in[conn_in[j + 1]]) || (isodata_in[conn_in[j]] < isovalue && isovalue <= isodata_in[conn_in[j + 1]]) || (isodata_in[conn_in[j]] >= isovalue && isovalue > isodata_in[conn_in[j + 1]]))
447 {
448 data_vertex1 = isodata_in[conn_in[j]];
449 data_vertex2 = isodata_in[conn_in[j + 1]];
450
451 intersec = VertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j + 1]], isovalue, data_vertex1, data_vertex2, conn_in[j], conn_in[j + 1]);
452
453 if (intersec.int_flag == 1)
454 {
455 /* Check for repeated intersections */
456 if (!test_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in, improper_topology))
457 {
458 if (sdata_in != NULL)
459 {
460 intersec.data_vertex_int.dimension = 1;
461
462 /* Map scalar data only if required */
463 if (isomap == true)
464 {
465 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, sdata_in[conn_in[j]], sdata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
466 }
467 }
468
469 else if (udata_in != NULL)
470 {
471 intersec.data_vertex_int.dimension = 3;
472 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, udata_in[conn_in[j]], udata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
473 intersec.data_vertex_int.v[1] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, vdata_in[conn_in[j]], vdata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
474 intersec.data_vertex_int.v[2] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[j + 1]], y_coord_in[conn_in[j]], y_coord_in[conn_in[j + 1]], z_coord_in[conn_in[j]], z_coord_in[conn_in[j + 1]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, wdata_in[conn_in[j]], wdata_in[conn_in[j + 1]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
475 }
476
477 intsec_vector.push_back(intersec);
478 }
479 }
480 }
481 }
482
483 else if (j == (num_conn_in - 1))
484 {
485 if ((isodata_in[conn_in[j]] <= isovalue && isovalue < isodata_in[conn_in[elem_in[i]]]) || (isodata_in[conn_in[j]] > isovalue && isovalue >= isodata_in[conn_in[elem_in[i]]]) || (isodata_in[conn_in[j]] < isovalue && isovalue <= isodata_in[conn_in[elem_in[i]]]) || (isodata_in[conn_in[j]] >= isovalue && isovalue > isodata_in[conn_in[elem_in[i]]]))
486 {
487 data_vertex1 = isodata_in[conn_in[j]];
488 data_vertex2 = isodata_in[conn_in[elem_in[i]]];
489
490 intersec = VertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[elem_in[i]]], isovalue, data_vertex1, data_vertex2, conn_in[j], conn_in[elem_in[i]]);
491
492 if (intersec.int_flag == 1)
493 {
494 /* Check for repeated intersections */
495 if (!test_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in, improper_topology))
496 {
497 if (sdata_in != NULL)
498 {
499 intersec.data_vertex_int.dimension = 1;
500
501 /* Map scalar data only if required */
502 if (isomap == true)
503 {
504 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, sdata_in[conn_in[j]], sdata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
505 }
506 }
507
508 else if (udata_in != NULL)
509 {
510 intersec.data_vertex_int.dimension = 3;
511 intersec.data_vertex_int.v[0] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, udata_in[conn_in[j]], udata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
512 intersec.data_vertex_int.v[1] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, vdata_in[conn_in[j]], vdata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
513 intersec.data_vertex_int.v[2] = map_to_isosurface(x_coord_in[conn_in[j]], x_coord_in[conn_in[elem_in[i]]], y_coord_in[conn_in[j]], y_coord_in[conn_in[elem_in[i]]], z_coord_in[conn_in[j]], z_coord_in[conn_in[elem_in[i]]], intersec.intersection.x, intersec.intersection.y, intersec.intersection.z, wdata_in[conn_in[j]], wdata_in[conn_in[elem_in[i]]], intersec.intersection_at_vertex1, intersec.intersection_at_vertex2);
514 }
515
516 intsec_vector.push_back(intersec);
517 }
518 }
519 }
520 }
521 }
522 }
523 }
524
525 return intsec_vector;
526}
527
528int assign_int_index(ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, int edge_vertex1, int edge_vertex2)
529{
530 size_t i;
531 int index;
532
533 index = 0;
534
535 for (i = 0; i < intsec_vector.size(); i++)
536 {
537 if ((edge_vertex1 == intsec_vector[i].vertex1) && (edge_vertex2 == intsec_vector[i].vertex2))
538 {
539 index = (int)i;
540 break;
541 }
542
543 if ((edge_vertex2 == intsec_vector[i].vertex1) && (edge_vertex1 == intsec_vector[i].vertex2))
544 {
545 index = (int)i;
546 break;
547 }
548 }
549
550 return index;
551}
552
553bool find_intersection(ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, int &edge_vertex1, int &edge_vertex2, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool improper_topology, int &int_index)
554{
555 int i;
556 bool int_found;
557
558 int_found = false;
559
560 for (i = 0; i < ssize_t(intsec_vector.size()); i++)
561 {
562 /***********************************************/
563 /* Case 1: Local topology of the cell is "proper" */
564 /***********************************************/
565
566 /* Edge intersections share both vertices */
567 if ((edge_vertex1 == intsec_vector[i].vertex1) && (edge_vertex2 == intsec_vector[i].vertex2))
568 {
569 int_found = true;
570 int_index = i;
571 break;
572 }
573
574 /* Edge intersections share both vertices (swapped) */
575 else if ((edge_vertex2 == intsec_vector[i].vertex1) && (edge_vertex1 == intsec_vector[i].vertex2))
576 {
577 int_found = true;
578 int_index = i;
579 break;
580 }
581
582 /*************************************************/
583 /* Case 2: Local topology of the cell is "improper" */
584 /*************************************************/
585
586 // Check for T-vertices
587 else if (improper_topology)
588 {
589 if (edge_vertex1 == intsec_vector[i].vertex1 || edge_vertex2 == intsec_vector[i].vertex2 || edge_vertex2 == intsec_vector[i].vertex1 || edge_vertex1 == intsec_vector[i].vertex2)
590 {
591 int vtx1;
592 int vtx2;
593
594 if (edge_vertex1 == intsec_vector[i].vertex1 || edge_vertex2 == intsec_vector[i].vertex2)
595 {
596 vtx1 = intsec_vector[i].vertex1;
597 vtx2 = intsec_vector[i].vertex2;
598 }
599
600 if (edge_vertex1 == intsec_vector[i].vertex2 || edge_vertex2 == intsec_vector[i].vertex1)
601 {
602 vtx1 = intsec_vector[i].vertex2;
603 vtx2 = intsec_vector[i].vertex1;
604 }
605
606 // Check direction of edges
607 EDGE_VECTOR d1;
608 d1.x = x_coord_in[vtx1] - x_coord_in[vtx2];
609 d1.y = y_coord_in[vtx1] - y_coord_in[vtx2];
610 d1.z = z_coord_in[vtx1] - z_coord_in[vtx2];
611
612 EDGE_VECTOR d2;
613 d2.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
614 d2.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
615 d2.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
616
617 double length1 = vec_length(d1);
618 double length2 = vec_length(d2);
619
620 if ((length1 > 0) && (length2 > 0))
621 {
622 double cosangle = dot_product(d1, d2) / (length1 * length2);
623
624 // Cell is topologically "improper"
625 // Two edges have the same direction (and share a common vertex)
626 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
627 {
628 if (edge_vertex1 == intsec_vector[i].vertex1)
629 edge_vertex2 = intsec_vector[i].vertex2;
630 else if (edge_vertex2 == intsec_vector[i].vertex2)
631 edge_vertex1 = intsec_vector[i].vertex1;
632 else if (edge_vertex1 == intsec_vector[i].vertex2)
633 edge_vertex2 = intsec_vector[i].vertex1;
634 else if (edge_vertex2 == intsec_vector[i].vertex1)
635 edge_vertex1 = intsec_vector[i].vertex2;
636
637 int_found = true;
638 int_index = i;
639 break;
640 }
641 }
642 }
643 }
644 }
645
646 return int_found;
647}
648
649void find_current_face(CONTOUR &contour, ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, int &edge_vertex1, int &edge_vertex2, float &data_vertex1, float &data_vertex2, float *isodata_in, int *elem_in, int *conn_in, int *index_list, int *polygon_list, int num_coord_in, int num_conn_in, int num_elem_in, int &ring_counter, int &current_face, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool improper_topology, bool &abort_tracing_isocontour)
650{
651 bool adjacent_vertices;
652 bool T_vertices;
653
654 int i;
655 int j;
656 int k;
657 int face_flag;
658 int copy_current_face;
659 int neighbour_face1;
660 int neighbour_face2;
661 int next_vertex;
662 int previous_vertex;
663 //int vertex_index;
664 int new_edge_vertex;
665
666 ITERATOR it;
667
668 face_flag = 0;
669 copy_current_face = current_face;
670
671 /******************************************************************************/
672 /* Find the next face of the polyhedron to continue tracing the convex contour */
673 /******************************************************************************/
674
675 if ((edge_vertex1 < num_coord_in - 1) && (edge_vertex2 < num_coord_in - 1))
676 {
677 for (i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
678 {
679 adjacent_vertices = false;
680 neighbour_face1 = polygon_list[i];
681 for (j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
682 {
683 neighbour_face2 = polygon_list[j];
684 if (face_flag == 0)
685 {
686 /*********************************************************/
687 /* Search among the elements that contain both vertices */
688 /*********************************************************/
689
690 if (neighbour_face1 == neighbour_face2)
691 {
692 if (neighbour_face1 < num_elem_in - 1)
693 {
694 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
695 {
696 if (conn_in[k] == edge_vertex1)
697 {
698 //vertex_index = k;
699 if (k == elem_in[neighbour_face1])
700 {
701 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
702 next_vertex = conn_in[k + 1];
703 }
704
705 else if (k < elem_in[neighbour_face1 + 1] - 1)
706 {
707 previous_vertex = conn_in[k - 1];
708 next_vertex = conn_in[k + 1];
709 }
710
711 else if (k == elem_in[neighbour_face1 + 1] - 1)
712 {
713 previous_vertex = conn_in[k - 1];
714 next_vertex = conn_in[elem_in[neighbour_face1]];
715 }
716
717 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
718 {
719 adjacent_vertices = true;
720 break;
721 }
722 }
723 }
724 }
725
726 else
727 {
728 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
729 {
730 if (conn_in[k] == edge_vertex1)
731 {
732 //vertex_index = k;
733 if (k == elem_in[neighbour_face1])
734 {
735 previous_vertex = conn_in[num_conn_in - 1];
736 next_vertex = conn_in[k + 1];
737 }
738
739 else if (k < num_conn_in - 1)
740 {
741 previous_vertex = conn_in[k - 1];
742 next_vertex = conn_in[k + 1];
743 }
744
745 else if (k == num_conn_in - 1)
746 {
747 previous_vertex = conn_in[k - 1];
748 next_vertex = conn_in[elem_in[neighbour_face1]];
749 }
750
751 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
752 {
753 adjacent_vertices = true;
754 break;
755 }
756 }
757 }
758 }
759
760 /**************************************************/
761 /* Test if the element has already been processed */
762 /**************************************************/
763
764 if (adjacent_vertices == true)
765 {
766 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
767
768 /* The current face that contains the vertices has not been processed */
769 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
770 {
771 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
772 contour.polyhedron_faces.push_back(neighbour_face1);
773 current_face = neighbour_face1;
774 ring_counter++;
775 face_flag = 1;
776 break;
777 }
778 }
779 }
780 }
781 }
782
783 if (face_flag == 1)
784 {
785 break;
786 }
787 }
788
789 /**********************************/
790 /* No element has been found yet */
791 /**********************************/
792
793 if (face_flag == 0)
794 {
795 /**********************************************/
796 /* Case 1: The cell has an "improper" topology */
797 /**********************************************/
798
799 if (improper_topology)
800 {
801 T_vertices = false;
802
803 /**************************************************/
804 /* Search among the elements that share vertex1 */
805 /**************************************************/
806
807 if (!T_vertices)
808 {
809 for (i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
810 {
811 neighbour_face1 = polygon_list[i];
812 if (neighbour_face1 != current_face)
813 {
814 if (neighbour_face1 < num_elem_in - 1)
815 {
816 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
817 {
818 if (conn_in[k] == edge_vertex1)
819 {
820 //vertex_index = k;
821 if (k == elem_in[neighbour_face1])
822 {
823 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
824 next_vertex = conn_in[k + 1];
825 }
826
827 else if (k < elem_in[neighbour_face1 + 1] - 1)
828 {
829 previous_vertex = conn_in[k - 1];
830 next_vertex = conn_in[k + 1];
831 }
832
833 else if (k == elem_in[neighbour_face1 + 1] - 1)
834 {
835 previous_vertex = conn_in[k - 1];
836 next_vertex = conn_in[elem_in[neighbour_face1]];
837 }
838
839 /* Check for edges that have the same direction and share a common vertex */
840 if (previous_vertex != edge_vertex2 && next_vertex != edge_vertex2)
841 {
842 // Check direction of edges
843 EDGE_VECTOR d1;
844 d1.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
845 d1.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
846 d1.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
847
848 EDGE_VECTOR d2;
849 d2.x = x_coord_in[edge_vertex1] - x_coord_in[previous_vertex];
850 d2.y = y_coord_in[edge_vertex1] - y_coord_in[previous_vertex];
851 d2.z = z_coord_in[edge_vertex1] - z_coord_in[previous_vertex];
852
853 EDGE_VECTOR d3;
854 d3.x = x_coord_in[edge_vertex1] - x_coord_in[next_vertex];
855 d3.y = y_coord_in[edge_vertex1] - y_coord_in[next_vertex];
856 d3.z = z_coord_in[edge_vertex1] - z_coord_in[next_vertex];
857
858 double length1 = vec_length(d1);
859 double length2 = vec_length(d2);
860 double length3 = vec_length(d3);
861
862 if ((length1 > 0) && (length2 > 0))
863 {
864 double cosangle = dot_product(d1, d2) / (length1 * length2);
865
866 // Cell is topologically "improper"
867 // Two edges have the same direction (and share a common vertex)
868 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
869 {
870 T_vertices = true;
871 new_edge_vertex = previous_vertex;
872 break;
873 }
874 }
875
876 if ((length1 > 0) && (length3 > 0))
877 {
878 double cosangle = dot_product(d1, d3) / (length1 * length3);
879
880 // Cell is topologically "improper"
881 // Two edges have the same direction (and share a common vertex)
882 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
883 {
884 T_vertices = true;
885 new_edge_vertex = next_vertex;
886 break;
887 }
888 }
889 }
890 }
891 }
892 }
893
894 else
895 {
896 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
897 {
898 if (conn_in[k] == edge_vertex1)
899 {
900 //vertex_index = k;
901 if (k == elem_in[neighbour_face1])
902 {
903 previous_vertex = conn_in[num_conn_in - 1];
904 next_vertex = conn_in[k + 1];
905 }
906
907 else if (k < num_conn_in - 1)
908 {
909 previous_vertex = conn_in[k - 1];
910 next_vertex = conn_in[k + 1];
911 }
912
913 else if (k == num_conn_in - 1)
914 {
915 previous_vertex = conn_in[k - 1];
916 next_vertex = conn_in[elem_in[neighbour_face1]];
917 }
918
919 /* Check for edges that have the same direction and share a common vertex */
920 if (previous_vertex != edge_vertex2 && next_vertex != edge_vertex2)
921 {
922 // Check direction of edges
923 EDGE_VECTOR d1;
924 d1.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
925 d1.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
926 d1.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
927
928 EDGE_VECTOR d2;
929 d2.x = x_coord_in[edge_vertex1] - x_coord_in[previous_vertex];
930 d2.y = y_coord_in[edge_vertex1] - y_coord_in[previous_vertex];
931 d2.z = z_coord_in[edge_vertex1] - z_coord_in[previous_vertex];
932
933 EDGE_VECTOR d3;
934 d3.x = x_coord_in[edge_vertex1] - x_coord_in[next_vertex];
935 d3.y = y_coord_in[edge_vertex1] - y_coord_in[next_vertex];
936 d3.z = z_coord_in[edge_vertex1] - z_coord_in[next_vertex];
937
938 double length1 = vec_length(d1);
939 double length2 = vec_length(d2);
940 double length3 = vec_length(d3);
941
942 if ((length1 > 0) && (length2 > 0))
943 {
944 double cosangle = dot_product(d1, d2) / (length1 * length2);
945
946 // Cell is topologically "improper"
947 // Two edges have the same direction (and share a common vertex)
948 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
949 {
950 T_vertices = true;
951 new_edge_vertex = previous_vertex;
952 break;
953 }
954 }
955
956 if ((length1 > 0) && (length3 > 0))
957 {
958 double cosangle = dot_product(d1, d3) / (length1 * length3);
959
960 // Cell is topologically "improper"
961 // Two edges have the same direction (and share a common vertex)
962 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
963 {
964 T_vertices = true;
965 new_edge_vertex = next_vertex;
966 break;
967 }
968 }
969 }
970 }
971 }
972 }
973 }
974
975 /**************************************************/
976 /* Test if the element has already been processed */
977 /**************************************************/
978
979 if (T_vertices == true)
980 {
981 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
982
983 /* The current face that contains the vertices has not been processed */
984 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
985 {
986 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
987 contour.polyhedron_faces.push_back(neighbour_face1);
988 current_face = neighbour_face1;
989
990 /* Update edge vertices */
991 edge_vertex2 = new_edge_vertex;
992 data_vertex2 = isodata_in[edge_vertex2];
993 ring_counter++;
994 face_flag = 1;
995 break;
996 }
997 }
998 }
999 }
1000
1001 /**************************************************/
1002 /* Search among the elements that share vertex2 */
1003 /**************************************************/
1004
1005 if (!T_vertices)
1006 {
1007 for (j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
1008 {
1009 neighbour_face2 = polygon_list[j];
1010 if (neighbour_face2 != current_face)
1011 {
1012 if (neighbour_face2 < num_elem_in - 1)
1013 {
1014 for (k = elem_in[neighbour_face2]; k < elem_in[neighbour_face2 + 1]; k++)
1015 {
1016 if (conn_in[k] == edge_vertex2)
1017 {
1018 //vertex_index = k;
1019 if (k == elem_in[neighbour_face2])
1020 {
1021 previous_vertex = conn_in[elem_in[neighbour_face2 + 1] - 1];
1022 next_vertex = conn_in[k + 1];
1023 }
1024
1025 else if (k < elem_in[neighbour_face2 + 1] - 1)
1026 {
1027 previous_vertex = conn_in[k - 1];
1028 next_vertex = conn_in[k + 1];
1029 }
1030
1031 else if (k == elem_in[neighbour_face2 + 1] - 1)
1032 {
1033 previous_vertex = conn_in[k - 1];
1034 next_vertex = conn_in[elem_in[neighbour_face2]];
1035 }
1036
1037 /* Check for edges that have the same direction and share a common vertex */
1038 if (previous_vertex != edge_vertex1 && next_vertex != edge_vertex1)
1039 {
1040 // Check direction of edges
1041 EDGE_VECTOR d1;
1042 d1.x = x_coord_in[edge_vertex2] - x_coord_in[edge_vertex1];
1043 d1.y = y_coord_in[edge_vertex2] - y_coord_in[edge_vertex1];
1044 d1.z = z_coord_in[edge_vertex2] - z_coord_in[edge_vertex1];
1045
1046 EDGE_VECTOR d2;
1047 d2.x = x_coord_in[edge_vertex2] - x_coord_in[previous_vertex];
1048 d2.y = y_coord_in[edge_vertex2] - y_coord_in[previous_vertex];
1049 d2.z = z_coord_in[edge_vertex2] - z_coord_in[previous_vertex];
1050
1051 EDGE_VECTOR d3;
1052 d3.x = x_coord_in[edge_vertex2] - x_coord_in[next_vertex];
1053 d3.y = y_coord_in[edge_vertex2] - y_coord_in[next_vertex];
1054 d3.z = z_coord_in[edge_vertex2] - z_coord_in[next_vertex];
1055
1056 double length1 = vec_length(d1);
1057 double length2 = vec_length(d2);
1058 double length3 = vec_length(d3);
1059
1060 if ((length1 > 0) && (length2 > 0))
1061 {
1062 double cosangle = dot_product(d1, d2) / (length1 * length2);
1063
1064 // Cell is topologically "improper"
1065 // Two edges have the same direction (and share a common vertex)
1066 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1067 {
1068 T_vertices = true;
1069 new_edge_vertex = previous_vertex;
1070 break;
1071 }
1072 }
1073
1074 if ((length1 > 0) && (length3 > 0))
1075 {
1076 double cosangle = dot_product(d1, d3) / (length1 * length3);
1077
1078 // Cell is topologically "improper"
1079 // Two edges have the same direction (and share a common vertex)
1080 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1081 {
1082 T_vertices = true;
1083 new_edge_vertex = next_vertex;
1084 break;
1085 }
1086 }
1087 }
1088 }
1089 }
1090 }
1091
1092 else
1093 {
1094 for (k = elem_in[neighbour_face2]; k < num_conn_in; k++)
1095 {
1096 if (conn_in[k] == edge_vertex2)
1097 {
1098 //vertex_index = k;
1099 if (k == elem_in[neighbour_face2])
1100 {
1101 previous_vertex = conn_in[num_conn_in - 1];
1102 next_vertex = conn_in[k + 1];
1103 }
1104
1105 else if (k < num_conn_in - 1)
1106 {
1107 previous_vertex = conn_in[k - 1];
1108 next_vertex = conn_in[k + 1];
1109 }
1110
1111 else if (k == num_conn_in - 1)
1112 {
1113 previous_vertex = conn_in[k - 1];
1114 next_vertex = conn_in[elem_in[neighbour_face2]];
1115 }
1116
1117 /* Check for edges that have the same direction and share a common vertex */
1118 if (previous_vertex != edge_vertex1 && next_vertex != edge_vertex1)
1119 {
1120 // Check direction of edges
1121 EDGE_VECTOR d1;
1122 d1.x = x_coord_in[edge_vertex2] - x_coord_in[edge_vertex1];
1123 d1.y = y_coord_in[edge_vertex2] - y_coord_in[edge_vertex1];
1124 d1.z = z_coord_in[edge_vertex2] - z_coord_in[edge_vertex1];
1125
1126 EDGE_VECTOR d2;
1127 d2.x = x_coord_in[edge_vertex2] - x_coord_in[previous_vertex];
1128 d2.y = y_coord_in[edge_vertex2] - y_coord_in[previous_vertex];
1129 d2.z = z_coord_in[edge_vertex2] - z_coord_in[previous_vertex];
1130
1131 EDGE_VECTOR d3;
1132 d3.x = x_coord_in[edge_vertex2] - x_coord_in[next_vertex];
1133 d3.y = y_coord_in[edge_vertex2] - y_coord_in[next_vertex];
1134 d3.z = z_coord_in[edge_vertex2] - z_coord_in[next_vertex];
1135
1136 double length1 = vec_length(d1);
1137 double length2 = vec_length(d2);
1138 double length3 = vec_length(d3);
1139
1140 if ((length1 > 0) && (length2 > 0))
1141 {
1142 double cosangle = dot_product(d1, d2) / (length1 * length2);
1143
1144 // Cell is topologically "improper"
1145 // Two edges have the same direction (and share a common vertex)
1146 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1147 {
1148 T_vertices = true;
1149 new_edge_vertex = previous_vertex;
1150 break;
1151 }
1152 }
1153
1154 if ((length1 > 0) && (length3 > 0))
1155 {
1156 double cosangle = dot_product(d1, d3) / (length1 * length3);
1157
1158 // Cell is topologically "improper"
1159 // Two edges have the same direction (and share a common vertex)
1160 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1161 {
1162 T_vertices = true;
1163 new_edge_vertex = next_vertex;
1164 break;
1165 }
1166 }
1167 }
1168 }
1169 }
1170 }
1171 }
1172
1173 /**************************************************/
1174 /* Test if the element has already been processed */
1175 /**************************************************/
1176
1177 if (T_vertices == true)
1178 {
1179 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face2);
1180
1181 /* The current face that contains the vertices has not been processed */
1182 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
1183 {
1184 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
1185 contour.polyhedron_faces.push_back(neighbour_face2);
1186 current_face = neighbour_face2;
1187
1188 /* Update edge vertices */
1189 edge_vertex1 = new_edge_vertex;
1190 data_vertex1 = isodata_in[edge_vertex1];
1191 ring_counter++;
1192 face_flag = 1;
1193 break;
1194 }
1195 }
1196 }
1197 }
1198 }
1199
1200 /*******************************************************************/
1201 /* Case 2: All intersection faces have been processed at least once */
1202 /*******************************************************************/
1203
1204 else
1205 {
1206 for (i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
1207 {
1208 adjacent_vertices = false;
1209 neighbour_face1 = polygon_list[i];
1210 for (j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
1211 {
1212 neighbour_face2 = polygon_list[j];
1213 if (face_flag == 0)
1214 {
1215 /********************************************************/
1216 /* Search among the elements that contain both vertices */
1217 /********************************************************/
1218
1219 if (neighbour_face1 == neighbour_face2)
1220 {
1221 if (neighbour_face1 < num_elem_in - 1)
1222 {
1223 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
1224 {
1225 if (conn_in[k] == edge_vertex1)
1226 {
1227 //vertex_index = k;
1228 if (k == elem_in[neighbour_face1])
1229 {
1230 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
1231 next_vertex = conn_in[k + 1];
1232 }
1233
1234 else if (k < elem_in[neighbour_face1 + 1] - 1)
1235 {
1236 previous_vertex = conn_in[k - 1];
1237 next_vertex = conn_in[k + 1];
1238 }
1239
1240 else if (k == elem_in[neighbour_face1 + 1] - 1)
1241 {
1242 previous_vertex = conn_in[k - 1];
1243 next_vertex = conn_in[elem_in[neighbour_face1]];
1244 }
1245
1246 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1247 {
1248 adjacent_vertices = true;
1249 break;
1250 }
1251 }
1252 }
1253 }
1254
1255 else
1256 {
1257 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
1258 {
1259 if (conn_in[k] == edge_vertex1)
1260 {
1261 //vertex_index = k;
1262 if (k == elem_in[neighbour_face1])
1263 {
1264 previous_vertex = conn_in[num_conn_in - 1];
1265 next_vertex = conn_in[k + 1];
1266 }
1267
1268 else if (k < num_conn_in - 1)
1269 {
1270 previous_vertex = conn_in[k - 1];
1271 next_vertex = conn_in[k + 1];
1272 }
1273
1274 else if (k == num_conn_in - 1)
1275 {
1276 previous_vertex = conn_in[k - 1];
1277 next_vertex = conn_in[elem_in[neighbour_face1]];
1278 }
1279
1280 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1281 {
1282 adjacent_vertices = true;
1283 break;
1284 }
1285 }
1286 }
1287 }
1288
1289 if (adjacent_vertices == true)
1290 {
1291 if (neighbour_face1 != copy_current_face)
1292 {
1293 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
1294 contour.polyhedron_faces.push_back(neighbour_face1);
1295 current_face = neighbour_face1;
1296 ring_counter++;
1297 face_flag = 1;
1298 break;
1299 }
1300 }
1301 }
1302 }
1303 }
1304
1305 if (face_flag == 1)
1306 {
1307 break;
1308 }
1309 }
1310 }
1311
1312 /* A new element to continue tracing the isocontour could not be found */
1313 if (current_face == copy_current_face)
1314 {
1315 abort_tracing_isocontour = true;
1316 }
1317 }
1318 }
1319
1320 else if ((edge_vertex1 < num_coord_in - 1) && (edge_vertex2 == num_coord_in - 1))
1321 {
1322 for (i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
1323 {
1324 adjacent_vertices = false;
1325 neighbour_face1 = polygon_list[i];
1326 for (j = index_list[edge_vertex2]; j < num_conn_in; j++)
1327 {
1328 neighbour_face2 = polygon_list[j];
1329 if (face_flag == 0)
1330 {
1331 /********************************************************/
1332 /* Search among the elements that contain both vertices */
1333 /********************************************************/
1334
1335 if (neighbour_face1 == neighbour_face2)
1336 {
1337 if (neighbour_face1 < num_elem_in - 1)
1338 {
1339 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
1340 {
1341 if (conn_in[k] == edge_vertex1)
1342 {
1343 //vertex_index = k;
1344 if (k == elem_in[neighbour_face1])
1345 {
1346 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
1347 next_vertex = conn_in[k + 1];
1348 }
1349
1350 else if (k < elem_in[neighbour_face1 + 1] - 1)
1351 {
1352 previous_vertex = conn_in[k - 1];
1353 next_vertex = conn_in[k + 1];
1354 }
1355
1356 else if (k == elem_in[neighbour_face1 + 1] - 1)
1357 {
1358 previous_vertex = conn_in[k - 1];
1359 next_vertex = conn_in[elem_in[neighbour_face1]];
1360 }
1361
1362 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1363 {
1364 adjacent_vertices = true;
1365 break;
1366 }
1367 }
1368 }
1369 }
1370
1371 else
1372 {
1373 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
1374 {
1375 if (conn_in[k] == edge_vertex1)
1376 {
1377 //vertex_index = k;
1378 if (k == elem_in[neighbour_face1])
1379 {
1380 previous_vertex = conn_in[num_conn_in - 1];
1381 next_vertex = conn_in[k + 1];
1382 }
1383
1384 else if (k < num_conn_in - 1)
1385 {
1386 previous_vertex = conn_in[k - 1];
1387 next_vertex = conn_in[k + 1];
1388 }
1389
1390 else if (k == num_conn_in - 1)
1391 {
1392 previous_vertex = conn_in[k - 1];
1393 next_vertex = conn_in[elem_in[neighbour_face1]];
1394 }
1395
1396 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1397 {
1398 adjacent_vertices = true;
1399 break;
1400 }
1401 }
1402 }
1403 }
1404
1405 /*************************************************/
1406 /* Test if the element has already been processed */
1407 /*************************************************/
1408
1409 if (adjacent_vertices == true)
1410 {
1411 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
1412
1413 /* The current face that contains the vertices has not been processed */
1414 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
1415 {
1416 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
1417 contour.polyhedron_faces.push_back(neighbour_face1);
1418 current_face = neighbour_face1;
1419 ring_counter++;
1420 face_flag = 1;
1421 break;
1422 }
1423 }
1424 }
1425 }
1426 }
1427
1428 if (face_flag == 1)
1429 {
1430 break;
1431 }
1432 }
1433
1434 /**********************************/
1435 /* No element has been found yet */
1436 /**********************************/
1437
1438 if (face_flag == 0)
1439 {
1440 /**********************************************/
1441 /* Case 1: The cell has an "improper" topology */
1442 /**********************************************/
1443
1444 if (improper_topology)
1445 {
1446 T_vertices = false;
1447
1448 /*************************************************/
1449 /* Search among the elements that share vertex1 */
1450 /*************************************************/
1451
1452 if (!T_vertices)
1453 {
1454 for (i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
1455 {
1456 neighbour_face1 = polygon_list[i];
1457 if (neighbour_face1 != current_face)
1458 {
1459 if (neighbour_face1 < num_elem_in - 1)
1460 {
1461 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
1462 {
1463 if (conn_in[k] == edge_vertex1)
1464 {
1465 //vertex_index = k;
1466 if (k == elem_in[neighbour_face1])
1467 {
1468 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
1469 next_vertex = conn_in[k + 1];
1470 }
1471
1472 else if (k < elem_in[neighbour_face1 + 1] - 1)
1473 {
1474 previous_vertex = conn_in[k - 1];
1475 next_vertex = conn_in[k + 1];
1476 }
1477
1478 else if (k == elem_in[neighbour_face1 + 1] - 1)
1479 {
1480 previous_vertex = conn_in[k - 1];
1481 next_vertex = conn_in[elem_in[neighbour_face1]];
1482 }
1483
1484 /* Check for edges that have the same direction and share a common vertex */
1485 if (previous_vertex != edge_vertex2 && next_vertex != edge_vertex2)
1486 {
1487 // Check direction of edges
1488 EDGE_VECTOR d1;
1489 d1.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
1490 d1.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
1491 d1.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
1492
1493 EDGE_VECTOR d2;
1494 d2.x = x_coord_in[edge_vertex1] - x_coord_in[previous_vertex];
1495 d2.y = y_coord_in[edge_vertex1] - y_coord_in[previous_vertex];
1496 d2.z = z_coord_in[edge_vertex1] - z_coord_in[previous_vertex];
1497
1498 EDGE_VECTOR d3;
1499 d3.x = x_coord_in[edge_vertex1] - x_coord_in[next_vertex];
1500 d3.y = y_coord_in[edge_vertex1] - y_coord_in[next_vertex];
1501 d3.z = z_coord_in[edge_vertex1] - z_coord_in[next_vertex];
1502
1503 double length1 = vec_length(d1);
1504 double length2 = vec_length(d2);
1505 double length3 = vec_length(d3);
1506
1507 if ((length1 > 0) && (length2 > 0))
1508 {
1509 double cosangle = dot_product(d1, d2) / (length1 * length2);
1510
1511 // Cell is topologically "improper"
1512 // Two edges have the same direction (and share a common vertex)
1513 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1514 {
1515 T_vertices = true;
1516 new_edge_vertex = previous_vertex;
1517 break;
1518 }
1519 }
1520
1521 if ((length1 > 0) && (length3 > 0))
1522 {
1523 double cosangle = dot_product(d1, d3) / (length1 * length3);
1524
1525 // Cell is topologically "improper"
1526 // Two edges have the same direction (and share a common vertex)
1527 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1528 {
1529 T_vertices = true;
1530 new_edge_vertex = next_vertex;
1531 break;
1532 }
1533 }
1534 }
1535 }
1536 }
1537 }
1538
1539 else
1540 {
1541 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
1542 {
1543 if (conn_in[k] == edge_vertex1)
1544 {
1545 //vertex_index = k;
1546 if (k == elem_in[neighbour_face1])
1547 {
1548 previous_vertex = conn_in[num_conn_in - 1];
1549 next_vertex = conn_in[k + 1];
1550 }
1551
1552 else if (k < num_conn_in - 1)
1553 {
1554 previous_vertex = conn_in[k - 1];
1555 next_vertex = conn_in[k + 1];
1556 }
1557
1558 else if (k == num_conn_in - 1)
1559 {
1560 previous_vertex = conn_in[k - 1];
1561 next_vertex = conn_in[elem_in[neighbour_face1]];
1562 }
1563
1564 /* Check for edges that have the same direction and share a common vertex */
1565 if (previous_vertex != edge_vertex2 && next_vertex != edge_vertex2)
1566 {
1567 // Check direction of edges
1568 EDGE_VECTOR d1;
1569 d1.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
1570 d1.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
1571 d1.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
1572
1573 EDGE_VECTOR d2;
1574 d2.x = x_coord_in[edge_vertex1] - x_coord_in[previous_vertex];
1575 d2.y = y_coord_in[edge_vertex1] - y_coord_in[previous_vertex];
1576 d2.z = z_coord_in[edge_vertex1] - z_coord_in[previous_vertex];
1577
1578 EDGE_VECTOR d3;
1579 d3.x = x_coord_in[edge_vertex1] - x_coord_in[next_vertex];
1580 d3.y = y_coord_in[edge_vertex1] - y_coord_in[next_vertex];
1581 d3.z = z_coord_in[edge_vertex1] - z_coord_in[next_vertex];
1582
1583 double length1 = vec_length(d1);
1584 double length2 = vec_length(d2);
1585 double length3 = vec_length(d3);
1586
1587 if ((length1 > 0) && (length2 > 0))
1588 {
1589 double cosangle = dot_product(d1, d2) / (length1 * length2);
1590
1591 // Cell is topologically "improper"
1592 // Two edges have the same direction (and share a common vertex)
1593 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1594 {
1595 T_vertices = true;
1596 new_edge_vertex = previous_vertex;
1597 break;
1598 }
1599 }
1600
1601 if ((length1 > 0) && (length3 > 0))
1602 {
1603 double cosangle = dot_product(d1, d3) / (length1 * length3);
1604
1605 // Cell is topologically "improper"
1606 // Two edges have the same direction (and share a common vertex)
1607 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1608 {
1609 T_vertices = true;
1610 new_edge_vertex = next_vertex;
1611 break;
1612 }
1613 }
1614 }
1615 }
1616 }
1617 }
1618 }
1619
1620 if (T_vertices == true)
1621 {
1622 /* Test if the element has already been processed */
1623 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
1624
1625 /* The current face that contains the vertices has not been processed */
1626 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
1627 {
1628 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
1629 contour.polyhedron_faces.push_back(neighbour_face1);
1630 current_face = neighbour_face1;
1631
1632 /* Update edge vertices */
1633 edge_vertex2 = new_edge_vertex;
1634 data_vertex2 = isodata_in[edge_vertex2];
1635 ring_counter++;
1636 face_flag = 1;
1637 break;
1638 }
1639 }
1640 }
1641 }
1642
1643 /*************************************************/
1644 /* Search among the elements that share vertex2 */
1645 /*************************************************/
1646
1647 if (!T_vertices)
1648 {
1649 for (j = index_list[edge_vertex2]; j < num_conn_in; j++)
1650 {
1651 neighbour_face2 = polygon_list[j];
1652 if (neighbour_face2 != current_face)
1653 {
1654 if (neighbour_face2 < num_elem_in - 1)
1655 {
1656 for (k = elem_in[neighbour_face2]; k < elem_in[neighbour_face2 + 1]; k++)
1657 {
1658 if (conn_in[k] == edge_vertex2)
1659 {
1660 //vertex_index = k;
1661 if (k == elem_in[neighbour_face2])
1662 {
1663 previous_vertex = conn_in[elem_in[neighbour_face2 + 1] - 1];
1664 next_vertex = conn_in[k + 1];
1665 }
1666
1667 else if (k < elem_in[neighbour_face2 + 1] - 1)
1668 {
1669 previous_vertex = conn_in[k - 1];
1670 next_vertex = conn_in[k + 1];
1671 }
1672
1673 else if (k == elem_in[neighbour_face2 + 1] - 1)
1674 {
1675 previous_vertex = conn_in[k - 1];
1676 next_vertex = conn_in[elem_in[neighbour_face2]];
1677 }
1678
1679 /* Check for edges that have the same direction and share a common vertex */
1680 if (previous_vertex != edge_vertex1 && next_vertex != edge_vertex1)
1681 {
1682 // Check direction of edges
1683 EDGE_VECTOR d1;
1684 d1.x = x_coord_in[edge_vertex2] - x_coord_in[edge_vertex1];
1685 d1.y = y_coord_in[edge_vertex2] - y_coord_in[edge_vertex1];
1686 d1.z = z_coord_in[edge_vertex2] - z_coord_in[edge_vertex1];
1687
1688 EDGE_VECTOR d2;
1689 d2.x = x_coord_in[edge_vertex2] - x_coord_in[previous_vertex];
1690 d2.y = y_coord_in[edge_vertex2] - y_coord_in[previous_vertex];
1691 d2.z = z_coord_in[edge_vertex2] - z_coord_in[previous_vertex];
1692
1693 EDGE_VECTOR d3;
1694 d3.x = x_coord_in[edge_vertex2] - x_coord_in[next_vertex];
1695 d3.y = y_coord_in[edge_vertex2] - y_coord_in[next_vertex];
1696 d3.z = z_coord_in[edge_vertex2] - z_coord_in[next_vertex];
1697
1698 double length1 = vec_length(d1);
1699 double length2 = vec_length(d2);
1700 double length3 = vec_length(d3);
1701
1702 if ((length1 > 0) && (length2 > 0))
1703 {
1704 double cosangle = dot_product(d1, d2) / (length1 * length2);
1705
1706 // Cell is topologically "improper"
1707 // Two edges have the same direction (and share a common vertex)
1708 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1709 {
1710 T_vertices = true;
1711 new_edge_vertex = previous_vertex;
1712 break;
1713 }
1714 }
1715
1716 if ((length1 > 0) && (length3 > 0))
1717 {
1718 double cosangle = dot_product(d1, d3) / (length1 * length3);
1719
1720 // Cell is topologically "improper"
1721 // Two edges have the same direction (and share a common vertex)
1722 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1723 {
1724 T_vertices = true;
1725 new_edge_vertex = next_vertex;
1726 break;
1727 }
1728 }
1729 }
1730 }
1731 }
1732 }
1733
1734 else
1735 {
1736 for (k = elem_in[neighbour_face2]; k < num_conn_in; k++)
1737 {
1738 if (conn_in[k] == edge_vertex2)
1739 {
1740 //vertex_index = k;
1741 if (k == elem_in[neighbour_face2])
1742 {
1743 previous_vertex = conn_in[num_conn_in - 1];
1744 next_vertex = conn_in[k + 1];
1745 }
1746
1747 else if (k < num_conn_in - 1)
1748 {
1749 previous_vertex = conn_in[k - 1];
1750 next_vertex = conn_in[k + 1];
1751 }
1752
1753 else if (k == num_conn_in - 1)
1754 {
1755 previous_vertex = conn_in[k - 1];
1756 next_vertex = conn_in[elem_in[neighbour_face2]];
1757 }
1758
1759 /* Check for edges that have the same direction and share a common vertex */
1760 if (previous_vertex != edge_vertex1 && next_vertex != edge_vertex1)
1761 {
1762 // Check direction of edges
1763 EDGE_VECTOR d1;
1764 d1.x = x_coord_in[edge_vertex2] - x_coord_in[edge_vertex1];
1765 d1.y = y_coord_in[edge_vertex2] - y_coord_in[edge_vertex1];
1766 d1.z = z_coord_in[edge_vertex2] - z_coord_in[edge_vertex1];
1767
1768 EDGE_VECTOR d2;
1769 d2.x = x_coord_in[edge_vertex2] - x_coord_in[previous_vertex];
1770 d2.y = y_coord_in[edge_vertex2] - y_coord_in[previous_vertex];
1771 d2.z = z_coord_in[edge_vertex2] - z_coord_in[previous_vertex];
1772
1773 EDGE_VECTOR d3;
1774 d3.x = x_coord_in[edge_vertex2] - x_coord_in[next_vertex];
1775 d3.y = y_coord_in[edge_vertex2] - y_coord_in[next_vertex];
1776 d3.z = z_coord_in[edge_vertex2] - z_coord_in[next_vertex];
1777
1778 double length1 = vec_length(d1);
1779 double length2 = vec_length(d2);
1780 double length3 = vec_length(d3);
1781
1782 if ((length1 > 0) && (length2 > 0))
1783 {
1784 double cosangle = dot_product(d1, d2) / (length1 * length2);
1785
1786 // Cell is topologically "improper"
1787 // Two edges have the same direction (and share a common vertex)
1788 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1789 {
1790 T_vertices = true;
1791 new_edge_vertex = previous_vertex;
1792 break;
1793 }
1794 }
1795
1796 if ((length1 > 0) && (length3 > 0))
1797 {
1798 double cosangle = dot_product(d1, d3) / (length1 * length3);
1799
1800 // Cell is topologically "improper"
1801 // Two edges have the same direction (and share a common vertex)
1802 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
1803 {
1804 T_vertices = true;
1805 new_edge_vertex = next_vertex;
1806 break;
1807 }
1808 }
1809 }
1810 }
1811 }
1812 }
1813 }
1814
1815 /*************************************************/
1816 /* Test if the element has already been processed */
1817 /*************************************************/
1818
1819 if (T_vertices == true)
1820 {
1821 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face2);
1822
1823 /* The current face that contains the vertices has not been processed */
1824 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
1825 {
1826 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
1827 contour.polyhedron_faces.push_back(neighbour_face2);
1828 current_face = neighbour_face2;
1829
1830 /* Update edge vertices */
1831 edge_vertex1 = new_edge_vertex;
1832 data_vertex1 = isodata_in[edge_vertex1];
1833 ring_counter++;
1834 face_flag = 1;
1835 break;
1836 }
1837 }
1838 }
1839 }
1840 }
1841
1842 /******************************************************************/
1843 /* Case 2: All intersection faces have been processed at least once */
1844 /******************************************************************/
1845
1846 else
1847 {
1848 for (i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
1849 {
1850 adjacent_vertices = false;
1851 neighbour_face1 = polygon_list[i];
1852 for (j = index_list[edge_vertex2]; j < num_conn_in; j++)
1853 {
1854 neighbour_face2 = polygon_list[j];
1855 if (face_flag == 0)
1856 {
1857 /* Search among the elements that contain both vertices */
1858 if (neighbour_face1 == neighbour_face2)
1859 {
1860 if (neighbour_face1 < num_elem_in - 1)
1861 {
1862 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
1863 {
1864 if (conn_in[k] == edge_vertex1)
1865 {
1866 //vertex_index = k;
1867 if (k == elem_in[neighbour_face1])
1868 {
1869 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
1870 next_vertex = conn_in[k + 1];
1871 }
1872
1873 else if (k < elem_in[neighbour_face1 + 1] - 1)
1874 {
1875 previous_vertex = conn_in[k - 1];
1876 next_vertex = conn_in[k + 1];
1877 }
1878
1879 else if (k == elem_in[neighbour_face1 + 1] - 1)
1880 {
1881 previous_vertex = conn_in[k - 1];
1882 next_vertex = conn_in[elem_in[neighbour_face1]];
1883 }
1884
1885 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1886 {
1887 adjacent_vertices = true;
1888 break;
1889 }
1890 }
1891 }
1892 }
1893
1894 else
1895 {
1896 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
1897 {
1898 if (conn_in[k] == edge_vertex1)
1899 {
1900 //vertex_index = k;
1901 if (k == elem_in[neighbour_face1])
1902 {
1903 previous_vertex = conn_in[num_conn_in - 1];
1904 next_vertex = conn_in[k + 1];
1905 }
1906
1907 else if (k < num_conn_in - 1)
1908 {
1909 previous_vertex = conn_in[k - 1];
1910 next_vertex = conn_in[k + 1];
1911 }
1912
1913 else if (k == num_conn_in - 1)
1914 {
1915 previous_vertex = conn_in[k - 1];
1916 next_vertex = conn_in[elem_in[neighbour_face1]];
1917 }
1918
1919 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1920 {
1921 adjacent_vertices = true;
1922 break;
1923 }
1924 }
1925 }
1926 }
1927
1928 if (adjacent_vertices == true)
1929 {
1930 if (neighbour_face1 != copy_current_face)
1931 {
1932 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
1933 contour.polyhedron_faces.push_back(neighbour_face1);
1934 current_face = neighbour_face1;
1935 ring_counter++;
1936 face_flag = 1;
1937 break;
1938 }
1939 }
1940 }
1941 }
1942 }
1943
1944 if (face_flag == 1)
1945 {
1946 break;
1947 }
1948 }
1949 }
1950
1951 /*A new element to continue tracing the isocontour could not be found */
1952 if (current_face == copy_current_face)
1953 {
1954 abort_tracing_isocontour = true;
1955 }
1956 }
1957 }
1958
1959 else if ((edge_vertex1 == num_coord_in - 1) && (edge_vertex2 < num_coord_in - 1))
1960 {
1961 for (i = index_list[edge_vertex1]; i < num_conn_in; i++)
1962 {
1963 adjacent_vertices = false;
1964 neighbour_face1 = polygon_list[i];
1965 for (j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
1966 {
1967 neighbour_face2 = polygon_list[j];
1968 if (face_flag == 0)
1969 {
1970 /* Search among the elements that contain both vertices */
1971 if (neighbour_face1 == neighbour_face2)
1972 {
1973 if (neighbour_face1 < num_elem_in - 1)
1974 {
1975 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
1976 {
1977 if (conn_in[k] == edge_vertex1)
1978 {
1979 //vertex_index = k;
1980 if (k == elem_in[neighbour_face1])
1981 {
1982 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
1983 next_vertex = conn_in[k + 1];
1984 }
1985
1986 else if (k < elem_in[neighbour_face1 + 1] - 1)
1987 {
1988 previous_vertex = conn_in[k - 1];
1989 next_vertex = conn_in[k + 1];
1990 }
1991
1992 else if (k == elem_in[neighbour_face1 + 1] - 1)
1993 {
1994 previous_vertex = conn_in[k - 1];
1995 next_vertex = conn_in[elem_in[neighbour_face1]];
1996 }
1997
1998 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
1999 {
2000 adjacent_vertices = true;
2001 break;
2002 }
2003 }
2004 }
2005 }
2006
2007 else
2008 {
2009 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
2010 {
2011 if (conn_in[k] == edge_vertex1)
2012 {
2013 //vertex_index = k;
2014 if (k == elem_in[neighbour_face1])
2015 {
2016 previous_vertex = conn_in[num_conn_in - 1];
2017 next_vertex = conn_in[k + 1];
2018 }
2019
2020 else if (k < num_conn_in - 1)
2021 {
2022 previous_vertex = conn_in[k - 1];
2023 next_vertex = conn_in[k + 1];
2024 }
2025
2026 else if (k == num_conn_in - 1)
2027 {
2028 previous_vertex = conn_in[k - 1];
2029 next_vertex = conn_in[elem_in[neighbour_face1]];
2030 }
2031
2032 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
2033 {
2034 adjacent_vertices = true;
2035 break;
2036 }
2037 }
2038 }
2039 }
2040
2041 if (adjacent_vertices == true)
2042 {
2043 /* Test if the element has already been processed */
2044 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
2045
2046 /* The current face that contains the vertices has not been processed */
2047 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
2048 {
2049 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
2050 contour.polyhedron_faces.push_back(neighbour_face1);
2051 current_face = neighbour_face1;
2052 ring_counter++;
2053 face_flag = 1;
2054 break;
2055 }
2056 }
2057 }
2058 }
2059 }
2060
2061 if (face_flag == 1)
2062 {
2063 break;
2064 }
2065 }
2066
2067 /* All intersection faces have been processed at least once or the cell has an "improper" topology */
2068 if (face_flag == 0)
2069 {
2070 if (improper_topology)
2071 {
2072 T_vertices = false;
2073 if (!T_vertices)
2074 {
2075 /* Search among elements that share vertex1 */
2076 for (i = index_list[edge_vertex1]; i < num_conn_in; i++)
2077 {
2078 neighbour_face1 = polygon_list[i];
2079 if (neighbour_face1 != current_face)
2080 {
2081 if (neighbour_face1 < num_elem_in - 1)
2082 {
2083 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
2084 {
2085 if (conn_in[k] == edge_vertex1)
2086 {
2087 //vertex_index = k;
2088 if (k == elem_in[neighbour_face1])
2089 {
2090 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
2091 next_vertex = conn_in[k + 1];
2092 }
2093
2094 else if (k < elem_in[neighbour_face1 + 1] - 1)
2095 {
2096 previous_vertex = conn_in[k - 1];
2097 next_vertex = conn_in[k + 1];
2098 }
2099
2100 else if (k == elem_in[neighbour_face1 + 1] - 1)
2101 {
2102 previous_vertex = conn_in[k - 1];
2103 next_vertex = conn_in[elem_in[neighbour_face1]];
2104 }
2105
2106 /* Check for edges that have the same direction and share a common vertex */
2107 if (previous_vertex != edge_vertex2 && next_vertex != edge_vertex2)
2108 {
2109 // Check direction of edges
2110 EDGE_VECTOR d1;
2111 d1.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
2112 d1.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
2113 d1.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
2114
2115 EDGE_VECTOR d2;
2116 d2.x = x_coord_in[edge_vertex1] - x_coord_in[previous_vertex];
2117 d2.y = y_coord_in[edge_vertex1] - y_coord_in[previous_vertex];
2118 d2.z = z_coord_in[edge_vertex1] - z_coord_in[previous_vertex];
2119
2120 EDGE_VECTOR d3;
2121 d3.x = x_coord_in[edge_vertex1] - x_coord_in[next_vertex];
2122 d3.y = y_coord_in[edge_vertex1] - y_coord_in[next_vertex];
2123 d3.z = z_coord_in[edge_vertex1] - z_coord_in[next_vertex];
2124
2125 double length1 = vec_length(d1);
2126 double length2 = vec_length(d2);
2127 double length3 = vec_length(d3);
2128
2129 if ((length1 > 0) && (length2 > 0))
2130 {
2131 double cosangle = dot_product(d1, d2) / (length1 * length2);
2132
2133 // Cell is topologically "improper"
2134 // Two edges have the same direction (and share a common vertex)
2135 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2136 {
2137 T_vertices = true;
2138 new_edge_vertex = previous_vertex;
2139 break;
2140 }
2141 }
2142
2143 if ((length1 > 0) && (length3 > 0))
2144 {
2145 double cosangle = dot_product(d1, d3) / (length1 * length3);
2146
2147 // Cell is topologically "improper"
2148 // Two edges have the same direction (and share a common vertex)
2149 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2150 {
2151 T_vertices = true;
2152 new_edge_vertex = next_vertex;
2153 break;
2154 }
2155 }
2156 }
2157 }
2158 }
2159 }
2160
2161 else
2162 {
2163 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
2164 {
2165 if (conn_in[k] == edge_vertex1)
2166 {
2167 //vertex_index = k;
2168 if (k == elem_in[neighbour_face1])
2169 {
2170 previous_vertex = conn_in[num_conn_in - 1];
2171 next_vertex = conn_in[k + 1];
2172 }
2173
2174 else if (k < num_conn_in - 1)
2175 {
2176 previous_vertex = conn_in[k - 1];
2177 next_vertex = conn_in[k + 1];
2178 }
2179
2180 else if (k == num_conn_in - 1)
2181 {
2182 previous_vertex = conn_in[k - 1];
2183 next_vertex = conn_in[elem_in[neighbour_face1]];
2184 }
2185
2186 /* Check for edges that have the same direction and share a common vertex */
2187 if (previous_vertex != edge_vertex2 && next_vertex != edge_vertex2)
2188 {
2189 // Check direction of edges
2190 EDGE_VECTOR d1;
2191 d1.x = x_coord_in[edge_vertex1] - x_coord_in[edge_vertex2];
2192 d1.y = y_coord_in[edge_vertex1] - y_coord_in[edge_vertex2];
2193 d1.z = z_coord_in[edge_vertex1] - z_coord_in[edge_vertex2];
2194
2195 EDGE_VECTOR d2;
2196 d2.x = x_coord_in[edge_vertex1] - x_coord_in[previous_vertex];
2197 d2.y = y_coord_in[edge_vertex1] - y_coord_in[previous_vertex];
2198 d2.z = z_coord_in[edge_vertex1] - z_coord_in[previous_vertex];
2199
2200 EDGE_VECTOR d3;
2201 d3.x = x_coord_in[edge_vertex1] - x_coord_in[next_vertex];
2202 d3.y = y_coord_in[edge_vertex1] - y_coord_in[next_vertex];
2203 d3.z = z_coord_in[edge_vertex1] - z_coord_in[next_vertex];
2204
2205 double length1 = vec_length(d1);
2206 double length2 = vec_length(d2);
2207 double length3 = vec_length(d3);
2208
2209 if ((length1 > 0) && (length2 > 0))
2210 {
2211 double cosangle = dot_product(d1, d2) / (length1 * length2);
2212
2213 // Cell is topologically "improper"
2214 // Two edges have the same direction (and share a common vertex)
2215 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2216 {
2217 T_vertices = true;
2218 new_edge_vertex = previous_vertex;
2219 break;
2220 }
2221 }
2222
2223 if ((length1 > 0) && (length3 > 0))
2224 {
2225 double cosangle = dot_product(d1, d3) / (length1 * length3);
2226
2227 // Cell is topologically "improper"
2228 // Two edges have the same direction (and share a common vertex)
2229 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2230 {
2231 T_vertices = true;
2232 new_edge_vertex = next_vertex;
2233 break;
2234 }
2235 }
2236 }
2237 }
2238 }
2239 }
2240 }
2241
2242 if (T_vertices == true)
2243 {
2244 /* Test if the element has already been processed */
2245 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
2246
2247 /* The current face that contains the vertices has not been processed */
2248 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
2249 {
2250 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
2251 contour.polyhedron_faces.push_back(neighbour_face1);
2252 current_face = neighbour_face1;
2253
2254 /* Update edge vertices */
2255 edge_vertex2 = new_edge_vertex;
2256 data_vertex2 = isodata_in[edge_vertex2];
2257 ring_counter++;
2258 face_flag = 1;
2259 break;
2260 }
2261 }
2262 }
2263 }
2264
2265 if (!T_vertices)
2266 {
2267 /* Search among elements that share vertex2 */
2268 for (j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
2269 {
2270 neighbour_face2 = polygon_list[j];
2271 if (neighbour_face2 != current_face)
2272 {
2273 if (neighbour_face2 < num_elem_in - 1)
2274 {
2275 for (k = elem_in[neighbour_face2]; k < elem_in[neighbour_face2 + 1]; k++)
2276 {
2277 if (conn_in[k] == edge_vertex2)
2278 {
2279 //vertex_index = k;
2280 if (k == elem_in[neighbour_face2])
2281 {
2282 previous_vertex = conn_in[elem_in[neighbour_face2 + 1] - 1];
2283 next_vertex = conn_in[k + 1];
2284 }
2285
2286 else if (k < elem_in[neighbour_face2 + 1] - 1)
2287 {
2288 previous_vertex = conn_in[k - 1];
2289 next_vertex = conn_in[k + 1];
2290 }
2291
2292 else if (k == elem_in[neighbour_face2 + 1] - 1)
2293 {
2294 previous_vertex = conn_in[k - 1];
2295 next_vertex = conn_in[elem_in[neighbour_face2]];
2296 }
2297
2298 /* Check for edges that have the same direction and share a common vertex */
2299 if (previous_vertex != edge_vertex1 && next_vertex != edge_vertex1)
2300 {
2301 // Check direction of edges
2302 EDGE_VECTOR d1;
2303 d1.x = x_coord_in[edge_vertex2] - x_coord_in[edge_vertex1];
2304 d1.y = y_coord_in[edge_vertex2] - y_coord_in[edge_vertex1];
2305 d1.z = z_coord_in[edge_vertex2] - z_coord_in[edge_vertex1];
2306
2307 EDGE_VECTOR d2;
2308 d2.x = x_coord_in[edge_vertex2] - x_coord_in[previous_vertex];
2309 d2.y = y_coord_in[edge_vertex2] - y_coord_in[previous_vertex];
2310 d2.z = z_coord_in[edge_vertex2] - z_coord_in[previous_vertex];
2311
2312 EDGE_VECTOR d3;
2313 d3.x = x_coord_in[edge_vertex2] - x_coord_in[next_vertex];
2314 d3.y = y_coord_in[edge_vertex2] - y_coord_in[next_vertex];
2315 d3.z = z_coord_in[edge_vertex2] - z_coord_in[next_vertex];
2316
2317 double length1 = vec_length(d1);
2318 double length2 = vec_length(d2);
2319 double length3 = vec_length(d3);
2320
2321 if ((length1 > 0) && (length2 > 0))
2322 {
2323 double cosangle = dot_product(d1, d2) / (length1 * length2);
2324
2325 // Cell is topologically "improper"
2326 // Two edges have the same direction (and share a common vertex)
2327 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2328 {
2329 T_vertices = true;
2330 new_edge_vertex = previous_vertex;
2331 break;
2332 }
2333 }
2334
2335 if ((length1 > 0) && (length3 > 0))
2336 {
2337 double cosangle = dot_product(d1, d3) / (length1 * length3);
2338
2339 // Cell is topologically "improper"
2340 // Two edges have the same direction (and share a common vertex)
2341 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2342 {
2343 T_vertices = true;
2344 new_edge_vertex = next_vertex;
2345 break;
2346 }
2347 }
2348 }
2349 }
2350 }
2351 }
2352
2353 else
2354 {
2355 for (k = elem_in[neighbour_face2]; k < num_conn_in; k++)
2356 {
2357 if (conn_in[k] == edge_vertex2)
2358 {
2359 //vertex_index = k;
2360 if (k == elem_in[neighbour_face2])
2361 {
2362 previous_vertex = conn_in[num_conn_in - 1];
2363 next_vertex = conn_in[k + 1];
2364 }
2365
2366 else if (k < num_conn_in - 1)
2367 {
2368 previous_vertex = conn_in[k - 1];
2369 next_vertex = conn_in[k + 1];
2370 }
2371
2372 else if (k == num_conn_in - 1)
2373 {
2374 previous_vertex = conn_in[k - 1];
2375 next_vertex = conn_in[elem_in[neighbour_face2]];
2376 }
2377
2378 /* Check for edges that have the same direction and share a common vertex */
2379 if (previous_vertex != edge_vertex1 && next_vertex != edge_vertex1)
2380 {
2381 // Check direction of edges
2382 EDGE_VECTOR d1;
2383 d1.x = x_coord_in[edge_vertex2] - x_coord_in[edge_vertex1];
2384 d1.y = y_coord_in[edge_vertex2] - y_coord_in[edge_vertex1];
2385 d1.z = z_coord_in[edge_vertex2] - z_coord_in[edge_vertex1];
2386
2387 EDGE_VECTOR d2;
2388 d2.x = x_coord_in[edge_vertex2] - x_coord_in[previous_vertex];
2389 d2.y = y_coord_in[edge_vertex2] - y_coord_in[previous_vertex];
2390 d2.z = z_coord_in[edge_vertex2] - z_coord_in[previous_vertex];
2391
2392 EDGE_VECTOR d3;
2393 d3.x = x_coord_in[edge_vertex2] - x_coord_in[next_vertex];
2394 d3.y = y_coord_in[edge_vertex2] - y_coord_in[next_vertex];
2395 d3.z = z_coord_in[edge_vertex2] - z_coord_in[next_vertex];
2396
2397 double length1 = vec_length(d1);
2398 double length2 = vec_length(d2);
2399 double length3 = vec_length(d3);
2400
2401 if ((length1 > 0) && (length2 > 0))
2402 {
2403 double cosangle = dot_product(d1, d2) / (length1 * length2);
2404
2405 // Cell is topologically "improper"
2406 // Two edges have the same direction (and share a common vertex)
2407 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2408 {
2409 T_vertices = true;
2410 new_edge_vertex = previous_vertex;
2411 break;
2412 }
2413 }
2414
2415 if ((length1 > 0) && (length3 > 0))
2416 {
2417 double cosangle = dot_product(d1, d3) / (length1 * length3);
2418
2419 // Cell is topologically "improper"
2420 // Two edges have the same direction (and share a common vertex)
2421 if (fabs(cosangle - 1.0) < 0.00001 && cosangle > 0)
2422 {
2423 T_vertices = true;
2424 new_edge_vertex = next_vertex;
2425 break;
2426 }
2427 }
2428 }
2429 }
2430 }
2431 }
2432 }
2433
2434 if (T_vertices == true)
2435 {
2436 /* Test if the element has already been processed */
2437 it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face2);
2438
2439 /* The current face that contains the vertices has not been processed */
2440 if (it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
2441 {
2442 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
2443 contour.polyhedron_faces.push_back(neighbour_face2);
2444 current_face = neighbour_face2;
2445
2446 /* Update edge vertices */
2447 edge_vertex1 = new_edge_vertex;
2448 data_vertex1 = isodata_in[edge_vertex1];
2449 ring_counter++;
2450 face_flag = 1;
2451 break;
2452 }
2453 }
2454 }
2455 }
2456 }
2457
2458 else
2459 {
2460 for (i = index_list[edge_vertex1]; i < num_conn_in; i++)
2461 {
2462 adjacent_vertices = false;
2463 neighbour_face1 = polygon_list[i];
2464 for (j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
2465 {
2466 neighbour_face2 = polygon_list[j];
2467 if (face_flag == 0)
2468 {
2469 /* Search among the elements that contain both vertices */
2470 if (neighbour_face1 == neighbour_face2)
2471 {
2472 if (neighbour_face1 < num_elem_in - 1)
2473 {
2474 for (k = elem_in[neighbour_face1]; k < elem_in[neighbour_face1 + 1]; k++)
2475 {
2476 if (conn_in[k] == edge_vertex1)
2477 {
2478 //vertex_index = k;
2479 if (k == elem_in[neighbour_face1])
2480 {
2481 previous_vertex = conn_in[elem_in[neighbour_face1 + 1] - 1];
2482 next_vertex = conn_in[k + 1];
2483 }
2484
2485 else if (k < elem_in[neighbour_face1 + 1] - 1)
2486 {
2487 previous_vertex = conn_in[k - 1];
2488 next_vertex = conn_in[k + 1];
2489 }
2490
2491 else if (k == elem_in[neighbour_face1 + 1] - 1)
2492 {
2493 previous_vertex = conn_in[k - 1];
2494 next_vertex = conn_in[elem_in[neighbour_face1]];
2495 }
2496
2497 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
2498 {
2499 adjacent_vertices = true;
2500 break;
2501 }
2502 }
2503 }
2504 }
2505
2506 else
2507 {
2508 for (k = elem_in[neighbour_face1]; k < num_conn_in; k++)
2509 {
2510 if (conn_in[k] == edge_vertex1)
2511 {
2512 //vertex_index = k;
2513 if (k == elem_in[neighbour_face1])
2514 {
2515 previous_vertex = conn_in[num_conn_in - 1];
2516 next_vertex = conn_in[k + 1];
2517 }
2518
2519 else if (k < num_conn_in - 1)
2520 {
2521 previous_vertex = conn_in[k - 1];
2522 next_vertex = conn_in[k + 1];
2523 }
2524
2525 else if (k == num_conn_in - 1)
2526 {
2527 previous_vertex = conn_in[k - 1];
2528 next_vertex = conn_in[elem_in[neighbour_face1]];
2529 }
2530
2531 if (previous_vertex == edge_vertex2 || next_vertex == edge_vertex2)
2532 {
2533 adjacent_vertices = true;
2534 break;
2535 }
2536 }
2537 }
2538 }
2539
2540 if (adjacent_vertices == true)
2541 {
2542 if (neighbour_face1 != copy_current_face)
2543 {
2544 contour.ring.push_back(assign_int_index(intsec_vector, edge_vertex1, edge_vertex2));
2545 contour.polyhedron_faces.push_back(neighbour_face1);
2546 current_face = neighbour_face1;
2547 ring_counter++;
2548 face_flag = 1;
2549 break;
2550 }
2551 }
2552 }
2553 }
2554 }
2555
2556 if (face_flag == 1)
2557 {
2558 break;
2559 }
2560 }
2561 }
2562
2563 /* A new element to continue tracing the isocontour could not be found */
2564 if (current_face == copy_current_face)
2565 {
2566 abort_tracing_isocontour = true;
2567 }
2568 }
2569 }
2570}
2571
2572void generate_isocontour(ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, float data_vertex1, float data_vertex2, int edge_vertex1, int edge_vertex2, int &new_edge_vertex1, int &new_edge_vertex2, int *elem_in, int *conn_in, float isovalue, int num_elem_in, int num_conn_in, int current_face, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool improper_topology, bool &abort_tracing_isocontour, CONTOUR contour, int num_of_rings, int ring_end)
2573{
2574 //bool new_int_found;
2575
2576 int i;
2577 int int_index;
2578 int index_v1;
2579 int index_v2;
2580 int index_flag_v1;
2581 int index_flag_v2;
2582 int counter1;
2583 int counter2;
2584
2585 ITERATOR it;
2586 ITERATOR it2;
2587
2588 /************************************/
2589 /* Test configuration of the vertices */
2590 /************************************/
2591
2592 //new_int_found = false;
2593 index_flag_v1 = 0;
2594 index_flag_v2 = 0;
2595 index_v1 = elem_in[current_face];
2596 index_v2 = elem_in[current_face];
2597
2598 /***********************************************************************/
2599 /* Locate index values of the edge vertices in the array of connectivities */
2600 /***********************************************************************/
2601
2602 if (current_face < num_elem_in - 1)
2603 {
2604 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
2605 {
2606 if (index_flag_v1 == 0)
2607 {
2608 if (edge_vertex1 != conn_in[index_v1])
2609 {
2610 index_v1++;
2611 }
2612
2613 else
2614 {
2615 index_flag_v1 = 1;
2616 }
2617 }
2618
2619 if (index_flag_v2 == 0)
2620 {
2621 if (edge_vertex2 != conn_in[index_v2])
2622 {
2623 index_v2++;
2624 }
2625
2626 else
2627 {
2628 index_flag_v2 = 1;
2629 }
2630 }
2631 }
2632
2633 /* This should not happen */
2634 if (index_v1 == elem_in[current_face + 1] || index_v2 == elem_in[current_face + 1])
2635 {
2636 abort_tracing_isocontour = true;
2637 }
2638 }
2639
2640 else
2641 {
2642 for (i = elem_in[current_face]; i < num_conn_in; i++)
2643 {
2644 if (index_flag_v1 == 0)
2645 {
2646 if (edge_vertex1 != conn_in[index_v1])
2647 {
2648 index_v1++;
2649 }
2650
2651 else
2652 {
2653 index_flag_v1 = 1;
2654 }
2655 }
2656
2657 if (index_flag_v2 == 0)
2658 {
2659 if (edge_vertex2 != conn_in[index_v2])
2660 {
2661 index_v2++;
2662 }
2663
2664 else
2665 {
2666 index_flag_v2 = 1;
2667 }
2668 }
2669 }
2670
2671 /* This should not happen */
2672 if (index_v1 == num_conn_in || index_v2 == num_conn_in)
2673 {
2674 abort_tracing_isocontour = true;
2675 }
2676 }
2677
2678 counter1 = index_v1;
2679 counter2 = index_v2;
2680
2681 if (!abort_tracing_isocontour)
2682 {
2683 /*******************************/
2684 /* Determine Tracing Direction */
2685 /*******************************/
2686
2687 /*********************************************************/
2688 /* Configuration 1: Data vertex 1 (+) --- Data vertex 2 (-) */
2689 /*********************************************************/
2690
2691 if (((data_vertex1 > isovalue) && (isovalue >= data_vertex2)) || ((data_vertex1 >= isovalue) && (isovalue > data_vertex2)))
2692 {
2693 if (current_face < num_elem_in - 1)
2694 {
2695 /***************************/
2696 /* Case 1: Clockwise (CW) */
2697 /***************************/
2698
2699 if (index_v1 > index_v2)
2700 {
2701 if ((index_v1 != elem_in[current_face + 1] - 1) || (index_v2 != elem_in[current_face]))
2702 {
2703 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
2704 {
2705 counter1++;
2706 counter2++;
2707
2708 if (counter1 == elem_in[current_face + 1])
2709 {
2710 counter1 = elem_in[current_face];
2711 }
2712
2713 if (counter2 == elem_in[current_face + 1])
2714 {
2715 counter2 = elem_in[current_face];
2716 }
2717
2718 /* Test if the a new intersection has been encountered */
2719 new_edge_vertex1 = conn_in[counter1];
2720 new_edge_vertex2 = conn_in[counter2];
2721
2722 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
2723 {
2724 if (num_of_rings == 0)
2725 {
2726 it = find(contour.ring.begin(), contour.ring.end(), int_index);
2727 if (it == contour.ring.end() || it == contour.ring.begin())
2728 {
2729 break;
2730 }
2731 }
2732
2733 else
2734 {
2735 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
2736 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
2737 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
2738 {
2739 break;
2740 }
2741 }
2742 }
2743 }
2744 }
2745
2746 else if ((index_v1 == elem_in[current_face + 1] - 1) && (index_v2 == elem_in[current_face]))
2747 {
2748 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
2749 {
2750 counter1--;
2751 counter2--;
2752
2753 if (counter1 == elem_in[current_face] - 1)
2754 {
2755 counter1 = elem_in[current_face + 1] - 1;
2756 }
2757
2758 if (counter2 == elem_in[current_face] - 1)
2759 {
2760 counter2 = elem_in[current_face + 1] - 1;
2761 }
2762
2763 /* Test if the a new intersection has been encountered */
2764 new_edge_vertex1 = conn_in[counter1];
2765 new_edge_vertex2 = conn_in[counter2];
2766
2767 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
2768 {
2769 if (num_of_rings == 0)
2770 {
2771 it = find(contour.ring.begin(), contour.ring.end(), int_index);
2772 if (it == contour.ring.end() || it == contour.ring.begin())
2773 {
2774 break;
2775 }
2776 }
2777
2778 else
2779 {
2780 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
2781 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
2782 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
2783 {
2784 break;
2785 }
2786 }
2787 }
2788 }
2789 }
2790 }
2791
2792 /************************************/
2793 /* Case 2: Counterclockwise (CCW) */
2794 /************************************/
2795
2796 else if (index_v1 < index_v2)
2797 {
2798 if ((index_v1 != elem_in[current_face]) || (index_v2 != elem_in[current_face + 1] - 1))
2799 {
2800 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
2801 {
2802 counter1--;
2803 counter2--;
2804
2805 if (counter1 == elem_in[current_face] - 1)
2806 {
2807 counter1 = elem_in[current_face + 1] - 1;
2808 }
2809
2810 if (counter2 == elem_in[current_face] - 1)
2811 {
2812 counter2 = elem_in[current_face + 1] - 1;
2813 }
2814
2815 /* Test if the a new intersection has been encountered */
2816 new_edge_vertex1 = conn_in[counter1];
2817 new_edge_vertex2 = conn_in[counter2];
2818
2819 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
2820 {
2821 if (num_of_rings == 0)
2822 {
2823 it = find(contour.ring.begin(), contour.ring.end(), int_index);
2824 if (it == contour.ring.end() || it == contour.ring.begin())
2825 {
2826 break;
2827 }
2828 }
2829
2830 else
2831 {
2832 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
2833 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
2834 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
2835 {
2836 break;
2837 }
2838 }
2839 }
2840 }
2841 }
2842
2843 else if ((index_v1 == elem_in[current_face]) && (index_v2 == elem_in[current_face + 1] - 1))
2844 {
2845 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
2846 {
2847 counter1++;
2848 counter2++;
2849
2850 if (counter1 == elem_in[current_face + 1])
2851 {
2852 counter1 = elem_in[current_face];
2853 }
2854
2855 if (counter2 == elem_in[current_face + 1])
2856 {
2857 counter2 = elem_in[current_face];
2858 }
2859
2860 /* Test if the a new intersection has been encountered */
2861 new_edge_vertex1 = conn_in[counter1];
2862 new_edge_vertex2 = conn_in[counter2];
2863
2864 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
2865 {
2866 if (num_of_rings == 0)
2867 {
2868 it = find(contour.ring.begin(), contour.ring.end(), int_index);
2869 if (it == contour.ring.end() || it == contour.ring.begin())
2870 {
2871 break;
2872 }
2873 }
2874
2875 else
2876 {
2877 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
2878 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
2879 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
2880 {
2881 break;
2882 }
2883 }
2884 }
2885 }
2886 }
2887 }
2888 }
2889
2890 else
2891 {
2892 /***************************/
2893 /* Case 1: Clockwise (CW) */
2894 /***************************/
2895
2896 if (index_v1 > index_v2)
2897 {
2898 if ((index_v1 != num_conn_in - 1) || (index_v2 != elem_in[current_face]))
2899 {
2900 for (i = elem_in[current_face]; i < num_conn_in; i++)
2901 {
2902 counter1++;
2903 counter2++;
2904
2905 if (counter1 == num_conn_in)
2906 {
2907 counter1 = elem_in[current_face];
2908 }
2909
2910 if (counter2 == num_conn_in)
2911 {
2912 counter2 = elem_in[current_face];
2913 }
2914
2915 /* Test if the a new intersection has been encountered */
2916 new_edge_vertex1 = conn_in[counter1];
2917 new_edge_vertex2 = conn_in[counter2];
2918
2919 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
2920 {
2921 if (num_of_rings == 0)
2922 {
2923 it = find(contour.ring.begin(), contour.ring.end(), int_index);
2924 if (it == contour.ring.end() || it == contour.ring.begin())
2925 {
2926 break;
2927 }
2928 }
2929
2930 else
2931 {
2932 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
2933 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
2934 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
2935 {
2936 break;
2937 }
2938 }
2939 }
2940 }
2941 }
2942
2943 else if ((index_v1 == num_conn_in - 1) && (index_v2 == elem_in[current_face]))
2944 {
2945 for (i = elem_in[current_face]; i < num_conn_in; i++)
2946 {
2947 counter1--;
2948 counter2--;
2949
2950 if (counter1 == elem_in[current_face] - 1)
2951 {
2952 counter1 = num_conn_in - 1;
2953 }
2954
2955 if (counter2 == elem_in[current_face] - 1)
2956 {
2957 counter2 = num_conn_in - 1;
2958 }
2959
2960 /* Test if the a new intersection has been encountered */
2961 new_edge_vertex1 = conn_in[counter1];
2962 new_edge_vertex2 = conn_in[counter2];
2963
2964 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
2965 {
2966 if (num_of_rings == 0)
2967 {
2968 it = find(contour.ring.begin(), contour.ring.end(), int_index);
2969 if (it == contour.ring.end() || it == contour.ring.begin())
2970 {
2971 break;
2972 }
2973 }
2974
2975 else
2976 {
2977 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
2978 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
2979 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
2980 {
2981 break;
2982 }
2983 }
2984 }
2985 }
2986 }
2987 }
2988
2989 /************************************/
2990 /* Case 2: Counterclockwise (CCW) */
2991 /************************************/
2992
2993 else if (index_v1 < index_v2)
2994 {
2995 if ((index_v1 != elem_in[current_face]) || (index_v2 != num_conn_in - 1))
2996 {
2997 for (i = elem_in[current_face]; i < num_conn_in; i++)
2998 {
2999 counter1--;
3000 counter2--;
3001
3002 if (counter1 == elem_in[current_face] - 1)
3003 {
3004 counter1 = num_conn_in - 1;
3005 }
3006
3007 if (counter2 == elem_in[current_face] - 1)
3008 {
3009 counter2 = num_conn_in - 1;
3010 }
3011
3012 /* Test if the a new intersection has been encountered */
3013 new_edge_vertex1 = conn_in[counter1];
3014 new_edge_vertex2 = conn_in[counter2];
3015
3016 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3017 {
3018 if (num_of_rings == 0)
3019 {
3020 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3021 if (it == contour.ring.end() || it == contour.ring.begin())
3022 {
3023 break;
3024 }
3025 }
3026
3027 else
3028 {
3029 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3030 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3031 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3032 {
3033 break;
3034 }
3035 }
3036 }
3037 }
3038 }
3039
3040 else if ((index_v1 == elem_in[current_face]) && (index_v2 == num_conn_in - 1))
3041 {
3042 for (i = elem_in[current_face]; i < num_conn_in; i++)
3043 {
3044 counter1++;
3045 counter2++;
3046
3047 if (counter1 == num_conn_in)
3048 {
3049 counter1 = elem_in[current_face];
3050 }
3051
3052 if (counter2 == num_conn_in)
3053 {
3054 counter2 = elem_in[current_face];
3055 }
3056
3057 /* Test if the a new intersection has been encountered */
3058 new_edge_vertex1 = conn_in[counter1];
3059 new_edge_vertex2 = conn_in[counter2];
3060
3061 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3062 {
3063 if (num_of_rings == 0)
3064 {
3065 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3066 if (it == contour.ring.end() || it == contour.ring.begin())
3067 {
3068 break;
3069 }
3070 }
3071
3072 else
3073 {
3074 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3075 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3076 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3077 {
3078 break;
3079 }
3080 }
3081 }
3082 }
3083 }
3084 }
3085 }
3086
3087 if ((edge_vertex1 == new_edge_vertex1 && edge_vertex2 == new_edge_vertex2) || (edge_vertex1 == new_edge_vertex2 && edge_vertex2 == new_edge_vertex1))
3088 {
3089 abort_tracing_isocontour = true;
3090 }
3091 }
3092
3093 /*********************************************************/
3094 /* Configuration 2: Data vertex 1 (-) --- Data vertex 2 (+) */
3095 /*********************************************************/
3096
3097 else if (((data_vertex2 > isovalue) && (isovalue >= data_vertex1)) || ((data_vertex2 >= isovalue) && (isovalue > data_vertex1)))
3098 {
3099 if (current_face < num_elem_in - 1)
3100 {
3101 /***************************/
3102 /* Case 1: Clockwise (CW) */
3103 /***************************/
3104
3105 if (index_v1 < index_v2)
3106 {
3107 if ((index_v1 != elem_in[current_face]) || (index_v2 != elem_in[current_face + 1] - 1))
3108 {
3109 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
3110 {
3111 counter1++;
3112 counter2++;
3113
3114 if (counter1 == elem_in[current_face + 1])
3115 {
3116 counter1 = elem_in[current_face];
3117 }
3118
3119 if (counter2 == elem_in[current_face + 1])
3120 {
3121 counter2 = elem_in[current_face];
3122 }
3123
3124 /* Test if the a new intersection has been encountered */
3125 new_edge_vertex1 = conn_in[counter1];
3126 new_edge_vertex2 = conn_in[counter2];
3127
3128 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3129 {
3130 if (num_of_rings == 0)
3131 {
3132 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3133 if (it == contour.ring.end() || it == contour.ring.begin())
3134 {
3135 break;
3136 }
3137 }
3138
3139 else
3140 {
3141 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3142 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3143 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3144 {
3145 break;
3146 }
3147 }
3148 }
3149 }
3150 }
3151
3152 else if ((index_v1 == elem_in[current_face]) && (index_v2 == elem_in[current_face + 1] - 1))
3153 {
3154 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
3155 {
3156 counter1--;
3157 counter2--;
3158
3159 if (counter1 == elem_in[current_face] - 1)
3160 {
3161 counter1 = elem_in[current_face + 1] - 1;
3162 }
3163
3164 if (counter2 == elem_in[current_face] - 1)
3165 {
3166 counter2 = elem_in[current_face + 1] - 1;
3167 }
3168
3169 /* Test if the a new intersection has been encountered */
3170 new_edge_vertex1 = conn_in[counter1];
3171 new_edge_vertex2 = conn_in[counter2];
3172
3173 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3174 {
3175 if (num_of_rings == 0)
3176 {
3177 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3178 if (it == contour.ring.end() || it == contour.ring.begin())
3179 {
3180 break;
3181 }
3182 }
3183
3184 else
3185 {
3186 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3187 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3188 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3189 {
3190 break;
3191 }
3192 }
3193 }
3194 }
3195 }
3196 }
3197
3198 /************************************/
3199 /* Case 2: Counterclockwise (CCW) */
3200 /************************************/
3201
3202 else if (index_v1 > index_v2)
3203 {
3204 if ((index_v1 != elem_in[current_face + 1] - 1) || (index_v2 != elem_in[current_face]))
3205 {
3206 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
3207 {
3208 counter1--;
3209 counter2--;
3210
3211 if (counter1 == elem_in[current_face] - 1)
3212 {
3213 counter1 = elem_in[current_face + 1] - 1;
3214 }
3215
3216 if (counter2 == elem_in[current_face] - 1)
3217 {
3218 counter2 = elem_in[current_face + 1] - 1;
3219 }
3220
3221 /* Test if the a new intersection has been encountered */
3222 new_edge_vertex1 = conn_in[counter1];
3223 new_edge_vertex2 = conn_in[counter2];
3224
3225 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3226 {
3227 if (num_of_rings == 0)
3228 {
3229 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3230 if (it == contour.ring.end() || it == contour.ring.begin())
3231 {
3232 break;
3233 }
3234 }
3235
3236 else
3237 {
3238 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3239 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3240 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3241 {
3242 break;
3243 }
3244 }
3245 }
3246 }
3247 }
3248
3249 else if ((index_v1 == elem_in[current_face + 1] - 1) && (index_v2 == elem_in[current_face]))
3250 {
3251 for (i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
3252 {
3253 counter1++;
3254 counter2++;
3255
3256 if (counter1 == elem_in[current_face + 1])
3257 {
3258 counter1 = elem_in[current_face];
3259 }
3260
3261 if (counter2 == elem_in[current_face + 1])
3262 {
3263 counter2 = elem_in[current_face];
3264 }
3265
3266 /* Test if the a new intersection has been encountered */
3267 new_edge_vertex1 = conn_in[counter1];
3268 new_edge_vertex2 = conn_in[counter2];
3269
3270 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3271 {
3272 if (num_of_rings == 0)
3273 {
3274 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3275 if (it == contour.ring.end() || it == contour.ring.begin())
3276 {
3277 break;
3278 }
3279 }
3280
3281 else
3282 {
3283 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3284 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3285 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3286 {
3287 break;
3288 }
3289 }
3290 }
3291 }
3292 }
3293 }
3294 }
3295
3296 else
3297 {
3298 /***************************/
3299 /* Case 1: Clockwise (CW) */
3300 /***************************/
3301
3302 if (index_v1 < index_v2)
3303 {
3304 if ((index_v1 != elem_in[current_face]) || (index_v2 != num_conn_in - 1))
3305 {
3306 for (i = elem_in[current_face]; i < num_conn_in; i++)
3307 {
3308 counter1++;
3309 counter2++;
3310
3311 if (counter1 == num_conn_in)
3312 {
3313 counter1 = elem_in[current_face];
3314 }
3315
3316 if (counter2 == num_conn_in)
3317 {
3318 counter2 = elem_in[current_face];
3319 }
3320
3321 /* Test if the a new intersection has been encountered */
3322 new_edge_vertex1 = conn_in[counter1];
3323 new_edge_vertex2 = conn_in[counter2];
3324
3325 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3326 {
3327 if (num_of_rings == 0)
3328 {
3329 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3330 if (it == contour.ring.end() || it == contour.ring.begin())
3331 {
3332 break;
3333 }
3334 }
3335
3336 else
3337 {
3338 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3339 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3340 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3341 {
3342 break;
3343 }
3344 }
3345 }
3346 }
3347 }
3348
3349 else if ((index_v1 == elem_in[current_face]) && (index_v2 == num_conn_in - 1))
3350 {
3351 for (i = elem_in[current_face]; i < num_conn_in; i++)
3352 {
3353 counter1--;
3354 counter2--;
3355
3356 if (counter1 == elem_in[current_face] - 1)
3357 {
3358 counter1 = num_conn_in - 1;
3359 }
3360
3361 if (counter2 == elem_in[current_face] - 1)
3362 {
3363 counter2 = num_conn_in - 1;
3364 }
3365
3366 /* Test if the a new intersection has been encountered */
3367 new_edge_vertex1 = conn_in[counter1];
3368 new_edge_vertex2 = conn_in[counter2];
3369
3370 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3371 {
3372 if (num_of_rings == 0)
3373 {
3374 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3375 if (it == contour.ring.end() || it == contour.ring.begin())
3376 {
3377 break;
3378 }
3379 }
3380
3381 else
3382 {
3383 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3384 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3385 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3386 {
3387 break;
3388 }
3389 }
3390 }
3391 }
3392 }
3393 }
3394
3395 /************************************/
3396 /* Case 2: Counterclockwise (CCW) */
3397 /************************************/
3398
3399 else if (index_v1 > index_v2)
3400 {
3401 if ((index_v1 != num_conn_in - 1) || (index_v2 != elem_in[current_face]))
3402 {
3403 for (i = elem_in[current_face]; i < num_conn_in; i++)
3404 {
3405 counter1--;
3406 counter2--;
3407
3408 if (counter1 == elem_in[current_face] - 1)
3409 {
3410 counter1 = num_conn_in - 1;
3411 }
3412
3413 if (counter2 == elem_in[current_face] - 1)
3414 {
3415 counter2 = num_conn_in - 1;
3416 }
3417
3418 /* Test if the a new intersection has been encountered */
3419 new_edge_vertex1 = conn_in[counter1];
3420 new_edge_vertex2 = conn_in[counter2];
3421
3422 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3423 {
3424 if (num_of_rings == 0)
3425 {
3426 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3427 if (it == contour.ring.end() || it == contour.ring.begin())
3428 {
3429 break;
3430 }
3431 }
3432
3433 else
3434 {
3435 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3436 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3437 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3438 {
3439 break;
3440 }
3441 }
3442 }
3443 }
3444 }
3445
3446 else if ((index_v1 == num_conn_in - 1) && (index_v2 == elem_in[current_face]))
3447 {
3448 for (i = elem_in[current_face]; i < num_conn_in; i++)
3449 {
3450 counter1++;
3451 counter2++;
3452
3453 if (counter1 == num_conn_in)
3454 {
3455 counter1 = elem_in[current_face];
3456 }
3457
3458 if (counter2 == num_conn_in)
3459 {
3460 counter2 = elem_in[current_face];
3461 }
3462
3463 /* Test if the a new intersection has been encountered */
3464 new_edge_vertex1 = conn_in[counter1];
3465 new_edge_vertex2 = conn_in[counter2];
3466
3467 if (find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2, x_coord_in, y_coord_in, z_coord_in, improper_topology, int_index))
3468 {
3469 if (num_of_rings == 0)
3470 {
3471 it = find(contour.ring.begin(), contour.ring.end(), int_index);
3472 if (it == contour.ring.end() || it == contour.ring.begin())
3473 {
3474 break;
3475 }
3476 }
3477
3478 else
3479 {
3480 it = find(contour.ring.begin(), contour.ring.begin() + ring_end, int_index);
3481 it2 = find(contour.ring.begin() + ring_end, contour.ring.end(), int_index);
3482 if (it == contour.ring.begin() + ring_end && (it2 == contour.ring.end() || it2 == contour.ring.begin() + ring_end))
3483 {
3484 break;
3485 }
3486 }
3487 }
3488 }
3489 }
3490 }
3491 }
3492
3493 if ((edge_vertex1 == new_edge_vertex1 && edge_vertex2 == new_edge_vertex2) || (edge_vertex1 == new_edge_vertex2 && edge_vertex2 == new_edge_vertex1))
3494 {
3495 abort_tracing_isocontour = true;
3496 }
3497 }
3498
3499 /***********************************************************************************************************************/
3500 /* Configuration 3: Degenerate case --> Data vertex 1 (+) --- Data vertex 2 (+) or Data vertex 1 (-) --- Data vertex 2 (-) */
3501 /* */
3502 /* This condition occurs only in topologically "improper" polyhedral cells. In this case it is not possible to continue */
3503 /* tracing the convex contour within the cell. This situation appears also in transition cells of multi-resolution grids where */
3504 /* the isopatch degenerates into a plane. */
3505 /***********************************************************************************************************************/
3506
3507 else if (((data_vertex1 < isovalue) && (isovalue > data_vertex2)) || ((data_vertex1 > isovalue) && (isovalue < data_vertex2)))
3508 {
3509 abort_tracing_isocontour = true;
3510 }
3511 }
3512}
3513
3514void generate_tesselation(TESSELATION &triangulation, CONTOUR contour, ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector)
3515{
3516 bool end_of_triangulation;
3517
3518 int i;
3519 int j;
3520 int polygon_index;
3521 int temp_polygon_index;
3522 int predecessor;
3523 int successor;
3524 int triangulation_counter;
3525 int normal_vector_counter;
3526
3527 TRIANGLE tesselation_triangle;
3528 TRIANGLE temp_triangle;
3529
3530 EDGE_VECTOR vector1;
3531 EDGE_VECTOR vector2;
3532 EDGE_VECTOR normal;
3533
3534 POLYGON polygon;
3535
3536 POLYGON_ITERATOR start;
3537 POLYGON_ITERATOR end;
3538 POLYGON_ITERATOR new_end;
3539
3540 /* Avoid Unnecessary Reallocations */
3541 polygon.reserve(15);
3542
3543 /*******************************************/
3544 /* Triangulate the convex contour(s) found */
3545 /*******************************************/
3546
3547 for (i = 0; i < ssize_t(contour.ring_index.size()); i++)
3548 {
3549 /* Make sure polygon and vertex containers are empty */
3550 polygon.clear();
3551
3552 /* Analyzed contour is not the last one */
3553 if (i < ssize_t(contour.ring_index.size()) - 1)
3554 {
3555 /**************************************************/
3556 /* Case A - --> The analyzed contour is a triangle */
3557 /**************************************************/
3558
3559 if ((contour.ring_index[i + 1] - contour.ring_index[i]) == 3)
3560 {
3561 /* Generate the triangle */
3562 tesselation_triangle.vertex1 = contour.ring[contour.ring_index[i]];
3563 tesselation_triangle.vertex2 = contour.ring[contour.ring_index[i] + 1];
3564 tesselation_triangle.vertex3 = contour.ring[contour.ring_index[i] + 2];
3565
3566 /* Define two edge vectors of the triangle */
3567 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3568 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3569 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3570
3571 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3572 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3573 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3574
3575 /* Calculate normal to the triangle */
3576 normal = cross_product(vector1, vector2);
3577 tesselation_triangle.normal.x = normal.x;
3578 tesselation_triangle.normal.y = normal.y;
3579 tesselation_triangle.normal.z = normal.z;
3580
3581 /* Store the triangle in the tesselation vector */
3582 triangulation.push_back(tesselation_triangle);
3583 }
3584
3585 /***********************************************************/
3586 /* Case B ---> The analyzed contour is an arbitrary polygon */
3587 /***********************************************************/
3588
3589 else if ((contour.ring_index[i + 1] - contour.ring_index[i]) > 3)
3590 {
3591 for (j = contour.ring_index[i]; j < contour.ring_index[i + 1]; j++)
3592 {
3593 /* Define a polygon data structure with the current ring */
3594 polygon.push_back(contour.ring[j]);
3595 }
3596
3597 /************************/
3598 /* Initiate Graham Scan */
3599 /************************/
3600
3601 polygon_index = 1;
3602 end_of_triangulation = false;
3603 triangulation_counter = 0;
3604
3605 do
3606 {
3607 if (polygon.size() > 3)
3608 {
3609 /***********************************************/
3610 /* Update node predecessors and successors */
3611 /***********************************************/
3612
3613 if (polygon_index == 1)
3614 {
3615 predecessor = polygon_index - 1;
3616 successor = polygon_index + 1;
3617 }
3618
3619 else if (polygon_index == 0)
3620 {
3621 predecessor = (int)polygon.size() - 1;
3622 successor = polygon_index + 1;
3623 }
3624
3625 else if (polygon_index == ssize_t(polygon.size() - 1))
3626 {
3627 predecessor = polygon_index - 1;
3628 successor = 0;
3629 }
3630
3631 else if (polygon_index == ssize_t(polygon.size()))
3632 {
3633 polygon_index = 0;
3634 predecessor = (int)polygon.size() - 1;
3635 successor = polygon_index + 1;
3636 }
3637
3638 else
3639 {
3640 predecessor = polygon_index - 1;
3641 successor = polygon_index + 1;
3642 }
3643
3644 /* Define a triangle */
3645 tesselation_triangle.vertex1 = polygon[predecessor];
3646 tesselation_triangle.vertex2 = polygon[polygon_index];
3647 tesselation_triangle.vertex3 = polygon[successor];
3648
3649 /* Define two edge vectors of the triangle */
3650 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3651 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3652 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3653
3654 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3655 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3656 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3657
3658 /* Store initial triangle */
3659 temp_polygon_index = polygon_index;
3660 temp_triangle.vertex1 = tesselation_triangle.vertex1;
3661 temp_triangle.vertex2 = tesselation_triangle.vertex2;
3662 temp_triangle.vertex3 = tesselation_triangle.vertex3;
3663
3664 /* Calculate normal to the triangle */
3665 normal = cross_product(vector1, vector2);
3666 tesselation_triangle.normal.x = normal.x;
3667 tesselation_triangle.normal.y = normal.y;
3668 tesselation_triangle.normal.z = normal.z;
3669
3670 temp_triangle.normal.x = normal.x;
3671 temp_triangle.normal.y = normal.y;
3672 temp_triangle.normal.z = normal.z;
3673
3674 normal_vector_counter = 0;
3675
3676 /*******************************************************************************************************/
3677 /* Special Case ---> Normal is equal to the null vector */
3678 /* */
3679 /* This result can only be obtained when two edges of the triangle are collinear (degenerate triangle) . */
3680 /* This implies that the cell is too small or that at least two adjacent intersections are almost coincident or */
3681 /* actually overlap. */
3682 /*******************************************************************************************************/
3683
3684 if (normal.x == 0 && normal.y == 0 && normal.z == 0)
3685 {
3686 do
3687 {
3688 normal_vector_counter++;
3689
3690 /* Update Scan */
3691 if (polygon_index + 1 <= ssize_t(polygon.size()))
3692 {
3693 polygon_index++;
3694 }
3695
3696 /***********************************************/
3697 /* Update node predecessors and successors */
3698 /***********************************************/
3699
3700 if (polygon_index == 0)
3701 {
3702 predecessor = (int)polygon.size() - 1;
3703 successor = polygon_index + 1;
3704 }
3705
3706 else if (polygon_index == ssize_t(polygon.size() - 1))
3707 {
3708 predecessor = polygon_index - 1;
3709 successor = 0;
3710 }
3711
3712 else if (polygon_index == ssize_t(polygon.size()))
3713 {
3714 polygon_index = 0;
3715 predecessor = (int)polygon.size() - 1;
3716 successor = polygon_index + 1;
3717 }
3718
3719 else if (polygon_index != 0 && polygon_index != ssize_t(polygon.size()) - 1 && polygon_index != ssize_t(polygon.size()))
3720 {
3721 predecessor = polygon_index - 1;
3722 successor = polygon_index + 1;
3723 }
3724
3725 /* Define a new triangle */
3726 tesselation_triangle.vertex1 = polygon[predecessor];
3727 tesselation_triangle.vertex2 = polygon[polygon_index];
3728 tesselation_triangle.vertex3 = polygon[successor];
3729
3730 /* Define two edge vectors of the triangle */
3731 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3732 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3733 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3734
3735 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3736 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3737 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3738
3739 /* Recalculate normal to the triangle */
3740 normal = cross_product(vector1, vector2);
3741 tesselation_triangle.normal.x = normal.x;
3742 tesselation_triangle.normal.y = normal.y;
3743 tesselation_triangle.normal.z = normal.z;
3744 } while ((normal.x == 0 && normal.y == 0 && normal.z == 0) && normal_vector_counter < ssize_t(polygon.size()));
3745 }
3746
3747 /* If normal still equals the null vector implies that the cell is too small */
3748 /* Proceed to store the last triangle in the tesselation vector */
3749 triangulation_counter = (int)polygon.size() + 1;
3750
3751 /* Avoid infinite loops with the triangulation counter */
3752 if (triangulation_counter > (int)polygon.size())
3753 {
3754 /* Store the triangle in the triangulation vector */
3755 triangulation.push_back(temp_triangle);
3756
3757 /* Cut ear from the polygon vector */
3758 start = polygon.begin();
3759 end = polygon.end();
3760 new_end = remove(start, end, polygon[temp_polygon_index]);
3761 polygon.erase(new_end, end);
3762
3763 /* Update scan */
3764 polygon_index = 1;
3765 }
3766 }
3767
3768 /* All ears have been cut; only one triangle is left */
3769 if (polygon.size() == 3)
3770 {
3771 /* Update node predecessors and successors */
3772 if (polygon_index == 1)
3773 {
3774 predecessor = polygon_index - 1;
3775 successor = polygon_index + 1;
3776 }
3777
3778 else if (polygon_index == 0)
3779 {
3780 predecessor = (int)polygon.size() - 1;
3781 successor = polygon_index + 1;
3782 }
3783
3784 else if (polygon_index == (int)polygon.size() - 1)
3785 {
3786 predecessor = polygon_index - 1;
3787 successor = 0;
3788 }
3789
3790 else if (polygon_index == (int)polygon.size())
3791 {
3792 polygon_index = 0;
3793 predecessor = (int)polygon.size() - 1;
3794 successor = polygon_index + 1;
3795 }
3796
3797 else
3798 {
3799 predecessor = polygon_index - 1;
3800 successor = polygon_index + 1;
3801 }
3802
3803 /* Store the last triangle in the tesselation vector */
3804 tesselation_triangle.vertex1 = polygon[predecessor];
3805 tesselation_triangle.vertex2 = polygon[polygon_index];
3806 tesselation_triangle.vertex3 = polygon[successor];
3807
3808 /* Define two edge vectors of the triangle */
3809 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3810 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3811 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3812
3813 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3814 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3815 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3816
3817 /* Calculate normal to the triangle */
3818 normal = cross_product(vector1, vector2);
3819 tesselation_triangle.normal.x = normal.x;
3820 tesselation_triangle.normal.y = normal.y;
3821 tesselation_triangle.normal.z = normal.z;
3822
3823 /* Store the triangle in the tesselation vector */
3824 triangulation.push_back(tesselation_triangle);
3825 end_of_triangulation = true;
3826 }
3827 } while (end_of_triangulation != true);
3828 }
3829 }
3830
3831 /* Analyzed contour is the last one */
3832 else
3833 {
3834 /**************************************************/
3835 /* Case A - --> The analyzed contour is a triangle */
3836 /**************************************************/
3837
3838 if ((contour.ring.size() - contour.ring_index[i]) == 3)
3839 {
3840 /* Generate the triangle */
3841 tesselation_triangle.vertex1 = contour.ring[contour.ring_index[i]];
3842 tesselation_triangle.vertex2 = contour.ring[contour.ring_index[i] + 1];
3843 tesselation_triangle.vertex3 = contour.ring[contour.ring_index[i] + 2];
3844
3845 /* Define two edge vectors of the triangle */
3846 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3847 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3848 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3849
3850 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3851 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3852 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3853
3854 /* Calculate normal to the triangle */
3855 normal = cross_product(vector1, vector2);
3856 tesselation_triangle.normal.x = normal.x;
3857 tesselation_triangle.normal.y = normal.y;
3858 tesselation_triangle.normal.z = normal.z;
3859
3860 /* Store the triangle in the tesselation vector */
3861 triangulation.push_back(tesselation_triangle);
3862 }
3863
3864 /***********************************************************/
3865 /* Case B ---> The analyzed contour is an arbitrary polygon */
3866 /***********************************************************/
3867
3868 else if ((contour.ring.size() - contour.ring_index[i]) > 3)
3869 {
3870 for (j = contour.ring_index[i]; j < ssize_t(contour.ring.size()); j++)
3871 {
3872 /* Define a polygon data structure with the current convex contour */
3873 polygon.push_back(contour.ring[j]);
3874 }
3875
3876 /************************/
3877 /* Initiate Graham Scan */
3878 /************************/
3879
3880 polygon_index = 1;
3881 end_of_triangulation = false;
3882 triangulation_counter = 0;
3883
3884 do
3885 {
3886 if (polygon.size() > 3)
3887 {
3888 /***********************************************/
3889 /* Update node predecessors and successors */
3890 /***********************************************/
3891
3892 if (polygon_index == 1)
3893 {
3894 predecessor = polygon_index - 1;
3895 successor = polygon_index + 1;
3896 }
3897
3898 else if (polygon_index == 0)
3899 {
3900 predecessor = (int)polygon.size() - 1;
3901 successor = polygon_index + 1;
3902 }
3903
3904 else if (polygon_index == (int)polygon.size() - 1)
3905 {
3906 predecessor = polygon_index - 1;
3907 successor = 0;
3908 }
3909
3910 else if (polygon_index == (int)polygon.size())
3911 {
3912 polygon_index = 0;
3913 predecessor = (int)polygon.size() - 1;
3914 successor = polygon_index + 1;
3915 }
3916
3917 else
3918 {
3919 predecessor = polygon_index - 1;
3920 successor = polygon_index + 1;
3921 }
3922
3923 /* Define a triangle */
3924 tesselation_triangle.vertex1 = polygon[predecessor];
3925 tesselation_triangle.vertex2 = polygon[polygon_index];
3926 tesselation_triangle.vertex3 = polygon[successor];
3927
3928 /* Define two edge vectors of the triangle */
3929 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3930 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3931 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3932
3933 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
3934 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
3935 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
3936
3937 /* Store initial triangle */
3938 temp_polygon_index = polygon_index;
3939 temp_triangle.vertex1 = tesselation_triangle.vertex1;
3940 temp_triangle.vertex2 = tesselation_triangle.vertex2;
3941 temp_triangle.vertex3 = tesselation_triangle.vertex3;
3942
3943 /* Calculate normal to the triangle */
3944 normal = cross_product(vector1, vector2);
3945 tesselation_triangle.normal.x = normal.x;
3946 tesselation_triangle.normal.y = normal.y;
3947 tesselation_triangle.normal.z = normal.z;
3948
3949 temp_triangle.normal.x = normal.x;
3950 temp_triangle.normal.y = normal.y;
3951 temp_triangle.normal.z = normal.z;
3952
3953 normal_vector_counter = 0;
3954
3955 /*******************************************************************************************************/
3956 /* Special Case ---> Normal is equal to the null vector */
3957 /* */
3958 /* This result can only be obtained when two edges of the triangle are collinear (degenerate triangle). */
3959 /* This implies that the cell is too small or that at least two adjacent intersections are almost coincident or */
3960 /* actually overlap. */
3961 /*******************************************************************************************************/
3962
3963 if (normal.x == 0 && normal.y == 0 && normal.z == 0)
3964 {
3965 do
3966 {
3967 normal_vector_counter++;
3968
3969 /* Update Scan */
3970 if (polygon_index + 1 <= ssize_t(polygon.size()))
3971 {
3972 polygon_index++;
3973 }
3974
3975 /***********************************************/
3976 /* Update node predecessors and successors */
3977 /***********************************************/
3978
3979 if (polygon_index == 1)
3980 {
3981 predecessor = polygon_index - 1;
3982 successor = polygon_index + 1;
3983 }
3984
3985 else if (polygon_index == 0)
3986 {
3987 predecessor = (int)polygon.size() - 1;
3988 successor = polygon_index + 1;
3989 }
3990
3991 else if (polygon_index == (int)polygon.size() - 1)
3992 {
3993 predecessor = polygon_index - 1;
3994 successor = 0;
3995 }
3996
3997 else if (polygon_index == (int)polygon.size())
3998 {
3999 polygon_index = 0;
4000 predecessor = (int)polygon.size() - 1;
4001 successor = polygon_index + 1;
4002 }
4003
4004 else
4005 {
4006 predecessor = polygon_index - 1;
4007 successor = polygon_index + 1;
4008 }
4009
4010 /* Define a new triangle */
4011 tesselation_triangle.vertex1 = polygon[predecessor];
4012 tesselation_triangle.vertex2 = polygon[polygon_index];
4013 tesselation_triangle.vertex3 = polygon[successor];
4014
4015 /* Define two edge vectors of the triangle */
4016 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
4017 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
4018 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
4019
4020 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
4021 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
4022 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
4023
4024 /* Recalculate normal to the triangle */
4025 normal = cross_product(vector1, vector2);
4026
4027 tesselation_triangle.normal.x = normal.x;
4028 tesselation_triangle.normal.y = normal.y;
4029 tesselation_triangle.normal.z = normal.z;
4030 } while ((normal.x == 0 && normal.y == 0 && normal.z == 0) && normal_vector_counter < ssize_t(polygon.size()));
4031 }
4032
4033 /* If normal still equals the null vector implies that the cell is too small */
4034 /* Proceed to store the last triangle in the tesselation vector */
4035 triangulation_counter = (int)polygon.size() + 1;
4036
4037 /* Avoid infinite loops with the triangulation counter */
4038 if (triangulation_counter > (int)polygon.size())
4039 {
4040 /* Store the triangle in the triangulation vector */
4041 triangulation.push_back(temp_triangle);
4042
4043 /* Cut ear from the polygon vector */
4044 start = polygon.begin();
4045 end = polygon.end();
4046 new_end = remove(start, end, polygon[temp_polygon_index]);
4047 polygon.erase(new_end, end);
4048
4049 /* Update scan */
4050 polygon_index = 1;
4051 }
4052 }
4053
4054 /* All ears have been cut; only one triangle is left */
4055 if (polygon.size() == 3)
4056 {
4057 /***********************************************/
4058 /* Update node predecessors and successors */
4059 /***********************************************/
4060
4061 if (polygon_index == 1)
4062 {
4063 predecessor = polygon_index - 1;
4064 successor = polygon_index + 1;
4065 }
4066
4067 if (polygon_index == 0)
4068 {
4069 predecessor = (int)polygon.size() - 1;
4070 successor = polygon_index + 1;
4071 }
4072
4073 else if (polygon_index == ssize_t(polygon.size() - 1))
4074 {
4075 predecessor = polygon_index - 1;
4076 successor = 0;
4077 }
4078
4079 else if (polygon_index == ssize_t(polygon.size()))
4080 {
4081 polygon_index = 0;
4082 predecessor = (int)polygon.size() - 1;
4083 successor = polygon_index + 1;
4084 }
4085
4086 else
4087 {
4088 predecessor = polygon_index - 1;
4089 successor = polygon_index + 1;
4090 }
4091
4092 /* Store the last triangle in the tesselation vector */
4093 tesselation_triangle.vertex1 = polygon[predecessor];
4094 tesselation_triangle.vertex2 = polygon[polygon_index];
4095 tesselation_triangle.vertex3 = polygon[successor];
4096
4097 /* Define two edge vectors of the triangle */
4098 vector1.x = intsec_vector[tesselation_triangle.vertex2].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
4099 vector1.y = intsec_vector[tesselation_triangle.vertex2].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
4100 vector1.z = intsec_vector[tesselation_triangle.vertex2].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
4101
4102 vector2.x = intsec_vector[tesselation_triangle.vertex3].intersection.x - intsec_vector[tesselation_triangle.vertex1].intersection.x;
4103 vector2.y = intsec_vector[tesselation_triangle.vertex3].intersection.y - intsec_vector[tesselation_triangle.vertex1].intersection.y;
4104 vector2.z = intsec_vector[tesselation_triangle.vertex3].intersection.z - intsec_vector[tesselation_triangle.vertex1].intersection.z;
4105
4106 /* Calculate normal to the triangle */
4107 normal = cross_product(vector1, vector2);
4108 tesselation_triangle.normal.x = normal.x;
4109 tesselation_triangle.normal.y = normal.y;
4110 tesselation_triangle.normal.z = normal.z;
4111
4112 /* Store the triangle in the tesselation vector */
4113 triangulation.push_back(tesselation_triangle);
4114 end_of_triangulation = true;
4115 }
4116 } while (end_of_triangulation != true);
4117 }
4118 }
4119 }
4120}
4121}
4122#endif
NULL
#define NULL
Definition: covise_list.h:22
start
GLuint start
Definition: khronos-glext.h:6343
y2
GLfixed GLfixed GLfixed y2
Definition: khronos-glext.h:11325
y1
GLfixed y1
Definition: khronos-glext.h:11325
end
GLuint GLuint end
Definition: khronos-glext.h:6343
v
const GLdouble * v
Definition: khronos-glext.h:6442
x1
GLuint GLfloat GLfloat GLfloat x1
Definition: khronos-glext.h:13144
index
GLuint index
Definition: khronos-glext.h:6722
x2
GLfixed GLfixed x2
Definition: khronos-glext.h:11325
v1
GLfloat GLfloat v1
Definition: khronos-glext.h:6753
v2
GLfloat GLfloat GLfloat v2
Definition: khronos-glext.h:6754
fabs
__host__ __device__ float2 fabs(float2 v)
Definition: cutil_math.h:1380
covise
list of all chemical elements
Definition: coConfig.h:27
covise::generate_tesselation
void generate_tesselation(TESSELATION &triangulation, CONTOUR contour, ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector)
Definition: IsoSurfaceGPMUtil.h:3514
covise::CONNECTIVITY_VECTOR
std::vector< int > CONNECTIVITY_VECTOR
Definition: IsoSurfaceGPMUtil.h:77
covise::cross_product
EDGE_VECTOR cross_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:59
covise::VertexInterpolate
ISOSURFACE_EDGE_INTERSECTION VertexInterpolate(float x1, float y1, float z1, float x2, float y2, float z2, float isovalue, float data1, float data2, int v1, int v2)
Definition: IsoSurfaceGPMUtil.h:108
covise::vec_length
double vec_length(EDGE_VECTOR &vector)
Definition: IsoSurfaceGPMUtil.h:103
covise::ITERATOR
std::vector< int >::iterator ITERATOR
Definition: IsoSurfaceGPMUtil.h:82
covise::find_current_face
void find_current_face(CONTOUR &contour, ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, int &edge_vertex1, int &edge_vertex2, float &data_vertex1, float &data_vertex2, float *isodata_in, int *elem_in, int *conn_in, int *index_list, int *polygon_list, int num_coord_in, int num_conn_in, int num_elem_in, int &ring_counter, int &current_face, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool improper_topology, bool &abort_tracing_isocontour)
Definition: IsoSurfaceGPMUtil.h:649
covise::map_to_isosurface
float map_to_isosurface(float coord_x1, float coord_x2, float coord_y1, float coord_y2, float coord_z1, float coord_z2, float coord_isox, float coord_isoy, float coord_isoz, float data_1, float data_2, bool int_vertex1, bool int_vertex2)
Definition: IsoSurfaceGPMUtil.h:304
covise::assign_int_index
int assign_int_index(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
Definition: CuttingSurfaceGPMUtil.h:469
covise::ISOSURFACE_EDGE_INTERSECTION_VECTOR
std::vector< ISOSURFACE_EDGE_INTERSECTION > ISOSURFACE_EDGE_INTERSECTION_VECTOR
Definition: IsoSurfaceGPMUtil.h:73
covise::test_intersection
bool test_intersection(ISOSURFACE_EDGE_INTERSECTION_VECTOR &intsec_vector, ISOSURFACE_EDGE_INTERSECTION &intsec, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool &improper_topology)
Definition: IsoSurfaceGPMUtil.h:209
covise::TESSELATION
std::vector< TRIANGLE > TESSELATION
Definition: IsoSurfaceGPMUtil.h:75
covise::generate_isocontour
void generate_isocontour(ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, float data_vertex1, float data_vertex2, int edge_vertex1, int edge_vertex2, int &new_edge_vertex1, int &new_edge_vertex2, int *elem_in, int *conn_in, float isovalue, int num_elem_in, int num_conn_in, int current_face, float *x_coord_in, float *y_coord_in, float *z_coord_in, bool improper_topology, bool &abort_tracing_isocontour, CONTOUR contour, int num_of_rings, int ring_end)
Definition: IsoSurfaceGPMUtil.h:2572
covise::dot_product
double dot_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:54
covise::POLYGON_ITERATOR
std::vector< int >::iterator POLYGON_ITERATOR
Definition: IsoSurfaceGPMUtil.h:81
covise::POLYGON
std::vector< int > POLYGON
Definition: IsoSurfaceGPMUtil.h:78
covise::find_intersection
bool find_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
Definition: CuttingSurfaceGPMUtil.h:512
covise::calculate_intersections
PLANE_EDGE_INTERSECTION_VECTOR calculate_intersections(float *udata_in, float *vdata_in, float *wdata_in, int num_elem_in, int *elem_in, int num_conn_in, int *conn_in, float *x_coord_in, float *y_coord_in, float *z_coord_in, float p, EDGE_VECTOR unit_normal_vector)
Definition: CuttingSurfaceGPMUtil.h:300
covise::PROCESSED_ELEMENTS
std::vector< int > PROCESSED_ELEMENTS
Definition: IsoSurfaceGPMUtil.h:79
covise::S_V_DATA
Definition: CuttingSurfaceGPMUtil.h:14
covise::S_V_DATA::dimension
int dimension
Definition: CuttingSurfaceGPMUtil.h:16
covise::S_V_DATA::v
float v[3]
Definition: CuttingSurfaceGPMUtil.h:15
covise::EDGE_VECTOR
Definition: CuttingSurfaceGPMUtil.h:20
covise::EDGE_VECTOR::x
float x
Definition: IsoSurfaceGPMUtil.h:36
covise::EDGE_VECTOR::x
double x
Definition: CuttingSurfaceGPMUtil.h:21
covise::EDGE_VECTOR::z
double z
Definition: CuttingSurfaceGPMUtil.h:23
covise::EDGE_VECTOR::y
float y
Definition: IsoSurfaceGPMUtil.h:37
covise::EDGE_VECTOR::y
double y
Definition: CuttingSurfaceGPMUtil.h:22
covise::EDGE_VECTOR::z
float z
Definition: IsoSurfaceGPMUtil.h:38
covise::CONTOUR
Definition: CuttingSurfaceGPMUtil.h:38
covise::CONTOUR::ring_index
std::vector< int > ring_index
Definition: IsoSurfaceGPMUtil.h:69
covise::CONTOUR::polyhedron_faces
std::vector< int > polyhedron_faces
Definition: IsoSurfaceGPMUtil.h:70
covise::CONTOUR::ring_index
vector< int > ring_index
Definition: CuttingSurfaceGPMUtil.h:40
covise::CONTOUR::ring
std::vector< int > ring
Definition: IsoSurfaceGPMUtil.h:68
covise::CONTOUR::polyhedron_faces
vector< int > polyhedron_faces
Definition: CuttingSurfaceGPMUtil.h:41
covise::CONTOUR::ring
vector< int > ring
Definition: CuttingSurfaceGPMUtil.h:39
covise::TRIANGLE
Definition: IsoSurfaceGPMUtil.h:42
covise::TRIANGLE::vertex3
int vertex3
Definition: IsoSurfaceGPMUtil.h:45
covise::TRIANGLE::normal
EDGE_VECTOR normal
Definition: IsoSurfaceGPMUtil.h:47
covise::TRIANGLE::vertex2
int vertex2
Definition: IsoSurfaceGPMUtil.h:44
covise::TRIANGLE::vertex1
int vertex1
Definition: IsoSurfaceGPMUtil.h:43
covise::ISOSURFACE_EDGE_INTERSECTION
Definition: IsoSurfaceGPMUtil.h:51
covise::ISOSURFACE_EDGE_INTERSECTION::intersection_at_vertex2
bool intersection_at_vertex2
Definition: IsoSurfaceGPMUtil.h:53
covise::ISOSURFACE_EDGE_INTERSECTION::vertex1
int vertex1
Definition: IsoSurfaceGPMUtil.h:56
covise::ISOSURFACE_EDGE_INTERSECTION::data_vertex_int
S_V_DATA data_vertex_int
Definition: IsoSurfaceGPMUtil.h:62
covise::ISOSURFACE_EDGE_INTERSECTION::intersection
EDGE_VECTOR intersection
Definition: IsoSurfaceGPMUtil.h:63
covise::ISOSURFACE_EDGE_INTERSECTION::data_vertex1
float data_vertex1
Definition: IsoSurfaceGPMUtil.h:59
covise::ISOSURFACE_EDGE_INTERSECTION::vertex2
int vertex2
Definition: IsoSurfaceGPMUtil.h:57
covise::ISOSURFACE_EDGE_INTERSECTION::data_vertex2
float data_vertex2
Definition: IsoSurfaceGPMUtil.h:60
covise::ISOSURFACE_EDGE_INTERSECTION::intersection_at_vertex1
bool intersection_at_vertex1
Definition: IsoSurfaceGPMUtil.h:52
covise::ISOSURFACE_EDGE_INTERSECTION::int_flag
int int_flag
Definition: IsoSurfaceGPMUtil.h:55
Generated by doxygen 1.9.4