COVISE Core
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros
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 /************************************************************/
25 namespace covise
26 {
27 
28 typedef struct
29 {
30  float v[3];
31  int dimension;
32 } S_V_DATA;
33 
34 typedef struct
35 {
36  float x;
37  float y;
38  float z;
39 } EDGE_VECTOR;
40 
41 typedef struct
42 {
43  int vertex1;
44  int vertex2;
45  int vertex3;
46 
47  EDGE_VECTOR normal;
48 } TRIANGLE;
49 
50 typedef 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 
66 typedef struct
67 {
68  std::vector<int> ring;
69  std::vector<int> ring_index;
70  std::vector<int> polyhedron_faces;
71 } CONTOUR;
72 
73 typedef std::vector<ISOSURFACE_EDGE_INTERSECTION> ISOSURFACE_EDGE_INTERSECTION_VECTOR;
74 
75 typedef std::vector<TRIANGLE> TESSELATION;
76 
77 typedef std::vector<int> CONNECTIVITY_VECTOR;
78 typedef std::vector<int> POLYGON;
79 typedef std::vector<int> PROCESSED_ELEMENTS;
80 
81 typedef std::vector<int>::iterator POLYGON_ITERATOR;
82 typedef std::vector<int>::iterator ITERATOR;
83 
84 float 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 
92 EDGE_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 
103 double vec_length(EDGE_VECTOR &vector)
104 {
105  return sqrt(vector.x * vector.x + vector.y * vector.y + vector.z * vector.z);
106 }
107 
108 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)
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 
209 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)
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 
304 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)
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) + pow(coord_y1 - coord_y2, 2) + pow(coord_z1 - coord_z2, 2));
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) + pow(coord_y1 - coord_isoy, 2) + pow(coord_z1 - coord_isoz, 2));
333  mapped_value = data_1 + ((data_2 - data_1) / dist_x1x2) * dist_x1xiso;
334  }
335  }
336 
337  return mapped_value;
338 }
339 
340 ISOSURFACE_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 
528 int assign_int_index(ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector, int edge_vertex1, int edge_vertex2)
529 {
530  int 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 = i;
540  break;
541  }
542 
543  if ((edge_vertex2 == intsec_vector[i].vertex1) && (edge_vertex1 == intsec_vector[i].vertex2))
544  {
545  index = i;
546  break;
547  }
548  }
549 
550  return index;
551 }
552 
553 bool 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 < 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 
649 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)
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 
2572 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)
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 
3514 void 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 < 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 < 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 == polygon.size() - 1)
3626  {
3627  predecessor = polygon_index - 1;
3628  successor = 0;
3629  }
3630 
3631  else if (polygon_index == 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 <= 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 == polygon.size() - 1)
3707  {
3708  predecessor = polygon_index - 1;
3709  successor = 0;
3710  }
3711 
3712  else if (polygon_index == 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 != polygon.size() - 1 && polygon_index != 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 < 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 < 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 <= 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 < 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 == polygon.size() - 1)
4074  {
4075  predecessor = polygon_index - 1;
4076  successor = 0;
4077  }
4078 
4079  else if (polygon_index == 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
covise::ITERATOR
std::vector< int >::iterator ITERATOR
Definition: IsoSurfaceGPMUtil.h:82
covise::ISOSURFACE_EDGE_INTERSECTION::intersection_at_vertex2
bool intersection_at_vertex2
Definition: IsoSurfaceGPMUtil.h:53
covise::POLYGON
std::vector< int > POLYGON
Definition: IsoSurfaceGPMUtil.h:78
covise::ISOSURFACE_EDGE_INTERSECTION::data_vertex1
float data_vertex1
Definition: IsoSurfaceGPMUtil.h:59
covise::CONTOUR::ring
std::vector< int > ring
Definition: IsoSurfaceGPMUtil.h:68
covise::CONTOUR::ring_index
vector< int > ring_index
Definition: CuttingSurfaceGPMUtil.h:40
covise::ISOSURFACE_EDGE_INTERSECTION::vertex1
int vertex1
Definition: IsoSurfaceGPMUtil.h:56
covise::EDGE_VECTOR::x
double x
Definition: CuttingSurfaceGPMUtil.h:21
v1
GLfloat GLfloat v1
Definition: khronos-glext.h:6753
covise::TRIANGLE::vertex3
int vertex3
Definition: IsoSurfaceGPMUtil.h:45
covise::CONTOUR
Definition: CuttingSurfaceGPMUtil.h:37
covise::CONTOUR::polyhedron_faces
vector< int > polyhedron_faces
Definition: CuttingSurfaceGPMUtil.h:41
covise::vec_length
double vec_length(EDGE_VECTOR &vector)
Definition: IsoSurfaceGPMUtil.h:103
end
GLuint GLuint end
Definition: khronos-glext.h:6343
covise::ISOSURFACE_EDGE_INTERSECTION::int_flag
int int_flag
Definition: IsoSurfaceGPMUtil.h:55
fabs
__host__ __device__ float2 fabs(float2 v)
Definition: cutil_math.h:1380
covise::EDGE_VECTOR::y
double y
Definition: CuttingSurfaceGPMUtil.h:22
covise::ISOSURFACE_EDGE_INTERSECTION::data_vertex2
float data_vertex2
Definition: IsoSurfaceGPMUtil.h:60
covise::ISOSURFACE_EDGE_INTERSECTION_VECTOR
std::vector< ISOSURFACE_EDGE_INTERSECTION > ISOSURFACE_EDGE_INTERSECTION_VECTOR
Definition: IsoSurfaceGPMUtil.h:73
covise::TRIANGLE
Definition: IsoSurfaceGPMUtil.h:41
covise::cross_product
EDGE_VECTOR cross_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:59
covise::CONTOUR::polyhedron_faces
std::vector< int > polyhedron_faces
Definition: IsoSurfaceGPMUtil.h:70
covise::CONNECTIVITY_VECTOR
std::vector< int > CONNECTIVITY_VECTOR
Definition: IsoSurfaceGPMUtil.h:77
NULL
#define NULL
Definition: covise_list.h:22
covise::EDGE_VECTOR::z
double z
Definition: CuttingSurfaceGPMUtil.h:23
v
const GLdouble * v
Definition: khronos-glext.h:6442
covise::POLYGON_ITERATOR
std::vector< int >::iterator POLYGON_ITERATOR
Definition: IsoSurfaceGPMUtil.h:81
covise::TRIANGLE::vertex1
int vertex1
Definition: IsoSurfaceGPMUtil.h:43
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::S_V_DATA::v
float v[3]
Definition: CuttingSurfaceGPMUtil.h:15
covise::find_intersection
bool find_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
Definition: CuttingSurfaceGPMUtil.h:512
covise::EDGE_VECTOR::z
float z
Definition: IsoSurfaceGPMUtil.h:38
covise::ISOSURFACE_EDGE_INTERSECTION::data_vertex_int
S_V_DATA data_vertex_int
Definition: IsoSurfaceGPMUtil.h:62
start
GLuint start
Definition: khronos-glext.h:6343
covise::generate_tesselation
void generate_tesselation(TESSELATION &triangulation, CONTOUR contour, ISOSURFACE_EDGE_INTERSECTION_VECTOR intsec_vector)
Definition: IsoSurfaceGPMUtil.h:3514
y2
GLfixed GLfixed GLfixed y2
Definition: khronos-glext.h:11325
covise::CONTOUR::ring_index
std::vector< int > ring_index
Definition: IsoSurfaceGPMUtil.h:69
covise::TRIANGLE::vertex2
int vertex2
Definition: IsoSurfaceGPMUtil.h:44
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
index
GLuint index
Definition: khronos-glext.h:6722
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
x2
GLfixed GLfixed x2
Definition: khronos-glext.h:11325
covise::TESSELATION
std::vector< TRIANGLE > TESSELATION
Definition: IsoSurfaceGPMUtil.h:75
covise::TRIANGLE::normal
EDGE_VECTOR normal
Definition: IsoSurfaceGPMUtil.h:47
x1
GLuint GLfloat GLfloat GLfloat x1
Definition: khronos-glext.h:13144
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::ISOSURFACE_EDGE_INTERSECTION::intersection
EDGE_VECTOR intersection
Definition: IsoSurfaceGPMUtil.h:63
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::ISOSURFACE_EDGE_INTERSECTION::intersection_at_vertex1
bool intersection_at_vertex1
Definition: IsoSurfaceGPMUtil.h:52
covise::PROCESSED_ELEMENTS
std::vector< int > PROCESSED_ELEMENTS
Definition: IsoSurfaceGPMUtil.h:79
covise::ISOSURFACE_EDGE_INTERSECTION::vertex2
int vertex2
Definition: IsoSurfaceGPMUtil.h:57
covise::S_V_DATA
Definition: CuttingSurfaceGPMUtil.h:13
y1
GLfixed y1
Definition: khronos-glext.h:11325
covise::dot_product
double dot_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:54
covise::S_V_DATA::dimension
int dimension
Definition: CuttingSurfaceGPMUtil.h:16
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
Definition: IsoSurfaceGPMUtil.h:50
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::EDGE_VECTOR
Definition: CuttingSurfaceGPMUtil.h:19
covise::EDGE_VECTOR::y
float y
Definition: IsoSurfaceGPMUtil.h:37
covise::EDGE_VECTOR::x
float x
Definition: IsoSurfaceGPMUtil.h:36
v2
GLfloat GLfloat GLfloat v2
Definition: khronos-glext.h:6754
covise::CONTOUR::ring
vector< int > ring
Definition: CuttingSurfaceGPMUtil.h:39
Generated by   doxygen 1.8.5