COVISE Core
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
10namespace covise
11{
12
13typedef struct
14{
15 float v[3];
16 int dimension;
17} S_V_DATA;
18
19typedef struct
20{
21 double x;
22 double y;
23 double z;
24} EDGE_VECTOR;
25typedef struct
26{
27 bool intersection_at_vertex1;
28 bool intersection_at_vertex2;
29
30 int int_flag;
31 int vertex1;
32 int vertex2;
33
34 S_V_DATA data_vertex_int;
35 EDGE_VECTOR intersection;
36} PLANE_EDGE_INTERSECTION;
37typedef 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
51typedef std::vector<PLANE_EDGE_INTERSECTION> PLANE_EDGE_INTERSECTION_VECTOR;
52//typedef std::vector<PLANE_CELL_INTERSECTION> PLANE_CELL_INTERSECTION_VECTOR;
53
54double 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
59EDGE_VECTOR cross_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
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
70double length(EDGE_VECTOR &vector)
71{
72 return sqrt(vector.x * vector.x + vector.y * vector.y + vector.z * vector.z);
73}
74
75PLANE_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
221bool 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
300PLANE_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
469int 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
512bool 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
1569double 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
1590void 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 }
1697 tr_intVector result;
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}
NULL
#define NULL
Definition: covise_list.h:22
y2
GLfixed GLfixed GLfixed y2
Definition: khronos-glext.h:11325
result
GLuint64EXT * result
Definition: khronos-glext.h:12573
y1
GLfixed y1
Definition: khronos-glext.h:11325
size
GLsizeiptr size
Definition: khronos-glext.h:6610
end
GLuint GLuint end
Definition: khronos-glext.h:6343
v
const GLdouble * v
Definition: khronos-glext.h:6442
x1
GLuint GLfloat GLfloat GLfloat x1
Definition: khronos-glext.h:13144
index
GLuint index
Definition: khronos-glext.h:6722
t
GLdouble GLdouble t
Definition: khronos-glext.h:6449
x2
GLfixed GLfixed x2
Definition: khronos-glext.h:11325
p
GLfloat GLfloat p
Definition: khronos-glext.h:9861
len
GLenum GLsizei len
Definition: khronos-glext.h:7440
v1
GLfloat GLfloat v1
Definition: khronos-glext.h:6753
v2
GLfloat GLfloat GLfloat v2
Definition: khronos-glext.h:6754
tr_intVector
std::vector< int > tr_intVector
Definition: Triangulate.h:24
tr_vertexVector
std::vector< tr_vertex > tr_vertexVector
Definition: Triangulate.h:22
fabs
__host__ __device__ float2 fabs(float2 v)
Definition: cutil_math.h:1380
covise
list of all chemical elements
Definition: coConfig.h:27
covise::cross_product
EDGE_VECTOR cross_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:59
covise::PLANE_EDGE_INTERSECTION_VECTOR
std::vector< PLANE_EDGE_INTERSECTION > PLANE_EDGE_INTERSECTION_VECTOR
Definition: CuttingSurfaceGPMUtil.h:51
covise::test_plane_edge_intersection
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
covise::generate_capping_contour
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
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::PlaneEdgeVertexInterpolate
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
covise::length
double length(EDGE_VECTOR &vector)
Definition: CuttingSurfaceGPMUtil.h:70
covise::get_cosine
double get_cosine(EDGE_VECTOR &vector, EDGE_VECTOR &plane_base_x, EDGE_VECTOR &plane_base_y)
Definition: CuttingSurfaceGPMUtil.h:1569
covise::dot_product
double dot_product(EDGE_VECTOR &vector1, EDGE_VECTOR &vector2)
Definition: CuttingSurfaceGPMUtil.h:54
covise::find_intersection
bool find_intersection(PLANE_EDGE_INTERSECTION_VECTOR &intsec_vector, int edge_vertex1, int edge_vertex2)
Definition: CuttingSurfaceGPMUtil.h:512
covise::calculate_intersections
PLANE_EDGE_INTERSECTION_VECTOR calculate_intersections(float *udata_in, float *vdata_in, float *wdata_in, int num_elem_in, int *elem_in, int num_conn_in, int *conn_in, float *x_coord_in, float *y_coord_in, float *z_coord_in, float p, EDGE_VECTOR unit_normal_vector)
Definition: CuttingSurfaceGPMUtil.h:300
tr_vertex
Definition: Triangulate.h:16
tr_vertex::x
float x
Definition: Triangulate.h:17
tr_vertex::y
float y
Definition: Triangulate.h:18
tr_vertex::index
int index
Definition: Triangulate.h:19
Triangulate::Process
static bool Process(const tr_vertexVector &contour, tr_intVector &result)
Definition: Triangulate.h:31
covise::S_V_DATA
Definition: CuttingSurfaceGPMUtil.h:14
covise::S_V_DATA::dimension
int dimension
Definition: CuttingSurfaceGPMUtil.h:16
covise::S_V_DATA::v
float v[3]
Definition: CuttingSurfaceGPMUtil.h:15
covise::EDGE_VECTOR
Definition: CuttingSurfaceGPMUtil.h:20
covise::EDGE_VECTOR::x
double x
Definition: CuttingSurfaceGPMUtil.h:21
covise::EDGE_VECTOR::z
double z
Definition: CuttingSurfaceGPMUtil.h:23
covise::EDGE_VECTOR::y
double y
Definition: CuttingSurfaceGPMUtil.h:22
covise::PLANE_EDGE_INTERSECTION
Definition: CuttingSurfaceGPMUtil.h:26
covise::PLANE_EDGE_INTERSECTION::intersection_at_vertex1
bool intersection_at_vertex1
Definition: CuttingSurfaceGPMUtil.h:27
covise::PLANE_EDGE_INTERSECTION::intersection_at_vertex2
bool intersection_at_vertex2
Definition: CuttingSurfaceGPMUtil.h:28
covise::PLANE_EDGE_INTERSECTION::vertex1
int vertex1
Definition: CuttingSurfaceGPMUtil.h:31
covise::PLANE_EDGE_INTERSECTION::vertex2
int vertex2
Definition: CuttingSurfaceGPMUtil.h:32
covise::PLANE_EDGE_INTERSECTION::intersection
EDGE_VECTOR intersection
Definition: CuttingSurfaceGPMUtil.h:35
covise::PLANE_EDGE_INTERSECTION::int_flag
int int_flag
Definition: CuttingSurfaceGPMUtil.h:30
covise::PLANE_EDGE_INTERSECTION::data_vertex_int
S_V_DATA data_vertex_int
Definition: CuttingSurfaceGPMUtil.h:34
covise::CONTOUR
Definition: CuttingSurfaceGPMUtil.h:38
covise::CONTOUR::ring_index
vector< int > ring_index
Definition: CuttingSurfaceGPMUtil.h:40
covise::CONTOUR::polyhedron_faces
vector< int > polyhedron_faces
Definition: CuttingSurfaceGPMUtil.h:41
covise::CONTOUR::ring
vector< int > ring
Definition: CuttingSurfaceGPMUtil.h:39
Generated by doxygen 1.9.4