COVISE Core
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros
CuttingSurfaceGPMUtil.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 #include <do/Triangulate.h>
9 
10 namespace covise
11 {
12 
13 typedef struct
14 {
15  float v[3];
16  int dimension;
17 } S_V_DATA;
18 
19 typedef struct
20 {
21  double x;
22  double y;
23  double z;
24 } EDGE_VECTOR;
25 typedef struct
26 {
29 
30  int int_flag;
31  int vertex1;
32  int vertex2;
33 
37 typedef struct
38 {
39  vector<int> ring;
40  vector<int> ring_index;
41  vector<int> polyhedron_faces;
42 } CONTOUR;
43 
44 // typedef struct
45 // {
46 // vector<int> element_list;
47 // vector<int> connectivity_list;
48 // vector<float> plane_data;
49 // }PLANE_CELL_INTERSECTION ;
50 
51 typedef std::vector<PLANE_EDGE_INTERSECTION> PLANE_EDGE_INTERSECTION_VECTOR;
52 //typedef std::vector<PLANE_CELL_INTERSECTION> PLANE_CELL_INTERSECTION_VECTOR;
53 
54 double dot_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
55 {
56  return vector1.x * vector2.x + vector1.y * vector2.y + vector1.z * vector2.z;
57 }
58 
60 {
61  EDGE_VECTOR normal;
62 
63  normal.x = (vector1.y * vector2.z) - (vector2.y * vector1.z);
64  normal.y = (vector2.x * vector1.z) - (vector1.x * vector2.z);
65  normal.z = (vector1.x * vector2.y) - (vector2.x * vector1.y);
66 
67  return normal;
68 }
69 
70 double length(EDGE_VECTOR &vector)
71 {
72  return sqrt(vector.x * vector.x + vector.y * vector.y + vector.z * vector.z);
73 }
74 
75 PLANE_EDGE_INTERSECTION PlaneEdgeVertexInterpolate(float x1, float y1, float z1, float x2, float y2, float z2, S_V_DATA data1, S_V_DATA data2, int v1, int v2, double p, EDGE_VECTOR unit_normal_vector, double D1, double D2)
76 {
77  double t;
78  double D;
79 
80  PLANE_EDGE_INTERSECTION intersection;
81 
82  /******************/
83  /* Regular Cases */
84  /******************/
85 
86  /* Case 1: Vertex 1 --> Behind; Vertex 2 --> In Front */
87  if (D1 < 0 && D2 > 0)
88  {
89  /* Intersection Point in Slice Plane */
90  D = 0;
91 
92  /* Calculate Parameter "t" for Line Equation */
93  t = (D - x1 * unit_normal_vector.x - y1 * unit_normal_vector.y - z1 * unit_normal_vector.z - p) / (((x2 - x1) * unit_normal_vector.x) + ((y2 - y1) * unit_normal_vector.y) + ((z2 - z1) * unit_normal_vector.z));
94 
95  intersection.intersection_at_vertex1 = false;
96  intersection.intersection_at_vertex2 = false;
97  intersection.int_flag = 1;
98  intersection.vertex1 = v1;
99  intersection.vertex2 = v2;
100 
101  /* Calculate Intersection Coordinates */
102  intersection.intersection.x = x1 + t * (x2 - x1);
103  intersection.intersection.y = y1 + t * (y2 - y1);
104  intersection.intersection.z = z1 + t * (z2 - z1);
105 
106  /* Interpolate Intersection Data Value */
107  intersection.data_vertex_int.dimension = data1.dimension;
108  for (int i = 0; i < intersection.data_vertex_int.dimension; ++i)
109  intersection.data_vertex_int.v[i] = (float)(t * (data2.v[i] - data1.v[i]) + data1.v[i]);
110  }
111 
112  /* Case 2: Vertex 2 --> Behind; Vertex 1 --> In front */
113  if (D1 > 0 && D2 < 0)
114  {
115  /* Intersection Point in Slice Plane */
116  D = 0;
117 
118  /* Calculate Parameter "t" for Line Equation */
119  t = (D - x1 * unit_normal_vector.x - y1 * unit_normal_vector.y - z1 * unit_normal_vector.z - p) / (((x2 - x1) * unit_normal_vector.x) + ((y2 - y1) * unit_normal_vector.y) + ((z2 - z1) * unit_normal_vector.z));
120 
121  intersection.intersection_at_vertex1 = false;
122  intersection.intersection_at_vertex2 = false;
123  intersection.int_flag = 1;
124  intersection.vertex1 = v1;
125  intersection.vertex2 = v2;
126 
127  /* Calculate Intersection Coordinates */
128  intersection.intersection.x = x1 + t * (x2 - x1);
129  intersection.intersection.y = y1 + t * (y2 - y1);
130  intersection.intersection.z = z1 + t * (z2 - z1);
131 
132  /* Interpolate Intersection Data Value */
133  intersection.data_vertex_int.dimension = data1.dimension;
134  for (int i = 0; i < intersection.data_vertex_int.dimension; ++i)
135  intersection.data_vertex_int.v[i] = (float)(t * (data2.v[i] - data1.v[i]) + data1.v[i]);
136  }
137 
138  /***********************/
139  /* Degenerate Cases */
140  /***********************/
141 
142  /* Case 3: Vertex 1 --> Behind; Vertex 2 --> In Slice Plane */
143  if (D1 < 0 && D2 == 0)
144  {
145  /* Calculate Edge-Plane Intersection */
146  intersection.intersection_at_vertex1 = false;
147  intersection.intersection_at_vertex2 = true;
148  intersection.int_flag = 2;
149  intersection.vertex1 = v1;
150  intersection.vertex2 = v2;
151 
152  /* Calculate Intersection Coordinates */
153  intersection.intersection.x = x2;
154  intersection.intersection.y = y2;
155  intersection.intersection.z = z2;
156 
157  /* Interpolate Intersection Data Value */
158  intersection.data_vertex_int = data2;
159  }
160 
161  /* Case 4: Vertex 2 --> Behind; Vertex 1 --> In Slice Plane */
162  if (D1 == 0 && D2 < 0)
163  {
164  /* Calculate Edge-Plane Intersection */
165  intersection.intersection_at_vertex1 = true;
166  intersection.intersection_at_vertex2 = false;
167  intersection.int_flag = 2;
168  intersection.vertex1 = v1;
169  intersection.vertex2 = v2;
170 
171  /* Calculate Intersection Coordinates */
172  intersection.intersection.x = x1;
173  intersection.intersection.y = y1;
174  intersection.intersection.z = z1;
175 
176  /* Interpolate Intersection Data Value */
177  intersection.data_vertex_int = data1;
178  }
179 
180  /* Case 5: Vertex 1 --> In Front; Vertex 2 --> In Slice Plane */
181  if (D1 > 0 && D2 == 0)
182  {
183  /* Calculate Edge-Plane Intersection */
184  intersection.intersection_at_vertex1 = false;
185  intersection.intersection_at_vertex2 = true;
186  intersection.int_flag = 2;
187  intersection.vertex1 = v1;
188  intersection.vertex2 = v2;
189 
190  /* Calculate Intersection Coordinates */
191  intersection.intersection.x = x2;
192  intersection.intersection.y = y2;
193  intersection.intersection.z = z2;
194 
195  /* Interpolate Intersection Data Value */
196  intersection.data_vertex_int = data2;
197  }
198 
199  /* Case 6: Vertex 2 --> In Front; Vertex 1 --> In Slice Plane */
200  if (D1 == 0 && D2 > 0)
201  {
202  /* Calculate Edge-Plane Intersection */
203  intersection.intersection_at_vertex1 = true;
204  intersection.intersection_at_vertex2 = false;
205  intersection.int_flag = 2;
206  intersection.vertex1 = v1;
207  intersection.vertex2 = v2;
208 
209  /* Calculate Intersection Coordinates */
210  intersection.intersection.x = x1;
211  intersection.intersection.y = y1;
212  intersection.intersection.z = z1;
213 
214  /* Interpolate Intersection Data Value */
215  intersection.data_vertex_int = data1;
216  }
217 
218  return (intersection);
219 }
220 
221 bool test_plane_edge_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, PLANE_EDGE_INTERSECTION &intsec, float *x_coord_in, float *y_coord_in, float *z_coord_in)
222 {
223  /**************************************************/
224  /* Test if intersection has already been processed */
225  /**************************************************/
226 
227  for (vector<PLANE_EDGE_INTERSECTION>::iterator existent_intsec = intsec_vector.begin(); existent_intsec < intsec_vector.end(); ++existent_intsec)
228  {
229  if (existent_intsec->intersection_at_vertex1) // Intersection located at vertex 1
230  {
231  if (existent_intsec->vertex1 == intsec.vertex1 || existent_intsec->vertex1 == intsec.vertex2)
232  {
233  return true;
234  }
235  }
236  else if (existent_intsec->intersection_at_vertex2) // Intersection located at vertex 2
237  {
238  if (existent_intsec->vertex1 == intsec.vertex1 || existent_intsec->vertex1 == intsec.vertex2)
239  {
240  return true;
241  }
242  }
243  else // Intersection located between two vertices
244  {
245 
246  // // Check intersection coordinates
247  // if(intsec_vector[i].x == intsec.x && intsec_vector[i].y == intsec.y && intsec_vector[i].z == intsec.z)
248  // {
249  // return true;
250  // }
251  //
252  // // Determine a tolerance to avoid machine precision errors
253  // if(fabs(intsec_vector[i].x - intsec.x) < 0.000005 && fabs(intsec_vector[i].y - intsec.y) < 0.000005 && fabs(intsec_vector[i].z - intsec.z) < 0.000005)
254  // {
255  // return true;
256  // }
257 
258  // Check edge vertices
259  if (existent_intsec->vertex1 == intsec.vertex1 && existent_intsec->vertex2 == intsec.vertex2)
260  {
261  return true;
262  }
263 
264  // Check edge vertices (swapped)
265  if (existent_intsec->vertex1 == intsec.vertex2 && existent_intsec->vertex2 == intsec.vertex1)
266  {
267  return true;
268  }
269 
270  // Check T-vertices
271  if (existent_intsec->vertex1 == intsec.vertex1 || existent_intsec->vertex2 == intsec.vertex2 || existent_intsec->vertex1 == intsec.vertex2 || existent_intsec->vertex2 == intsec.vertex1)
272  {
273  // check direction of edges
274  EDGE_VECTOR d1;
275  d1.x = x_coord_in[existent_intsec->vertex1] - x_coord_in[existent_intsec->vertex2];
276  d1.y = y_coord_in[existent_intsec->vertex1] - y_coord_in[existent_intsec->vertex2];
277  d1.z = z_coord_in[existent_intsec->vertex1] - z_coord_in[existent_intsec->vertex2];
278  EDGE_VECTOR d2;
279  d2.x = x_coord_in[intsec.vertex1] - x_coord_in[intsec.vertex2];
280  d2.y = y_coord_in[intsec.vertex1] - y_coord_in[intsec.vertex2];
281  d2.z = z_coord_in[intsec.vertex1] - z_coord_in[intsec.vertex2];
282  double length1 = length(d1);
283  double length2 = length(d2);
284  if ((length1 > 0) && (length2 > 0))
285  {
286  double cosangle = dot_product(d1, d2) / (length1 * length2);
287  if (fabs(cosangle - 1.0) < 0.0001 || fabs(cosangle + 1.0) < 0.0001)
288  {
289  // ignore if the two edges have the same direction (and share a common vertex)
290  return true;
291  }
292  }
293  }
294  }
295  }
296 
297  return false;
298 }
299 
300 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)
301 {
302  int i;
303  int j;
304 
305  S_V_DATA data_vertex1;
306  S_V_DATA data_vertex2;
307  data_vertex1.dimension = (wdata_in == NULL) ? 1 : 3;
308  data_vertex2.dimension = (wdata_in == NULL) ? 1 : 3;
309 
310  double D1;
311  double D2;
312 
313  int end;
314  int index;
315 
316  EDGE_VECTOR cell_vertex_1;
317  EDGE_VECTOR cell_vertex_2;
318 
319  PLANE_EDGE_INTERSECTION intersec;
320  PLANE_EDGE_INTERSECTION_VECTOR intsec_vector;
321 
322  /* Analyze Polyhedron Faces */
323  for (i = 0; i < num_elem_in; i++)
324  {
325  if (i < num_elem_in - 1)
326  {
327  end = elem_in[i + 1] - 1;
328  }
329  else
330  {
331  end = num_conn_in - 1;
332  }
333 
334  for (j = elem_in[i]; j <= end; j++)
335  {
336  if (j < end)
337  {
338  index = j + 1;
339  }
340  else
341  {
342  index = elem_in[i];
343  }
344 
345  /* Vector to First Cell Vertex */
346  cell_vertex_1.x = x_coord_in[conn_in[j]];
347  cell_vertex_1.y = y_coord_in[conn_in[j]];
348  cell_vertex_1.z = z_coord_in[conn_in[j]];
349 
350  /* Vector to Second Cell Vertex */
351  cell_vertex_2.x = x_coord_in[conn_in[index]];
352  cell_vertex_2.y = y_coord_in[conn_in[index]];
353  cell_vertex_2.z = z_coord_in[conn_in[index]];
354 
355  /* Vertex Data Values */
356  //sdata_in->get_point_value(conn_in[j], &data_vertex1);
357  //sdata_in->get_point_value(conn_in[index], &data_vertex2);
358  data_vertex1.v[0] = udata_in[conn_in[j]];
359  data_vertex2.v[0] = udata_in[conn_in[index]];
360  if (wdata_in != NULL)
361  {
362  data_vertex1.v[1] = vdata_in[conn_in[j]];
363  data_vertex2.v[1] = vdata_in[conn_in[index]];
364  data_vertex1.v[2] = wdata_in[conn_in[j]];
365  data_vertex2.v[2] = wdata_in[conn_in[index]];
366  }
367 
368  /* Calculate Point-Plane Distance (Hessian Normal Form) */
369  D1 = dot_product(unit_normal_vector, cell_vertex_1) + p;
370  D2 = dot_product(unit_normal_vector, cell_vertex_2) + p;
371 
372  /******************/
373  /* Regular Cases */
374  /******************/
375 
376  /* Case 1: Vertex 1 --> Behind; Vertex 2 --> In Front */
377  if (D1 < 0 && D2 > 0)
378  {
379  /* Calculate Edge-Plane Intersection */
380  intersec = PlaneEdgeVertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[index]], y_coord_in[conn_in[index]], z_coord_in[conn_in[index]], data_vertex1, data_vertex2, conn_in[j], conn_in[index], p, unit_normal_vector, D1, D2);
381 
382  /* Discard Repeated Intersections */
383  if (!test_plane_edge_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in))
384  {
385  intsec_vector.push_back(intersec);
386  }
387  }
388 
389  /* Case 2: Vertex 2 --> Behind; Vertex 1 --> In front */
390  if (D1 > 0 && D2 < 0)
391  {
392  /* Calculate Edge-Plane Intersection */
393  intersec = PlaneEdgeVertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[index]], y_coord_in[conn_in[index]], z_coord_in[conn_in[index]], data_vertex1, data_vertex2, conn_in[j], conn_in[index], p, unit_normal_vector, D1, D2);
394 
395  /* Discard Repeated Intersections */
396  if (!test_plane_edge_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in))
397  {
398  intsec_vector.push_back(intersec);
399  }
400  }
401 
402  /**********************/
403  /* Degenerate Cases */
404  /**********************/
405 
406  /* Case 3: Vertex 1 --> Behind; Vertex 2 --> In Slice Plane */
407  if (D1 < 0 && D2 == 0)
408  {
409  /* Calculate Edge-Plane Intersection */
410  intersec = PlaneEdgeVertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[index]], y_coord_in[conn_in[index]], z_coord_in[conn_in[index]], data_vertex1, data_vertex2, conn_in[j], conn_in[index], p, unit_normal_vector, D1, D2);
411 
412  /* Discard Repeated Intersections */
413  if (!test_plane_edge_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in))
414  {
415  intsec_vector.push_back(intersec);
416  }
417  }
418 
419  /* Case 4: Vertex 2 --> Behind; Vertex 1 --> In Slice Plane */
420  if (D1 == 0 && D2 < 0)
421  {
422  /* Calculate Edge-Plane Intersection */
423  intersec = PlaneEdgeVertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[index]], y_coord_in[conn_in[index]], z_coord_in[conn_in[index]], data_vertex1, data_vertex2, conn_in[j], conn_in[index], p, unit_normal_vector, D1, D2);
424 
425  /* Discard Repeated Intersections */
426  if (!test_plane_edge_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in))
427  {
428  intsec_vector.push_back(intersec);
429  }
430  }
431 
432  /* Case 5: Vertex 1 --> In Front; Vertex 2 --> In Slice Plane */
433  if (D1 > 0 && D2 == 0)
434  {
435  /* Calculate Edge-Plane Intersection */
436  intersec = PlaneEdgeVertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[index]], y_coord_in[conn_in[index]], z_coord_in[conn_in[index]], data_vertex1, data_vertex2, conn_in[j], conn_in[index], p, unit_normal_vector, D1, D2);
437 
438  /* Discard Repeated Intersections */
439  if (!test_plane_edge_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in))
440  {
441  intsec_vector.push_back(intersec);
442  }
443  }
444 
445  /* Case 6: Vertex 2 --> In Front; Vertex 1 --> In Slice Plane */
446  if (D1 == 0 && D2 > 0)
447  {
448  /* Calculate Edge-Plane Intersection */
449  intersec = PlaneEdgeVertexInterpolate(x_coord_in[conn_in[j]], y_coord_in[conn_in[j]], z_coord_in[conn_in[j]], x_coord_in[conn_in[index]], y_coord_in[conn_in[index]], z_coord_in[conn_in[index]], data_vertex1, data_vertex2, conn_in[j], conn_in[index], p, unit_normal_vector, D1, D2);
450 
451  /* Discard Repeated Intersections */
452  if (!test_plane_edge_intersection(intsec_vector, intersec, x_coord_in, y_coord_in, z_coord_in))
453  {
454  intsec_vector.push_back(intersec);
455  }
456  }
457 
458  /* Case 7: Vertex 1 --> In Slice Plane; Vertex 2 --> In Slice Plane */
459  if (D1 == 0 && D2 == 0)
460  {
461  /* Discard intersections when both of them are located within the slice plane */
462  }
463  }
464  }
465 
466  return intsec_vector;
467 }
468 
469 int assign_int_index(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
470 {
471  int i;
472  int index;
473 
474  index = 0;
475 
476  for (i = 0; i < intsec_vector.size(); i++)
477  {
478  PLANE_EDGE_INTERSECTION intsec = intsec_vector[i];
479  if (intsec.intersection_at_vertex1 == false && intsec.intersection_at_vertex2 == false)
480  {
481  if ((edge_vertex1 == intsec.vertex1) && (edge_vertex2 == intsec.vertex2))
482  {
483  index = i;
484  }
485 
486  if ((edge_vertex2 == intsec.vertex1) && (edge_vertex1 == intsec.vertex2))
487  {
488  index = i;
489  }
490  }
491 
492  if (intsec_vector[i].intersection_at_vertex1 == true)
493  {
494  if (intsec.vertex1 == edge_vertex1 || intsec.vertex1 == edge_vertex2)
495  {
496  index = i;
497  }
498  }
499 
500  if (intsec.intersection_at_vertex2 == true)
501  {
502  if (intsec.vertex2 == edge_vertex1 || intsec.vertex2 == edge_vertex2)
503  {
504  index = i;
505  }
506  }
507  }
508 
509  return index;
510 }
511 
512 bool find_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
513 {
514  bool int_found;
515 
516  int int_found_flag;
517 
518  int_found_flag = 0;
519  int_found = false;
520 
521  for (vector<PLANE_EDGE_INTERSECTION>::iterator intsec = intsec_vector.begin(); intsec < intsec_vector.end(); ++intsec)
522  {
523  if (int_found_flag == 0)
524  {
525  if (intsec->intersection_at_vertex1 == false && intsec->intersection_at_vertex2 == false)
526  {
527  if ((edge_vertex1 == intsec->vertex1) && (edge_vertex2 == intsec->vertex2))
528  {
529  int_found = true;
530  int_found_flag = 1;
531  }
532 
533  if ((edge_vertex2 == intsec->vertex1) && (edge_vertex1 == intsec->vertex2))
534  {
535  int_found = true;
536  int_found_flag = 1;
537  }
538  }
539 
540  if (intsec->intersection_at_vertex1 == true)
541  {
542  if (intsec->vertex1 == edge_vertex1 || intsec->vertex1 == edge_vertex2)
543  {
544  int_found = true;
545  int_found_flag = 1;
546  }
547  }
548 
549  if (intsec->intersection_at_vertex2 == true)
550  {
551  if (intsec->vertex2 == edge_vertex1 || intsec->vertex2 == edge_vertex2)
552  {
553  int_found = true;
554  int_found_flag = 1;
555  }
556  }
557  }
558  }
559 
560  return int_found;
561 }
562 
563 // void find_current_face(CONTOUR &contour, int edge_vertex1, int edge_vertex2, int *index_list, int *polygon_list, int num_coord_in, int num_conn_in, int &current_face)
564 // {
565 // int i;
566 // int j;
567 // int face_flag;
568 // int copy_current_face;
569 // int neighbour_face1;
570 // int neighbour_face2;
571 //
572 // ITERATOR it;
573 //
574 // face_flag = 0;
575 // copy_current_face = current_face;
576 //
577 // /******************************************************************************/
578 // /* Find the next face of the polyhedron to continue tracing the convex contour */
579 // /******************************************************************************/
580 //
581 // if((edge_vertex1 < num_coord_in - 1) && (edge_vertex2 < num_coord_in - 1))
582 // {
583 // for(i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
584 // {
585 // neighbour_face1 = polygon_list[i];
586 // for(j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
587 // {
588 // neighbour_face2 = polygon_list[j];
589 // if(face_flag == 0)
590 // {
591 // /* Search among the elements that contain both vertices */
592 // if(neighbour_face1 == neighbour_face2)
593 // {
594 // /* Test if the element has already been processed */
595 // it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
596 //
597 // /* The current face that contains the vertices has not been processed */
598 // if(it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
599 // {
600 // contour.polyhedron_faces.push_back(neighbour_face1);
601 // current_face = neighbour_face1;
602 // face_flag = 1;
603 // }
604 // }
605 // }
606 // }
607 // }
608 //
609 // /* All intersection faces have been processed at least once */
610 // if(face_flag == 0)
611 // {
612 // for(i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
613 // {
614 // neighbour_face1 = polygon_list[i];
615 // for(j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
616 // {
617 // neighbour_face2 = polygon_list[j];
618 // if(face_flag == 0)
619 // {
620 // /* Search among the elements that contain both vertices */
621 // if(neighbour_face1 == neighbour_face2)
622 // {
623 // if(neighbour_face1 != copy_current_face)
624 // {
625 // contour.polyhedron_faces.push_back(neighbour_face1);
626 // current_face = neighbour_face1;
627 // face_flag = 1;
628 // }
629 // }
630 // }
631 // }
632 // }
633 // }
634 // }
635 //
636 // else
637 // {
638 // if((edge_vertex1 < num_coord_in - 1) && (edge_vertex2 == num_coord_in - 1))
639 // {
640 // for(i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
641 // {
642 // neighbour_face1 = polygon_list[i];
643 // for(j = index_list[edge_vertex2]; j < num_conn_in; j++)
644 // {
645 // neighbour_face2 = polygon_list[j];
646 // if(face_flag == 0)
647 // {
648 // /* Search among the elements that contain both vertices */
649 // if(neighbour_face1 == neighbour_face2)
650 // {
651 // /* Test if the element has already been processed */
652 // it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
653 //
654 // /* The current face that contains the vertices has not been processed */
655 // if(it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
656 // {
657 // contour.polyhedron_faces.push_back(neighbour_face1);
658 // current_face = neighbour_face1;
659 // face_flag = 1;
660 // }
661 // }
662 // }
663 // }
664 // }
665 //
666 // /* All intersection faces have been processed at least once */
667 // if(face_flag == 0)
668 // {
669 // for(i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
670 // {
671 // neighbour_face1 = polygon_list[i];
672 // for(j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
673 // {
674 // neighbour_face2 = polygon_list[j];
675 // if(face_flag == 0)
676 // {
677 // /* Search among the elements that contain both vertices */
678 // if(neighbour_face1 == neighbour_face2)
679 // {
680 // if(neighbour_face1 != copy_current_face)
681 // {
682 // contour.polyhedron_faces.push_back(neighbour_face1);
683 // current_face = neighbour_face1;
684 // face_flag = 1;
685 // }
686 // }
687 // }
688 // }
689 // }
690 // }
691 // }
692 //
693 // if((edge_vertex1 == num_coord_in - 1) && (edge_vertex2 < num_coord_in - 1))
694 // {
695 // for(i = index_list[edge_vertex1]; i < num_conn_in; i++)
696 // {
697 // neighbour_face1 = polygon_list[i];
698 // for(j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
699 // {
700 // neighbour_face2 = polygon_list[j];
701 // if(face_flag == 0)
702 // {
703 // /* Search among the elements that contain both vertices */
704 // if(neighbour_face1 == neighbour_face2)
705 // {
706 // /* Test if the element has already been processed */
707 // it = find(contour.polyhedron_faces.begin(), contour.polyhedron_faces.end(), neighbour_face1);
708 //
709 // /* The current face that contains the vertices has not been processed */
710 // if(it == contour.polyhedron_faces.end() || contour.polyhedron_faces.size() == 0)
711 // {
712 // contour.polyhedron_faces.push_back(neighbour_face1);
713 // current_face = neighbour_face1;
714 // face_flag = 1;
715 // }
716 // }
717 // }
718 // }
719 // }
720 //
721 // /* All intersection faces have been processed at least once */
722 // if(face_flag == 0)
723 // {
724 // for(i = index_list[edge_vertex1]; i < index_list[edge_vertex1 + 1]; i++)
725 // {
726 // neighbour_face1 = polygon_list[i];
727 // for(j = index_list[edge_vertex2]; j < index_list[edge_vertex2 + 1]; j++)
728 // {
729 // neighbour_face2 = polygon_list[j];
730 // if(face_flag == 0)
731 // {
732 // /* Search among the elements that contain both vertices */
733 // if(neighbour_face1 == neighbour_face2)
734 // {
735 // if(neighbour_face1 != copy_current_face)
736 // {
737 // contour.polyhedron_faces.push_back(neighbour_face1);
738 // current_face = neighbour_face1;
739 // face_flag = 1;
740 // }
741 // }
742 // }
743 // }
744 // }
745 // }
746 // }
747 // }
748 // }
749 
750 // bool test_angle(float component1a, float component2a, float component1b, float component2b)
751 // {
752 // float cross_prod;
753 // float distance1;
754 // float distance2;
755 //
756 // cross_prod = component1a*component2b - component1b*component2a;
757 //
758 // if(cross_prod == 0)
759 // {
760 // /* Calculate Manhattan-Distances */
761 // distance1 = fabs(component1a) + fabs(component2a);
762 // distance2 = fabs(component1b) + fabs(component2b);
763 // }
764 //
765 // return cross_prod > 0 || cross_prod == 0 && distance1 > distance2;
766 // }
767 //
768 //
769 // void intersection_quicksortYZ(PLANE_EDGE_INTERSECTION_VECTOR &new_intsec_vector, vector<int> &vertex_polynome, int start, int end)
770 // {
771 // int i;
772 // int j;
773 //
774 // PLANE_EDGE_INTERSECTION quicksort_pivot;
775 //
776 // /* Start Quicksort */
777 // i = start;
778 // j = end;
779 //
780 // quicksort_pivot = new_intsec_vector[(start + end)/2];
781 //
782 // while(i <= j)
783 // {
784 //
785 // while(test_angle(new_intsec_vector[i].y, new_intsec_vector[i].z, quicksort_pivot.y, quicksort_pivot.z))
786 // {
787 // i++;
788 // }
789 //
790 // while(test_angle(quicksort_pivot.y, quicksort_pivot.z, new_intsec_vector[j].y, new_intsec_vector[j].z))
791 // {
792 // j--;
793 // }
794 //
795 // if(i <= j)
796 // {
797 // std::swap(vertex_polynome[i++], vertex_polynome[j--]);
798 // std::swap(new_intsec_vector[i++], new_intsec_vector[j--]);
799 // }
800 // }
801 //
802 // if (start < j)
803 // {
804 // intersection_quicksortYZ(new_intsec_vector, vertex_polynome, start, j);
805 // }
806 //
807 // if (i < end)
808 // {
809 // intersection_quicksortYZ(new_intsec_vector, vertex_polynome, i, end);
810 // }
811 // }
812 //
813 //
814 //
815 // void generate_capping_contour(CONTOUR &capping_contour, PLANE_EDGE_INTERSECTION_VECTOR intsec_vector, EDGE_VECTOR unit_normal_vector, vector<float> &data_vector)
816 // {
817 // int i;
818 //
819 // int pivot_index;
820 //
821 //
822 // vector<int> vertex_polynome;
823 //
824 //
825 //
826 //
827 //
828 // PLANE_EDGE_INTERSECTION_VECTOR new_intsec_vector;
829 //
830 //
831 //
832 //
833 //
834 // /*******************************************************/
835 // /* Case 1 --> Intersections form a point, line, or triangle */
836 // /*******************************************************/
837 //
838 // if(intsec_vector.size() <= 3)
839 // {
840 // for(i = 0; i < intsec_vector.size(); i++)
841 // {
842 // capping_contour.ring.push_back(i);
843 // }
844 //
845 // capping_contour.ring_index.push_back(0);
846 // }
847 //
848 // /*****************************************************/
849 // /* Case 2 --> Intersections form an arbitrary polygon */
850 // /*****************************************************/
851 //
852 // if(intsec_vector.size() > 3)
853 // {
854 // pivot_index = 0;
855 // vertex_polynome.reserve(15);
856 // new_intsec_vector.reserve(15);
857 //
858 // for(i = 0; i < intsec_vector.size(); i++)
859 // {
860 // vertex_polynome.push_back(i);
861 // new_intsec_vector.push_back(intsec_vector[i]);
862 //
863 // }
864 //
865 // /*********************************/
866 // /* Find Edge of Intersection Hull */
867 // /*********************************/
868 //
869 // /* General Case */
870 // if(unit_normal_vector.z != 0)
871 // {
872 // /* G*/
873 // }
874 //
875 // /* Special Cases */
876 // if(unit_normal_vector.z == 0)
877 // {
878 // /* Special Case 1 --> Sampling Plane Parallel to YZ-Axis */
879 // if(unit_normal_vector.x != 0 && unit_normal_vector.y == 0)
880 // {
881 // for(i = 1; i < intsec_vector.size(); i++)
882 // {
883 // /* Search minimum Z-Coordinate */
884 // if(intsec_vector[pivot_index].z > intsec_vector[i].z || intsec_vector[pivot_index].z == intsec_vector[i].z && intsec_vector[pivot_index].y > intsec_vector[i].y)
885 // {
886 // pivot_index = i;
887 // }
888 // }
889 // }
890 //
891 // /* Substract Pivot */
892 // for(i = 0; i < new_intsec_vector.size(); i++)
893 // {
894 // new_intsec_vector[i].x = new_intsec_vector[i].x - new_intsec_vector[pivot_index].x;
895 // new_intsec_vector[i].y = new_intsec_vector[i].y - new_intsec_vector[pivot_index].y;
896 // new_intsec_vector[i].z = new_intsec_vector[i].z - new_intsec_vector[pivot_index].z;
897 // }
898 //
899 // std::swap(vertex_polynome[0], vertex_polynome[pivot_index]);
900 // std::swap(new_intsec_vector[0], new_intsec_vector[pivot_index]);
901 //
902 // intersection_quicksortYZ(new_intsec_vector, vertex_polynome, 1, new_intsec_vector.size() - 1);
903 //
904 //
905 // // capping_contour.polyhedron_faces.erase(capping_contour.polyhedron_faces.begin(), capping_contour.polyhedron_faces.end());
906 //
907 // for(i = 0; i < vertex_polynome.size(); i++)
908 // {
909 // capping_contour.ring.push_back(vertex_polynome[i]);
910 // }
911 //
912 // capping_contour.ring_index.push_back(0);
913 //
914 // // /* Special Case 2 --> Sampling Plane Parallel to ZX-Plane */
915 // // else if(unit_normal_vector.x == 0 && unit_normal_vector.y != 0)
916 // // {
917 // // }
918 // //
919 // // /* Special Case 3 --> Sampling Plane rotates along the Z-Axis */
920 // // else if(unit_normal_vector.x != 0 && unit_normal_vector.y != 0)
921 // // {
922 // // }
923 //
924 // }
925 // }
926 //
927 // /* Finish Data Output Vector */
928 // /* Note: data vector must not be ordered according to the Graham Scan */
929 // for(i = 0; i < intsec_vector.size(); i++)
930 // {
931 // data_vector.push_back(intsec_vector[i].data_vertex_int);
932 // }
933 // }
934 
935 // void generate_capping_contour(CONTOUR &capping_contour, PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int *index_list, int *polygon_list, int num_coord_in, int num_conn_in, int num_elem_in, int *elem_in, int *conn_in, vector<float> &data_vector)
936 // {
937 // bool extreme_edge_found;
938 //
939 // //char info[256];
940 //
941 // int i;
942 // int intersection_pointer1;
943 // int intersection_pointer2;
944 // int current_face;
945 // int previous_face;
946 // int edge_vertex1;
947 // int edge_vertex2;
948 // int new_edge_vertex1;
949 // int new_edge_vertex2;
950 // int index_v1;
951 // int index_v2;
952 // int index_flag_v1;
953 // int index_flag_v2;
954 // int counter1;
955 // int counter2;
956 // int counter3;
957 // int new_int_found;
958 // int loop_counter;
959 //
960 // double extreme_edge_magnitude;
961 // double base_edge_magnitude;
962 // double cosine;
963 // double lastCosine;
964 //
965 // EDGE_VECTOR extreme_edge;
966 // EDGE_VECTOR base_edge;
967 // EDGE_VECTOR test_vector1;
968 // EDGE_VECTOR test_vector2;
969 //
970 // vector<int> int_index_vector;
971 // vector<double> cosine_vector;
972 // vector<double> difference_vector;
973 //
974 // int smaller_position;
975 // bool angleChanged;
976 // double distance1;
977 // double distance2;
978 //
979 //
980 // /*******************************************************/
981 // /* Case 1 --> Intersections form a point, line, or triangle */
982 // /*******************************************************/
983 //
984 // if(intsec_vector.size() <= 3)
985 // {
986 // for(i = 0; i < intsec_vector.size(); i++)
987 // {
988 // capping_contour.ring.push_back(i);
989 // }
990 //
991 // capping_contour.ring_index.push_back(0);
992 // }
993 //
994 // /*****************************************************/
995 // /* Case 2 --> Intersections form an arbitrary polygon */
996 // /*****************************************************/
997 //
998 // if(intsec_vector.size() > 3)
999 // {
1000 // /**********************/
1001 // /* Find Extreme Edge */
1002 // /**********************/
1003 //
1004 // /* First Intersection as Pivot */
1005 // intersection_pointer1 = 0;
1006 // edge_vertex1 = intsec_vector[0].vertex1;
1007 // edge_vertex2 = intsec_vector[0].vertex2;
1008 // loop_counter = 0;
1009 // counter3 = 0;
1010 // previous_face = -1;
1011 // extreme_edge_found = false;
1012 //
1013 // do
1014 // {
1015 // do
1016 // {
1017 // if(loop_counter == 1)
1018 // {
1019 // previous_face = current_face;
1020 // }
1021 //
1022 // loop_counter++;
1023 //
1024 // /* Find Current Face */
1025 // find_current_face(capping_contour, edge_vertex1, edge_vertex2, index_list, polygon_list, num_coord_in, num_conn_in, current_face);
1026 //
1027 // /*****************************/
1028 // /* Test Vertex Configuration */
1029 // /*****************************/
1030 //
1031 // index_flag_v1 = 0;
1032 // index_flag_v2 = 0;
1033 // new_int_found = 0;
1034 // index_v1 = elem_in[current_face];
1035 // index_v2 = elem_in[current_face];
1036 //
1037 // /***********************************************************************/
1038 // /* Locate index values of the edge vertices in the array of connectivities */
1039 // /***********************************************************************/
1040 //
1041 // if(current_face < num_elem_in - 1)
1042 // {
1043 // for(i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
1044 // {
1045 // if(index_flag_v1 == 0)
1046 // {
1047 // if(edge_vertex1 != conn_in[index_v1])
1048 // {
1049 // index_v1++;
1050 // }
1051 //
1052 // else
1053 // {
1054 // index_flag_v1 = 1;
1055 // }
1056 // }
1057 //
1058 // if(index_flag_v2 == 0)
1059 // {
1060 // if(edge_vertex2 != conn_in[index_v2])
1061 // {
1062 // index_v2++;
1063 // }
1064 //
1065 // else
1066 // {
1067 // index_flag_v2 = 1;
1068 // }
1069 // }
1070 // }
1071 // }
1072 //
1073 // else
1074 // {
1075 // for(i = elem_in[current_face]; i < num_conn_in; i++)
1076 // {
1077 // if(index_flag_v1 == 0)
1078 // {
1079 // if(edge_vertex1 != conn_in[index_v1])
1080 // {
1081 // index_v1++;
1082 // }
1083 //
1084 // else
1085 // {
1086 // index_flag_v1 = 1;
1087 // }
1088 // }
1089 //
1090 // if(index_flag_v2 == 0)
1091 // {
1092 // if(edge_vertex2 != conn_in[index_v2])
1093 // {
1094 // index_v2++;
1095 // }
1096 //
1097 // else
1098 // {
1099 // index_flag_v2 = 1;
1100 // }
1101 // }
1102 // }
1103 // }
1104 //
1105 // counter1 = index_v1;
1106 // counter2 = index_v2;
1107 //
1108 // /***************************/
1109 // /* Locate Next Intersection */
1110 // /***************************/
1111 //
1112 // if(current_face < num_elem_in - 1)
1113 // {
1114 // /**************************/
1115 // /* Search Clockwise (CW) */
1116 // /**************************/
1117 //
1118 // if(index_v1 > index_v2)
1119 // {
1120 // if((index_v1 != elem_in[current_face + 1] - 1) || (index_v2 != elem_in[current_face]))
1121 // {
1122 // for(i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
1123 // {
1124 // counter1++;
1125 // counter2++;
1126 //
1127 // if(counter1 == elem_in[current_face + 1])
1128 // {
1129 // counter1 = elem_in[current_face];
1130 // }
1131 //
1132 // if(counter2 == elem_in[current_face + 1])
1133 // {
1134 // counter2 = elem_in[current_face];
1135 // }
1136 //
1137 // if(new_int_found == 0)
1138 // {
1139 // /* Test if the a new intersection has been encountered */
1140 // new_edge_vertex1 = conn_in[counter1];
1141 // new_edge_vertex2 = conn_in[counter2];
1142 //
1143 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1144 // {
1145 // new_int_found = 1;
1146 // }
1147 // }
1148 // }
1149 // }
1150 //
1151 // if((index_v1 == elem_in[current_face + 1] - 1) && (index_v2 == elem_in[current_face]))
1152 // {
1153 // for(i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
1154 // {
1155 // counter1--;
1156 // counter2--;
1157 //
1158 // if(counter1 == elem_in[current_face] - 1)
1159 // {
1160 // counter1 = elem_in[current_face + 1] - 1;
1161 // }
1162 //
1163 // if(counter2 == elem_in[current_face] - 1)
1164 // {
1165 // counter2 = elem_in[current_face + 1] - 1;
1166 // }
1167 //
1168 // if(new_int_found == 0)
1169 // {
1170 // /* Test if the a new intersection has been encountered */
1171 // new_edge_vertex1 = conn_in[counter1];
1172 // new_edge_vertex2 = conn_in[counter2];
1173 //
1174 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1175 // {
1176 // new_int_found = 1;
1177 // }
1178 // }
1179 // }
1180 // }
1181 // }
1182 //
1183 // /************************************/
1184 // /* Search Counterclockwise (CCW) */
1185 // /************************************/
1186 //
1187 // if(index_v1 < index_v2)
1188 // {
1189 // if((index_v1 != elem_in[current_face]) || (index_v2 != elem_in[current_face +1] - 1))
1190 // {
1191 // for(i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
1192 // {
1193 // counter1--;
1194 // counter2--;
1195 //
1196 // if(counter1 == elem_in[current_face] - 1)
1197 // {
1198 // counter1 = elem_in[current_face + 1] - 1;
1199 // }
1200 //
1201 // if(counter2 == elem_in[current_face] - 1)
1202 // {
1203 // counter2 = elem_in[current_face + 1] - 1;
1204 // }
1205 //
1206 // if(new_int_found == 0)
1207 // {
1208 // /* Test if the a new intersection has been encountered */
1209 // new_edge_vertex1 = conn_in[counter1];
1210 // new_edge_vertex2 = conn_in[counter2];
1211 //
1212 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1213 // {
1214 // new_int_found = 1;
1215 // }
1216 // }
1217 // }
1218 // }
1219 //
1220 // if((index_v1 == elem_in[current_face]) && (index_v2 == elem_in[current_face +1] - 1))
1221 // {
1222 // for(i = elem_in[current_face]; i < elem_in[current_face + 1]; i++)
1223 // {
1224 // counter1++;
1225 // counter2++;
1226 //
1227 // if(counter1 == elem_in[current_face + 1])
1228 // {
1229 // counter1 = elem_in[current_face];
1230 // }
1231 //
1232 // if(counter2 == elem_in[current_face + 1])
1233 // {
1234 // counter2 = elem_in[current_face];
1235 // }
1236 //
1237 // if(new_int_found == 0)
1238 // {
1239 // /* Test if the a new intersection has been encountered */
1240 // new_edge_vertex1 = conn_in[counter1];
1241 // new_edge_vertex2 = conn_in[counter2];
1242 //
1243 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1244 // {
1245 // new_int_found = 1;
1246 // }
1247 // }
1248 // }
1249 // }
1250 // }
1251 // }
1252 //
1253 // else
1254 // {
1255 // /**************************/
1256 // /* Search Clockwise (CW) */
1257 // /**************************/
1258 //
1259 // if(index_v1 > index_v2)
1260 // {
1261 // if((index_v1 != num_conn_in - 1) || (index_v2 != elem_in[current_face]))
1262 // {
1263 // for(i = elem_in[current_face]; i < num_conn_in; i++)
1264 // {
1265 // counter1++;
1266 // counter2++;
1267 //
1268 // if(counter1 == num_conn_in)
1269 // {
1270 // counter1 = elem_in[current_face];
1271 // }
1272 //
1273 // if(counter2 == num_conn_in)
1274 // {
1275 // counter2 = elem_in[current_face];
1276 // }
1277 //
1278 // if(new_int_found == 0)
1279 // {
1280 // /* Test if the a new intersection has been encountered */
1281 // new_edge_vertex1 = conn_in[counter1];
1282 // new_edge_vertex2 = conn_in[counter2];
1283 //
1284 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1285 // {
1286 // new_int_found = 1;
1287 // }
1288 // }
1289 // }
1290 // }
1291 //
1292 // if((index_v1 == num_conn_in - 1) && (index_v2 == elem_in[current_face]))
1293 // {
1294 // for(i = elem_in[current_face]; i < num_conn_in; i++)
1295 // {
1296 // counter1--;
1297 // counter2--;
1298 //
1299 // if(counter1 == elem_in[current_face] - 1)
1300 // {
1301 // counter1 = num_conn_in - 1;
1302 // }
1303 //
1304 // if(counter2 == elem_in[current_face] - 1)
1305 // {
1306 // counter2 = num_conn_in - 1;
1307 // }
1308 //
1309 // if(new_int_found == 0)
1310 // {
1311 // /* Test if the a new intersection has been encountered */
1312 // new_edge_vertex1 = conn_in[counter1];
1313 // new_edge_vertex2 = conn_in[counter2];
1314 //
1315 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1316 // {
1317 // new_int_found = 1;
1318 // }
1319 // }
1320 // }
1321 // }
1322 // }
1323 //
1324 // /************************************/
1325 // /* Search Counterclockwise (CCW) */
1326 // /************************************/
1327 //
1328 // if(index_v1 < index_v2)
1329 // {
1330 // if((index_v1 != elem_in[current_face]) || (index_v2 != num_conn_in - 1))
1331 // {
1332 // for(i = elem_in[current_face]; i < num_conn_in; i++)
1333 // {
1334 // counter1--;
1335 // counter2--;
1336 //
1337 // if(counter1 == elem_in[current_face] - 1)
1338 // {
1339 // counter1 = num_conn_in - 1;
1340 // }
1341 //
1342 // if(counter2 == elem_in[current_face] - 1)
1343 // {
1344 // counter2 = num_conn_in - 1;
1345 // }
1346 //
1347 // if(new_int_found == 0)
1348 // {
1349 // /* Test if the a new intersection has been encountered */
1350 // new_edge_vertex1 = conn_in[counter1];
1351 // new_edge_vertex2 = conn_in[counter2];
1352 //
1353 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1354 // {
1355 // new_int_found = 1;
1356 // }
1357 // }
1358 // }
1359 // }
1360 //
1361 // if((index_v1 == elem_in[current_face]) && (index_v2 == num_conn_in - 1))
1362 // {
1363 // for(i = elem_in[current_face]; i < num_conn_in; i++)
1364 // {
1365 // counter1++;
1366 // counter2++;
1367 //
1368 // if(counter1 == num_conn_in)
1369 // {
1370 // counter1 = elem_in[current_face];
1371 // }
1372 //
1373 // if(counter2 == num_conn_in)
1374 // {
1375 // counter2 = elem_in[current_face];
1376 // }
1377 //
1378 // if(new_int_found == 0)
1379 // {
1380 // /* Test if the a new intersection has been encountered */
1381 // new_edge_vertex1 = conn_in[counter1];
1382 // new_edge_vertex2 = conn_in[counter2];
1383 //
1384 // if(find_intersection(intsec_vector, new_edge_vertex1, new_edge_vertex2))
1385 // {
1386 // new_int_found = 1;
1387 // }
1388 // }
1389 // }
1390 // }
1391 // }
1392 // }
1393 //
1394 // intersection_pointer2 = assign_int_index(intsec_vector, new_edge_vertex1, new_edge_vertex2);
1395 // }
1396 // while(intersection_pointer1 == intersection_pointer2 && previous_face != current_face);
1397 //
1398 // /* Extreme Edge Found */
1399 // if(intersection_pointer1 != intersection_pointer2)
1400 // {
1401 // extreme_edge_found = true;
1402 // }
1403 //
1404 // /* Extreme Edge Not Found --> Select Different Pivot */
1405 // if(!extreme_edge_found)
1406 // {
1407 // counter3++;
1408 // intersection_pointer1 = counter3;
1409 // edge_vertex1 = intsec_vector[counter3].vertex1;
1410 // edge_vertex2 = intsec_vector[counter3].vertex2;
1411 // capping_contour.polyhedron_faces.erase(capping_contour.polyhedron_faces.begin(), capping_contour.polyhedron_faces.end());
1412 // previous_face = -1;
1413 // loop_counter = 0;
1414 // }
1415 // }
1416 // while(!extreme_edge_found && counter3 < intsec_vector.size());
1417 //
1418 // /* Define Extreme Edge Vector */
1419 // extreme_edge.x = intsec_vector[intersection_pointer2].x - intsec_vector[intersection_pointer1].x;
1420 // extreme_edge.y = intsec_vector[intersection_pointer2].y - intsec_vector[intersection_pointer1].y;
1421 // extreme_edge.z = intsec_vector[intersection_pointer2].z - intsec_vector[intersection_pointer1].z;
1422 // /* Define Extreme Edge Magnitude */
1423 // extreme_edge_magnitude = length(extreme_edge);
1424 //
1425 // /*****************/
1426 // /* Graham Scan */
1427 // /*****************/
1428 //
1429 // /* Store first two intersections of the convex hull */
1430 // capping_contour.ring.push_back(intersection_pointer1);
1431 // capping_contour.ring.push_back(intersection_pointer2);
1432 // capping_contour.ring_index.push_back(0);
1433 //
1434 // /* Initiate Scan */
1435 // for(i = 0; i < intsec_vector.size(); i++)
1436 // {
1437 // if(i != intersection_pointer1 && i != intersection_pointer2)
1438 // {
1439 // /* Scan Intersection Points */
1440 // base_edge.x = intsec_vector[i].x - intsec_vector[intersection_pointer1].x;
1441 // base_edge.y = intsec_vector[i].y - intsec_vector[intersection_pointer1].y;
1442 // base_edge.z = intsec_vector[i].z - intsec_vector[intersection_pointer1].z;
1443 // /* Define Vector Magnitude */
1444 // base_edge_magnitude = length(base_edge);
1445 //
1446 // if (extreme_edge_magnitude > 0 && base_edge_magnitude > 0)
1447 // {
1448 // /* Define Angle between Vectors */
1449 // cosine = dot_product(base_edge, extreme_edge)/(extreme_edge_magnitude*base_edge_magnitude);
1450 // }
1451 // else
1452 // {
1453 // cosine = 0.0;
1454 // }
1455 //
1456 // cosine_vector.push_back(cosine);
1457 //
1458 // /* Store Intersection Indices */
1459 // int_index_vector.push_back(i);
1460 // }
1461 // }
1462 //
1463 // /*******************/
1464 // /* Debugging Info */
1465 // /*******************/
1466 //
1467 // //sprintf(info, "Angle Vector Size = %d", cosine_vector.size());
1468 // //Covise::send_info(info);
1469 // //sprintf(info, "Intersection Index Vector Size = %d", int_index_vector.size());
1470 // //Covise::send_info(info);
1471 //
1472 // // calculate angle of the first two vertices in the ring
1473 // base_edge.x = intsec_vector[intersection_pointer2].x - intsec_vector[intersection_pointer1].x;
1474 // base_edge.y = intsec_vector[intersection_pointer2].y - intsec_vector[intersection_pointer1].y;
1475 // base_edge.z = intsec_vector[intersection_pointer2].z - intsec_vector[intersection_pointer1].z;
1476 // base_edge_magnitude = length(base_edge);
1477 // if (extreme_edge_magnitude > 0 && base_edge_magnitude > 0)
1478 // {
1479 // lastCosine = dot_product(base_edge, extreme_edge)/(extreme_edge_magnitude*base_edge_magnitude);
1480 // }
1481 // else
1482 // {
1483 // lastCosine = 0.0;
1484 // }
1485 //
1486 // angleChanged = false;
1487 // const double ANGLE_THRESHHOLD = 0.0001; // higher values may display slightly ?konkave? cells correctly
1488 // do
1489 // {
1490 //
1491 // smaller_position = 0;
1492 // for(i = 1; i < cosine_vector.size(); i++)
1493 // {
1494 // if(fabs(cosine_vector[smaller_position]-cosine_vector[i]) < ANGLE_THRESHHOLD)
1495 // {
1496 // test_vector1.x = intsec_vector[int_index_vector[smaller_position]].x - intsec_vector[intersection_pointer1].x;
1497 // test_vector1.y = intsec_vector[int_index_vector[smaller_position]].y - intsec_vector[intersection_pointer1].y;
1498 // test_vector1.z = intsec_vector[int_index_vector[smaller_position]].z - intsec_vector[intersection_pointer1].z;
1499 //
1500 // test_vector2.x = intsec_vector[int_index_vector[i]].x - intsec_vector[intersection_pointer1].x;
1501 // test_vector2.y = intsec_vector[int_index_vector[i]].y - intsec_vector[intersection_pointer1].y;
1502 // test_vector2.z = intsec_vector[int_index_vector[i]].z - intsec_vector[intersection_pointer1].z;
1503 //
1504 // distance1 = length(test_vector1);
1505 // distance2 = length(test_vector2);
1506 //
1507 // if((angleChanged && (distance1 < distance2)) || (!angleChanged && (distance1 > distance2)))
1508 // {
1509 // // as long as the contour goes in a straight line away from intersection_pointer1,
1510 // // prefer points closer to intersection_pointer1
1511 // smaller_position = i;
1512 // }
1513 // }
1514 // else if(cosine_vector[smaller_position] < cosine_vector[i]) // "cosine1 < cosine2" means "angle1 > angle2"
1515 // {
1516 // smaller_position = i;
1517 // }
1518 // }
1519 //
1520 // capping_contour.ring.push_back(int_index_vector[smaller_position]);
1521 // angleChanged = angleChanged || (fabs(lastCosine-cosine_vector[smaller_position]) >= ANGLE_THRESHHOLD);
1522 // lastCosine = cosine_vector[smaller_position];
1523 //
1524 // // erase
1525 // int_index_vector.erase(int_index_vector.begin() + smaller_position);
1526 // cosine_vector.erase(cosine_vector.begin() + smaller_position);
1527 // }
1528 // while(capping_contour.ring.size() < intsec_vector.size());
1529 // }
1530 //
1531 // /* Finish Data Output Vector */
1532 // /* Note: data vector must not be ordered according to the Graham Scan */
1533 // for(i = 0; i < intsec_vector.size(); i++)
1534 // {
1535 // data_vector.push_back(intsec_vector[i].data_vertex_int);
1536 // }
1537 //
1538 // /*******************/
1539 // /* Debugging Info */
1540 // /*******************/
1541 //
1542 // //sprintf(info, "Ring Size = %d", capping_contour.ring.size());
1543 // //Covise::send_info(info);
1544 // //sprintf(info, "Ring Index Size = %d", capping_contour.ring_index.size());
1545 // //Covise::send_info(info);
1546 // //sprintf(info, "Convex Hull of Intersections:");
1547 // //Covise::send_info(info);
1548 //
1549 // //for(i = 0; i < capping_contour.ring.size(); i++)
1550 // //{
1551 // // sprintf(info, "Intersection [%d ] = %d", i, capping_contour.ring[i]);
1552 // // Covise::send_info(info);
1553 // //}
1554 //
1555 // //for(i = 0; i < data_vector.size(); i++)
1556 // //{
1557 // // sprintf(info, "Intersection Data Value [%d ] = %f", i, data_vector[i]);
1558 // // Covise::send_info(info);
1559 // //}
1560 // }
1561 
1562 // Explanation of cosine:
1563 // 0 if vector has the same direction as plane_base_x
1564 // 1 if vector has the same direction as plane_base_y
1565 // 2 opposite to plane_base_x
1566 // 3 opposite to normal plane_base_y
1567 // Note:
1568 // plane_base_x has to be normalized
1569 double get_cosine(EDGE_VECTOR &vector, EDGE_VECTOR &plane_base_x, EDGE_VECTOR &plane_base_y)
1570 {
1571  double len = length(vector);
1572  if (len > 0.0)
1573  {
1574  double cosine = dot_product(vector, plane_base_x) / len;
1575  if (dot_product(vector, plane_base_y) > 0.0)
1576  {
1577  return 1.0 - cosine;
1578  }
1579  else
1580  {
1581  return 3.0 + cosine;
1582  }
1583  }
1584  else
1585  {
1586  return 0.0;
1587  }
1588 }
1589 
1590 void generate_capping_contour(CONTOUR &capping_contour, PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, EDGE_VECTOR &plane_base_x, EDGE_VECTOR &plane_base_y, vector<float> &u_data_vector, vector<float> &v_data_vector, vector<float> &w_data_vector)
1591 {
1592  int size = (int)intsec_vector.size(); // we need this quite often
1593  if (size <= 3)
1594  {
1595  for (int i = 0; i < size; ++i)
1596  {
1597  capping_contour.ring.push_back(i);
1598  }
1599  capping_contour.ring_index.push_back(0);
1600  }
1601  else
1602  {
1603 
1604  // calculate the average position of all vertices
1605  EDGE_VECTOR average;
1606  average.x = 0.0;
1607  average.y = 0.0;
1608  average.z = 0.0;
1609  for (vector<PLANE_EDGE_INTERSECTION>::iterator intsec = intsec_vector.begin(); intsec < intsec_vector.end(); ++intsec)
1610  {
1611  average.x += intsec->intersection.x;
1612  average.y += intsec->intersection.y;
1613  average.z += intsec->intersection.z;
1614  }
1615  average.x /= double(size);
1616  average.y /= double(size);
1617  average.z /= double(size);
1618 
1619  // initiate scan (calculate angles)
1620  vector<int> vertex_index_vector;
1621  vector<double> vertex_cosine_vector;
1622  for (int i = 0; i < size; ++i)
1623  {
1624  EDGE_VECTOR current_vector = intsec_vector[i].intersection;
1625  current_vector.x -= average.x;
1626  current_vector.y -= average.y;
1627  current_vector.z -= average.z;
1628  vertex_cosine_vector.push_back(get_cosine(current_vector, plane_base_x, plane_base_y));
1629  vertex_index_vector.push_back(i);
1630  }
1631 
1632  // collect all vertices in ascending order (angles)
1633  capping_contour.ring_index.push_back(0);
1634  do
1635  {
1636  int smaller_position = 0;
1637  for (int i = 1; i < vertex_cosine_vector.size(); ++i)
1638  {
1639  if (vertex_cosine_vector[i] < vertex_cosine_vector[smaller_position])
1640  {
1641  smaller_position = i;
1642  }
1643  }
1644 
1645  capping_contour.ring.push_back(vertex_index_vector[smaller_position]);
1646 
1647  // erase
1648  vertex_index_vector.erase(vertex_index_vector.begin() + smaller_position);
1649  vertex_cosine_vector.erase(vertex_cosine_vector.begin() + smaller_position);
1650  } while (capping_contour.ring.size() < size);
1651 
1652  // test if polygon is concave
1653  double lastCosine = 0.;
1654  bool concave = false;
1655  int concaveVertex;
1656  for (int i = 0; i <= capping_contour.ring.size(); ++i)
1657  {
1658  int vertex_index1 = capping_contour.ring[concaveVertex = (i) % capping_contour.ring.size()];
1659  int vertex_index2 = capping_contour.ring[(i + 1) % capping_contour.ring.size()];
1660  EDGE_VECTOR intersection1 = intsec_vector[vertex_index1].intersection;
1661  EDGE_VECTOR intersection2 = intsec_vector[vertex_index2].intersection;
1662  EDGE_VECTOR current_edge;
1663  current_edge.x = intersection2.x - intersection1.x;
1664  current_edge.y = intersection2.y - intersection1.y;
1665  current_edge.z = intersection2.z - intersection1.z;
1666  double cosine = get_cosine(current_edge, plane_base_x, plane_base_y);
1667  if (i != 0)
1668  {
1669  // check
1670  double diff = cosine - lastCosine;
1671  if (diff < 0.0)
1672  diff = 4.0 + diff;
1673  if (diff > 2.0)
1674  {
1675  concave = true;
1676  break;
1677  }
1678  }
1679  lastCosine = cosine;
1680  }
1681 
1682  if (concave)
1683  {
1684  // we don't want concave polygons -> triangulate polygon
1685  tr_vertexVector vertices2d;
1686  for (vector<int>::iterator index = capping_contour.ring.begin(); index < capping_contour.ring.end(); ++index)
1687  {
1688  EDGE_VECTOR intersection = intsec_vector[*index].intersection;
1689  // go to 2D
1690  tr_vertex vertex2D;
1691  vertex2D.x = (float)dot_product(intersection, plane_base_x);
1692  vertex2D.y = (float)dot_product(intersection, plane_base_y);
1693  vertex2D.index = *index;
1694  // add
1695  vertices2d.push_back(vertex2D);
1696  }
1698  Triangulate::Process(vertices2d, result);
1699  capping_contour.ring.clear();
1700  capping_contour.ring_index.clear();
1701  for (int i = 0; i < result.size(); i += 3)
1702  {
1703  capping_contour.ring.push_back(result[i + 0]);
1704  capping_contour.ring.push_back(result[i + 1]);
1705  capping_contour.ring.push_back(result[i + 2]);
1706  capping_contour.ring_index.push_back(i);
1707  }
1708  }
1709 
1710  } // intsec_vector.size() > 3
1711 
1712  // Data output
1713  for (vector<PLANE_EDGE_INTERSECTION>::iterator intsec = intsec_vector.begin(); intsec < intsec_vector.end(); ++intsec)
1714  {
1715  u_data_vector.push_back(intsec->data_vertex_int.v[0]);
1716  if (intsec->data_vertex_int.dimension == 3)
1717  {
1718  v_data_vector.push_back(intsec->data_vertex_int.v[1]);
1719  w_data_vector.push_back(intsec->data_vertex_int.v[2]);
1720  }
1721  }
1722 }
1723 }
GLenum GLsizei len
Definition: khronos-glext.h:7440
__host__ __device__ float2 fabs(float2 v)
Definition: cutil_math.h:1380
vector< int > ring_index
Definition: CuttingSurfaceGPMUtil.h:40
double x
Definition: CuttingSurfaceGPMUtil.h:21
std::vector< PLANE_EDGE_INTERSECTION > PLANE_EDGE_INTERSECTION_VECTOR
Definition: CuttingSurfaceGPMUtil.h:51
void generate_capping_contour(CONTOUR &capping_contour, PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, EDGE_VECTOR &plane_base_x, EDGE_VECTOR &plane_base_y, vector< float > &u_data_vector, vector< float > &v_data_vector, vector< float > &w_data_vector)
Definition: CuttingSurfaceGPMUtil.h:1590
GLfloat GLfloat GLfloat v2
Definition: khronos-glext.h:6754
Definition: CuttingSurfaceGPMUtil.h:37
vector< int > polyhedron_faces
Definition: CuttingSurfaceGPMUtil.h:41
PLANE_EDGE_INTERSECTION PlaneEdgeVertexInterpolate(float x1, float y1, float z1, float x2, float y2, float z2, S_V_DATA data1, S_V_DATA data2, int v1, int v2, double p, EDGE_VECTOR unit_normal_vector, double D1, double D2)
Definition: CuttingSurfaceGPMUtil.h:75
double y
Definition: CuttingSurfaceGPMUtil.h:22
EDGE_VECTOR intersection
Definition: CuttingSurfaceGPMUtil.h:35
int vertex2
Definition: CuttingSurfaceGPMUtil.h:32
GLfixed GLfixed GLfixed y2
Definition: khronos-glext.h:11325
#define NULL
Definition: covise_list.h:22
bool intersection_at_vertex2
Definition: CuttingSurfaceGPMUtil.h:28
Definition: Triangulate.h:12
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
double z
Definition: CuttingSurfaceGPMUtil.h:23
int vertex1
Definition: CuttingSurfaceGPMUtil.h:31
bool test_plane_edge_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, PLANE_EDGE_INTERSECTION &intsec, float *x_coord_in, float *y_coord_in, float *z_coord_in)
Definition: CuttingSurfaceGPMUtil.h:221
GLuint GLuint end
Definition: khronos-glext.h:6343
float v[3]
Definition: CuttingSurfaceGPMUtil.h:15
S_V_DATA data_vertex_int
Definition: CuttingSurfaceGPMUtil.h:34
GLsizeiptr size
Definition: khronos-glext.h:6610
bool find_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
Definition: CuttingSurfaceGPMUtil.h:512
float y
Definition: Triangulate.h:15
const GLdouble * v
Definition: khronos-glext.h:6442
GLfloat GLfloat p
Definition: khronos-glext.h:9861
int int_flag
Definition: CuttingSurfaceGPMUtil.h:30
std::vector< tr_vertex > tr_vertexVector
Definition: Triangulate.h:19
EDGE_VECTOR cross_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:59
double length(EDGE_VECTOR &vector)
Definition: CuttingSurfaceGPMUtil.h:70
double dot_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:54
GLuint index
Definition: khronos-glext.h:6722
GLfixed GLfixed x2
Definition: khronos-glext.h:11325
GLuint64EXT * result
Definition: khronos-glext.h:12573
Definition: CuttingSurfaceGPMUtil.h:13
std::vector< int > tr_intVector
Definition: Triangulate.h:21
GLuint GLfloat GLfloat GLfloat x1
Definition: khronos-glext.h:13144
static bool Process(const tr_vertexVector &contour, tr_intVector &result)
Definition: Triangulate.h:28
int dimension
Definition: CuttingSurfaceGPMUtil.h:16
int assign_int_index(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
Definition: CuttingSurfaceGPMUtil.h:469
double get_cosine(EDGE_VECTOR &vector, EDGE_VECTOR &plane_base_x, EDGE_VECTOR &plane_base_y)
Definition: CuttingSurfaceGPMUtil.h:1569
int index
Definition: Triangulate.h:16
bool intersection_at_vertex1
Definition: CuttingSurfaceGPMUtil.h:27
GLdouble GLdouble t
Definition: khronos-glext.h:6449
Definition: CuttingSurfaceGPMUtil.h:19
GLfloat GLfloat v1
Definition: khronos-glext.h:6753
float x
Definition: Triangulate.h:14
Definition: CuttingSurfaceGPMUtil.h:25
vector< int > ring
Definition: CuttingSurfaceGPMUtil.h:39
GLfixed y1
Definition: khronos-glext.h:11325